[gmx-users] Exploding poly-serine

[Warren Gallin]

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?
	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?
	3) If so, is there an obvious method for relaxing the system into a  
more stable state prior to the main MD run?

Warren Gallin