[gmx-users] LINCS warning in md run
Justin Lemkul
jalemkul at vt.edu
Thu Sep 20 14:04:39 CEST 2012
On 9/20/12 6:52 AM, reisingere at rostlab.informatik.tu-muenchen.de wrote:
> Hi,
> now I tried it without any restriction and still the LINCS warnings occur.
> Since it is always the hydrogen atom where the huge force lies on I had
> the idea just to remove the hydrogen atom to find out whether another atom
> will occur to have such a high force on or everything went fine. And
> indeed, without the hydrogen atom there is no problem during the MD run.
> Even when I fix all the other atoms.
>
> I can not see why this occurs. I already looked at the residue with pymol
> and I could see that there is no other residue or atom clashing with it.
>
I suspect what is happening is that the topology is not stable. If removal of
the H atom leads to stable trajectories, then it is experiencing some
interaction that leads to the crash. Likely what is happening is that it is
experiencing a very strong 1-4 interaction with the other phosphate oxygens and
is being pulled very hard away from the oxygen to which it is bonded, leading to
deformed geometry and the LINCS warnings. There are two possible solutions:
1. Modify the 1-4 interactions through an appropriate pair term.
2. Add exclusions to the topology between the H atom and each of the O atoms to
which it is not bonded.
I would recommend testing with a very simple system of TYP in water to validate
that the topology is stable and rule out other factors.
-Justin
> The residue looks like this:
>
> ATOM 955 N TYP 65 27.310 23.150 27.670 1.00 0.00
> N
> ATOM 956 CA TYP 65 26.870 23.320 26.300 1.00 0.00
> C
> ATOM 957 C TYP 65 27.920 23.760 25.300 1.00 0.00
> C
> ATOM 958 O TYP 65 28.060 23.210 24.210 1.00 0.00
> O
> ATOM 959 CB TYP 65 25.740 24.340 26.170 1.00 0.00
> C
> ATOM 960 CG TYP 65 24.530 23.740 26.830 1.00 0.00
> C
> ATOM 961 CD1 TYP 65 24.410 23.760 28.220 1.00 0.00
> C
> ATOM 962 CD2 TYP 65 23.510 23.180 26.050 1.00 0.00
> C
> ATOM 963 CE1 TYP 65 23.290 23.210 28.840 1.00 0.00
> C
> ATOM 964 CE2 TYP 65 22.380 22.630 26.660 1.00 0.00
> C
> ATOM 965 CZ TYP 65 22.260 22.630 28.070 1.00 0.00
> C
> ATOM 966 OH TYP 65 21.160 22.100 28.670 1.00 0.00
> O
> ATOM 967 P TYP 65 20.298 21.422 27.722 1.00 0.00
> P
> ATOM 968 OP1 TYP 65 21.003 21.245 26.467 1.00 0.00
> O
> ATOM 969 OP2 TYP 65 19.920 20.126 28.251 1.00 0.00
> O
> ATOM 970 OP3 TYP 65 19.106 22.218 27.498 1.00 0.00
> O
> ATOM 971 H TYP 65 27.210 23.914 28.325 1.00 0.00
> H
> ATOM 972 HA TYP 65 26.485 22.366 25.933 1.00 0.00
> H
> ATOM 973 HB1 TYP 65 25.524 24.543 25.119 1.00 0.00
> H
> ATOM 974 HB2 TYP 65 26.010 25.270 26.669 1.00 0.00
> H
> ATOM 975 HD1 TYP 65 25.154 24.234 28.821 1.00 0.00
> H
> ATOM 976 HD2 TYP 65 23.603 23.156 24.973 1.00 0.00
> H
> ATOM 977 HE1 TYP 65 23.200 23.222 29.916 1.00 0.00
> H
> ATOM 978 HE2 TYP 65 21.603 22.189 26.053 1.00 0.00
> H
> ATOM 979 H1P TYP 65 20.357 20.996 25.802 1.00 0.00
> H
>
>
> The topology entry in the aminoacid.rtp file looks like this:
>
> [ TYP ]
> [ atoms ]
> N N -0.516300 1
> CA CT 0.275503 2
> HA H1 0.008223 3
> CB CT -0.354052 4
> HB1 HC 0.110326 5
> HB2 HC 0.110326 6
> CG CA 0.119728 7
> CD1 CA -0.198938 8
> HD1 HA 0.137143 9
> CE1 CA -0.284884 10
> HE1 HA 0.177179 11
> CZ C 0.452616 12
> OH OS -0.534452 13
> H H 0.293600 14
> CE2 CA -0.284884 15
> HE2 HA 0.177179 16
> CD2 CA -0.198938 17
> HD2 HA 0.137143 18
> C C 0.536600 19
> O O -0.581900 20
> P P 1.393213 21
> OP1 OH -0.752821 22
> OP2 O2 -0.822464 23
> OP3 O2 -0.822464 24
> H1P HO 0.423316 25
>
>
> [ bonds ]
> N H
> N CA
> CA HA
> CA CB
> CA C
> CB HB1
> CB HB2
> CB CG
> CG CD1
> CG CD2
> CD1 HD1
> CD1 CE1
> CE1 HE1
> CE1 CZ
> CZ OS
> CZ CE2
> OS P
> CE2 HE2
> CE2 CD2
> CD2 HD2
> C O
> -C N
> P OP1
> P OP2
> P OP3
> H1P OP1
>
>
> [ impropers ]
> -C CA N H
> CA +N C O
> CG CE2 CD2 HD2
> CZ CD2 CE2 HE2
> CD1 CZ CE1 HE1
> CG CE1 CD1 HD1
> CD1 CD2 CG CB
> CE1 CE2 CZ OH
>
>
> The LINCS problems lie between the oxygen atom of the phosphate where a
> hydrogen atom is bound to (OP1) and the hydrogen atom bound to the oxygen
> atom of the phosphate (H1P).
>
> Do you see where the problem lies?
> Thank you,
> Eva
>
>>
>>
>> On 9/18/12 8:58 AM, reisingere at rostlab.informatik.tu-muenchen.de wrote:
>>> Okey,
>>> now I tried it without any fixed residues. But still the energy after
>>> the
>>> minimization is not very low and I still get the LINCS warnings.
>>>
>>> The mdp file I use for the minimization looks like this:
>>>
>>> define = -DPOSRES
>>
>> Restraints during minimization generally restrict motion in the same way
>> that
>> freezing does. Again, this is a potential barrier to sufficient
>> minimization.
>>
>>> integrator = steep
>>> emtol = 10
>>> nsteps = 1500
>>> nstenergy = 1
>>> energygrps = System
>>> coulombtype = PME
>>> rcoulomb = 0.9
>>> rvdw = 0.9
>>> rlist = 0.9
>>> fourierspacing = 0.12
>>> pme_order = 4
>>> ewald_rtol = 1e-5
>>> pbc = xyz
>>>
>>>
>>> the mdp file for the md run looks like this:
>>>
>>> define = -DPOSRES
>>> integrator = md
>>> dt = 0.001
>>> nsteps = 5000
>>> nstxout = 100
>>> nstvout = 0
>>> nstfout = 0
>>> nstlog = 1000
>>> nstxtcout = 500
>>> nstenergy = 5
>>> energygrps = Protein Non-Protein
>>> nstcalcenergy = 5
>>> nstlist = 10
>>> ns-type = Grid
>>> pbc = xyz
>>> rlist = 0.9
>>> coulombtype = PME
>>> rcoulomb = 0.9
>>> rvdw = 0.9
>>> fourierspacing = 0.12
>>> pme_order = 4
>>> ewald_rtol = 1e-5
>>> gen_vel = yes
>>> gen_temp = 200.0
>>> gen_seed = 9999
>>> constraints = all-bonds
>>> tcoupl = V-rescale
>>> tc-grps = Protein Non-Protein
>>> tau_t = 0.1 0.1
>>> ref_t = 298 298
>>> pcoupl = no
>>>
>>>
>>>
>>> The output of the minimization run is:
>>>
>>> Steepest Descents converged to machine precision in 15 steps,
>>> but did not reach the requested Fmax < 10.
>>> Potential Energy = -7.9280412e+05
>>> Maximum force = 1.0772942e+04 on atom 979
>>> Norm of force = 9.6685356e+01
>>>
>>
>> Note that this outcome is even worse than before. The maximum force is
>> now over
>> 10,000.
>>
>>>
>>>
>>>
>>> The output of the MD run is:
>>>
>>> Step 1031, time 1.031 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 0.001431, max 0.010707 (between atoms 976 and 979)
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>> 976 979 81.6 0.1272 0.0970 0.0960
>>>
>>>
>>
>> As you should expect. If you can't get a reasonable minimization, the MD
>> will
>> always crash.
>>
>>>
>>> The atom 979 is the hydrogen atom on the phosphate. There has to be one
>>> hydrogen atoms because it is protonated once. The other atom 976 is the
>>> oxygen atom where the hydrogen atom is bound to.
>>> The bounding parameters for this kind of binding were already there. I
>>> didn't add them.
>>>
>>> I already did it for another phosphorylation on another position in this
>>> structure. And here I also got many LINCS errors. And again the problem
>>> is
>>> the connection between the hydrogen atom and the oxygen atom.
>>>
>>> But I do not understand why.
>>
>> Remove all restraints/freezing/whatever in the .mdp file and try the EM
>> again.
>> If it still does not converge to a reasonable value, then there are two
>> potential problems:
>>
>> 1. The topology is flawed and thus the structure cannot be run stably
>> 2. The structure cannot accommodate phosphate in this location and due to
>> unresolvable clashes, the runs fail
>>
>> -Justin
>>
>> --
>> ========================================
>>
>> Justin A. Lemkul, Ph.D.
>> Research Scientist
>> Department of Biochemistry
>> Virginia Tech
>> Blacksburg, VA
>> jalemkul[at]vt.edu | (540) 231-9080
>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>>
>> ========================================
>> --
>> gmx-users mailing list gmx-users at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-request at gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>
>
--
========================================
Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
More information about the gromacs.org_gmx-users
mailing list