[gmx-users] Exploding poly-serine

[Justin A. Lemkul]

Warren Gallin wrote:
> Hi,
> 
>     I am trying to work my way through learning to use GROMACS 4.0.5 for 
> doing MD simulations of short peptides in solution.
> 
>     My problem is that several peptides, most strikingly a Serine 
> 10-mer, are exploding during the production run.
> 
>     I construct the serine polymer in extended form using TINKER, with 
> neither end capped, then use pdb2gmx, editconv, genbox and genion to set 
> the forcefield to OPLSAA/L, place it in a box of tip3p water, 
> neutralized in a .096 M NaCl solution.
> 
>     I am trying to track down why this is occurring, and the first thing 
> that seems suspicious is that the initial energy minimization run 
> reaches convergence with a rather large force remaining on one of the 
> atoms - the final statement from the EM run is:
> 
> 
> Steepest Descents converged to machine precision in 6793 steps,
> but did not reach the requested Fmax < 10.
> Potential Energy  = -4.9731194e+05
> Maximum force     =  7.1219537e+02 on atom 111
> Norm of force     =  6.5289087e+00
> 
> 
>     So I am thinking that with a large initial force on the polymer the 
> system might be unrecoverably unstable and this is propagating through 
> the subsequent steps of relaxing the water and the actual MD run to pop 
> up as an explosion.
> 
>     If I look at the run log for the MD run, everything seems stable 
> until the last step, at which point the temperature shoots up and the 
> serine polymer explodes.  Here are the last four steps from the run log:
> 
>            Step           Time         Lambda
>          208200      416.40002        0.00000
> 
>    Energies (kJ/mol)
>            Bond          Angle    Proper Dih. Ryckaert-Bell.          LJ-14
>     9.75158e+03    1.25586e+02    7.62176e+00   -6.45482e+01    1.16773e+02
>      Coulomb-14        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential
>     1.80001e+03    6.08872e+04   -4.19574e+05   -4.65971e+04   -3.93547e+05
>     Kinetic En.   Total Energy  Conserved En.    Temperature Pressure (bar)
>     7.77043e+04   -3.15843e+05    1.01418e+06    3.17587e+02   -2.87690e+02
> 
>            Step           Time         Lambda
>          208300      416.60002        0.00000
> 
>    Energies (kJ/mol)
>            Bond          Angle    Proper Dih. Ryckaert-Bell.          LJ-14
>     1.04345e+04    1.96147e+02    1.55806e+01   -3.31126e+01    1.24787e+02
>      Coulomb-14        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential
>     1.80413e+03    6.08422e+04   -4.19233e+05   -4.66893e+04   -3.92538e+05
>     Kinetic En.   Total Energy  Conserved En.    Temperature Pressure (bar)
>     7.77229e+04   -3.14815e+05    1.02114e+06    3.17663e+02   -2.94492e+02
> 
>            Step           Time         Lambda
>          208400      416.80002        0.00000
> 
>    Energies (kJ/mol)
>            Bond          Angle    Proper Dih. Ryckaert-Bell.          LJ-14
>     7.10106e+03    2.06975e+02    1.21357e+01   -2.52705e+01    1.16229e+02
>      Coulomb-14        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential
>     1.85247e+03    6.04882e+04   -4.18865e+05   -4.67103e+04   -3.95824e+05
>     Kinetic En.   Total Energy  Conserved En.    Temperature Pressure (bar)
>     7.58806e+04   -3.19943e+05    1.02333e+06    3.10133e+02   -2.90942e+02
> 
>            Step           Time         Lambda
>          208500      417.00002        0.00000
> 
>    Energies (kJ/mol)
>            Bond          Angle    Proper Dih. Ryckaert-Bell.          LJ-14
>     2.30968e+04    4.27637e+02    1.17359e+01   -5.59373e+01    1.39477e+02
>      Coulomb-14        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential
>     1.80067e+03    6.20751e+04   -4.20609e+05   -4.67193e+04   -3.79833e+05
>     Kinetic En.   Total Energy  Conserved En.    Temperature Pressure (bar)
>     1.03323e+05   -2.76511e+05    1.07569e+06    4.22292e+02    2.94766e+02
> 
>    So my questions are:
> 
>     1) Am I missing some obvious step in setting up a stable simulation?

You'll have to describe your equilibration protocol in detail, including your 
.mdp file, because energy minimization alone is not enough.

>     2) Is it true that the high internal force present at the end of the 
> initial energy minimization could be the root of the problem?

Possibly.  However, I routinely minimize my protein structures to Fmax < 1000 
(usually they reach about 700-900 Fmax), and I've not had a problem.

>     3) If so, is there an obvious method for relaxing the system into a 
> more stable state prior to the main MD run?

See my response to #1.  Thorough equilibration is required, but we don't know 
all what you've done.

-Justin

> 
> Warren Gallin
> _______________________________________________
> 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/search before posting!
> Please don't post (un)subscribe requests to the list. Use thewww 
> interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================