Re: Fwd: Atoms too fast/periodic cell too small with ABF protein-ligand

From: Francesco Pietra (chiendarret_at_gmail.com)
Date: Thu Jul 06 2017 - 01:15:50 CDT

Hi Giacomo:

adding "wrapAll no" after running with "wrapAll yes" will not unwrap the
> coordinates that were wrapped in the previous simulations.
>
If potentially useful, I could run an MD from scratch with "wrapAll no". I
have access again to the cluster and I imagine that a few ns would be
enough (previous 58.2 ns were for other reasons than ABF/FEP, as a
preliminary to an accelerated MD).

that's mostly a concern for a continuous over time definition of your
> centers of mass.
>
I fear I do not understand this fully

can you produce a minimal input that reproduces the problem and include it
> as a download link?
>
These preliminary assays are for a large ABF/FEP project turning around
this protein and ligands for which I received yesterday a grant. So, I
can't put it on the web, but I could send file(s) to you personally.

Actually, I am now somewhat in a trouble, because, after having
successfully parameterized the organic ligands, I took for granted that
ABF/FEP will run. On the other hand, I recognized later that the problem
can't be attacked with any current form of accelerated MD, because MD would
be unable to reposition the ligand. I need free energies for various poses
of the ligand.

Since 6A (as obtained for "r" by "cv printframe") is unrealistic, while 29A
(as measured directly from the centers of mass) is acceptable froma rough
examination of the model, I was wondering whether "cv printframe" is
measuring the ligand against the protein chain where it resides, not agains
the whole protein (it is a two chains proteins, and there are two other
ligands, very common ones in biochemical systems, therefore with a good
parameterization, and on which ABF/FEP will have nothing to do). If so, it
might well be that also Theta, Psi, etc are measured incorrectly by "cv
printframe".

I could use the sole chain where the ligand of interest resides, but that
would be an oversimplification, removing the system from reality too much.

thanks

francesco

On Wed, Jul 5, 2017 at 6:30 PM, Giacomo Fiorin <giacomo.fiorin_at_gmail.com>
wrote:

