[gmx-users] LINCS warning in md run
reisingere at rostlab.informatik.tu-muenchen.de
reisingere at rostlab.informatik.tu-muenchen.de
Thu Sep 20 14:15:09 CEST 2012
Hi Justin,
thank you a lot for your answer. I will try it.
Best,
Eva
>
>
> 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
>
> ========================================
> --
> 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
>
More information about the gromacs.org_gmx-users
mailing list