On 9/20/12 6:52 AM, reising...@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, reising...@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@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-requ...@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@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-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