> Hi Francesco, adding "wrapAll no" after running with "wrapAll yes" will
> not unwrap the coordinates that were wrapped in the previous simulations.
> But that's mostly a concern for a continuous over time definition of your
> centers of mass.
>
> As for the consistency, can you produce a minimal input that reproduces
> the problem and include it as a download link?
>
> Giacomo
>
> On Wed, Jul 5, 2017 at 12:03 PM, Francesco Pietra <chiendarret_at_gmail.com>
> wrote:
>
>> Hi Giacomo:
>> I have now continued the "wrapall on" 58.2ns MD for further 0.11ns under
>> "wrapall no". The system remains in order, except a few TIP3Ps out of the
>> box margins.
>> All values from "cv printframe", are substantially unchanged, except Psi,
>> which has changed sign, from +166.2 to -168.1. Also the directly measured
>> "r" is substantially unchanged.
>>
>> "r" changed from 5.97/27.6 to 5.92/27.9, where the larger value is from
>> direct measurements as follows:
>>
>> set selprot [atomselect top "protein"]
>> measure center $selprot
>> set sellig [atomselect top "segname SAA1"]
>> measure center $sellig
>> set dist [ veclength [vecsub [measure center $selprot] [measure center
>> $sellig]]]
>> 27.91017455323315
>>
>> I loaded to VMD the DCSs from NAMD, I suppose that PBCs are the same.
>> Therefore, I do not understand why "r" is so much different with the two
>> methods.
>>
>> If - as it is likely - there is a bug in what I did, and such a bug can
>> be felt, I would like to learn how to get out of this labyrinth.
>>
>> thanks a lot
>>
>> francesco
>>
>>
>> On Wed, Jun 28, 2017 at 6:35 PM, Giacomo Fiorin <giacomo.fiorin_at_gmail.com
>> > wrote:
>>
>>> Hi Francesco, wrapAll is not recommended for a system where more than
>>> one molecule is involved in the variables' definition:
>>> http://colvars.github.io/colvars-refman-namd/colvars-refman-
>>> namd.html#x1-160004.3
>>>
>>> Again, the computation of a center of mass for a group requires
>>> unwrapped coordinates. I would double check that the same exact
>>> configuration is being loaded into VMD, against a NAMD simulation with 0
>>> steps.
>>>
>>> Giacomo
>>>
>>>
>>>
>>> On Tue, Jun 27, 2017 at 3:46 AM, Francesco Pietra <chiendarret_at_gmail.com
>>> > wrote:
>>>
>>>> Hi Giacomo:
>>>> I am using for ABF the last frame (.xsc file) of MD equilibration
>>>> generated without PBCs enabled in NAMD config file. That config file
>>>> only had:
>>>>
>>>> wrapNearest no
>>>> wrapAll on
>>>>
>>>> If I measure directly the distance between the centers of mass of
>>>> protein and ligand for that frame:
>>>>
>>>>
>>>>> >Main< (francesco) 33 % set selprot [atomselect top "protein"]
>>>>> atomselect3
>>>>> >Main< (francesco) 34 % measure center $selprot
>>>>> 326.3038024902344 453.4432678222656 342.5987243652344
>>>>> >Main< (francesco) 35 % set sellig [atomselect top "segname SAA1"]
>>>>> atomselect4
>>>>> >Main< (francesco) 36 % measure center $sellig
>>>>> 329.97216796875 441.479736328125 317.9831237792969
>>>>>
>>>> >Main< (francesco) 37% set dist [veclength [vecsub [measure center
>>>> $selprot] [measure center $sellig]]]
>>>> 27.613597797130065
>>>> >Main< (francesco) %
>>>>
>>>>
>>>> I get a result (27.6) different from "cv printframe" by loading that
>>>> NAMD frame to VMD (18.9).
>>>>
>>>> Should any of these conditions be unsuitable for ABF, then, without a
>>>> complete reeducation, it will be difficult for me carrying out an ABF
>>>> for protein-ligand. However, I am strongly interested in that particular
>>>> protein-ligand system, I spend a lot of time in parameterizing the ligand.
>>>> I imagine that with FEP (not yet attempted) I'll be faced by samilar
>>>> problems.
>>>>
>>>> francesco
>>>>
>>>>
>>>>
>>>> On Mon, Jun 26, 2017 at 6:32 PM, Giacomo Fiorin <
>>>> giacomo.fiorin_at_gmail.com> wrote:
>>>>
>>>>> Hi Francesco, when distances between groups of atoms are involved, it
>>>>> is really important to check that the PBCs are the same in both VMD and
>>>>> NAMD.
>>>>>
>>>>> Giacomo
>>>>>
>>>>> On Mon, Jun 26, 2017 at 12:24 PM, Francesco Pietra <
>>>>> chiendarret_at_gmail.com> wrote:
>>>>>
>>>>>> I also run the same ABF (rmsd all atoms of the ligand less methyl
>>>>>> hydrogens) with colvar "r" = 18.9 from cv printframe (unlike the directly
>>>>>> measured 27.6 with previous ABFs). In all cases, I took the colvars values
>>>>>> from the last frame of the 58.2ns MD, not averages along the whole
>>>>>> simulation)
>>>>>>
>>>>>>
>>>>>> Main< (Conformation:5xxf-SAA1) 28 % cv printframe
>>>>>>> 0 1.88926419197828e+01 5.22845307973583e+01
>>>>>>> -3.39612142392650e+01 9.26730407864645e+00 1.01135055762437e+02
>>>>>>> -5.05311208874754e+01 0.00000000000000e+00 0.00000000000000e+00
>>>>>>>
>>>>>>> >Main< (Conformation:5xxf-SAA1) 29 %
>>>>>>>
>>>>>>> >Main< (Conformation:Bound_5xxf-SAA1-allH_except_methyl_H) 33 % cv
>>>>>>> printframelabels
>>>>>>> # step r Theta
>>>>>>> Phi Psi theta
>>>>>>> phi RMSD r_RMSD
>>>>>>>
>>>>>>
>>>>>> The immediate (step 0) error was again atoms moving too fast, this
>>>>>> time one heavy atom and one H-atom of the ligand.
>>>>>>
>>>>>> f
>>>>>>
>>>>>> ---------- Forwarded message ----------
>>>>>> From: Francesco Pietra <chiendarret_at_gmail.com>
>>>>>> Date: Mon, Jun 26, 2017 at 12:28 PM
>>>>>> Subject: Fwd: Atoms too fast/periodic cell too small with ABF
>>>>>> protein-ligand
>>>>>> To: NAMD <namd-l_at_ks.uiuc.edu>
>>>>>>
>>>>>>
>>>>>> I must add that the attempted ABF was "Conformation:Bound". For rmsd
>>>>>> for the ligand I had chosen all heavy atoms.
>>>>>> By choosing all atoms of the ligand, except hydrogens at the methyl
>>>>>> groups, the simulation also crashed at the first step, this time for two H
>>>>>> atoms of the protein, one (HG1) very far from the ligand, the other one
>>>>>> (HA) close to the ligand.
>>>>>>
>>>>>> With "Conformation:Unbound" the ABF went to completion without errors.
>>>>>>
>>>>>> Hope this helps suggesting what to do.
>>>>>>
>>>>>> As I said, there was no problem with MD equilibration for this
>>>>>> system, at the same ts=1.0fs, along a trajectory of during 58.2ns.
>>>>>>
>>>>>> francesco
>>>>>>
>>>>>>
>>>>>> ---------- Forwarded message ----------
>>>>>> From: Francesco Pietra <chiendarret_at_gmail.com>
>>>>>> Date: Sun, Jun 25, 2017 at 12:56 PM
>>>>>> Subject: Atoms too fast/periodic cell too small with ABF
>>>>>> protein-ligand
>>>>>> To: NAMD <namd-l_at_ks.uiuc.edu>
>>>>>>
>>>>>>
>>>>>> Hello:
>>>>>>
>>>>>> I am attempting a protein-ligand ABF, following the 2017 tutorial, by
>>>>>> using a 58.2ns problemless equilibrated system (ts=1.0 fs, bond restriction
>>>>>> on water only) in a TIP3P water box on a main pure-CPU cluster. Rather
>>>>>> large protein, organic ligand as accurately parameterized as I could, by
>>>>>> fitting torsions and water interaction.
>>>>>>
>>>>>> I am experiencing immediate "atoms moving too fast" (two H atoms of
>>>>>> the ligand) when using a linux 4-core cpu desktop, or "periodic cell has
>>>>>> become too small" on a linux GPU workstation. At the moment I have no
>>>>>> access to the cluster.
>>>>>>
>>>>>> I used ts=1.0 fs, i.e. no bond restriction, except for TIP3P water,
>>>>>> as in the equilibration.
>>>>>>
>>>>>> As a possible cause, that I was unable to verify, is the setting of
>>>>>> Euler and polar angle in the colvars definition. That is , I used the
>>>>>> following values from "cv printframe"
>>>>>>
>>>>>> >Main< (Conformation:5syf-SAA1) 28 % cv printframe
>>>>>>> 0 1.88926419197828e+01 5.22845307973583e+01
>>>>>>> -3.39612142392650e+01 9.26730407864645e+00 1.01135055762437e+02
>>>>>>> -5.05311208874754e+01 0.00000000000000e+00 0.00000000000000e+00
>>>>>>>
>>>>>>
>>>>>> by discarding the first two, and using the directly measured
>>>>>> intercenter distance of 27.6A in place of 18.89 from printframe, i.e., as
>>>>>> follows:
>>>>>>
>>>>>>
>>>>>> harmonic {
>>>>>> colvars r
>>>>>> forceConstant 0.0
>>>>>> centers 27.6 # OK measured
>>>>>> }
>>>>>>
>>>>>>
>>>>>> harmonic {
>>>>>> colvars Theta
>>>>>> forceConstant 0.0
>>>>>> centers 52.3 # from printframe
>>>>>> }
>>>>>>
>>>>>>
>>>>>> harmonic {
>>>>>> colvars Phi
>>>>>> forceConstant 0.0
>>>>>> centers -40.0 # from printframe
>>>>>> }
>>>>>>
>>>>>>
>>>>>> harmonic {
>>>>>> colvars Psi
>>>>>> forceConstant 0.0
>>>>>> centers 9.27 # from printframe
>>>>>> }
>>>>>>
>>>>>>
>>>>>> harmonic {
>>>>>> colvars theta
>>>>>> forceConstant 0.0
>>>>>> centers 101.1 # from printframe
>>>>>> }
>>>>>>
>>>>>>
>>>>>> harmonic {
>>>>>> colvars phi
>>>>>> forceConstant 0.0
>>>>>> centers -50.5 # from printframe
>>>>>> }
>>>>>>
>>>>>> Should this assignment of colvars be correct, where to look for the
>>>>>> causes of the instability of the system?
>>>>>> I must confess that I am the first time at such an ABF beyond the
>>>>>> tutorial
>>>>>>
>>>>>> Thanks for advice
>>>>>>
>>>>>> francesco pietra
>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> Giacomo Fiorin
>>>>> Associate Professor of Research, Temple University, Philadelphia, PA
>>>>> Contractor, National Institutes of Health, Bethesda, MD
>>>>> http://goo.gl/Q3TBQU
>>>>> https://github.com/giacomofiorin
>>>>>
>>>>
>>>>
>>>
>>>
>>> --
>>> Giacomo Fiorin
>>> Associate Professor of Research, Temple University, Philadelphia, PA
>>> Contractor, National Institutes of Health, Bethesda, MD
>>> http://goo.gl/Q3TBQU
>>> https://github.com/giacomofiorin
>>>
>>
>>
>
>
> --
> Giacomo Fiorin
> Associate Professor of Research, Temple University, Philadelphia, PA
> Contractor, National Institutes of Health, Bethesda, MD
> http://goo.gl/Q3TBQU
> https://github.com/giacomofiorin
>

This archive was generated by hypermail 2.1.6 : Mon Dec 31 2018 - 23:20:25 CST