[gmx-users] Re: MD Simulation Errors for Ethylene Glycol-Water System

[Bruce D. Ray]

On Tuesday, August 4, 2009, at 6:30:57 PM, Nancy <nancy5villa at gmail.com> wrote:
> I am trying to simulate a simple system of ethylene
glycol (ethane-1,2-diol) solvated in a water box.
> I converted the pdb
file to a mol2 file and used topolbuild 1.2.1 to generate topology
files:
> 
> $ topolbuild -n ethanediol -dir .../topolbuild1_2_1/dat/gromacs -ff gmx53a6

I prefer to use an absolute reference to the parameters directory, but
I also have my molecules files in directories outside of the topolbuild
directory.

> I the used editconf to enlarge the box.  I then solvated the molecule using the following command:

I presume that you checked the topology produced to be sure that parameters
were found for all atoms, bonds, angles, dihedrals, and impropers in the molecule
before proceeding.  This is always a necessary step with any automatic topology
generation program because the force field chosen might not have parameters
for everything in your molecule.  (In this case, I believe that ethylene glycol
should not present any problems, but always checking the topology generated
is a good habit to develop.)

> $ genbox -cp ethanediol_box.gro -cs spc216.gro -o ethanediol_solv.gro -box 3 3 3 -p ethanediol.top
> 
> where "ethanediol_box.gro" is the structure of the molecule.  I used grompp to generate the mdrun input file:
> 
> $ grompp -f grompp.mdp -c ethanediol_solv.gro -p ethanediol.top -o ethanediol.tpr
> 
> grompp.mdp is the following:
> 
> ========================================
> title                    = Ethanediol
> cpp                      = /lib/cpp
> include                  = -I../top
> define                   = 
> integrator               = md
> dt                       = 0.002
> nsteps                   = 50000
> nstxout                  = 1
> nstvout                  = 5
> 
nstlog                   = 5
> nstenergy                = 10
> nstxtcout                = 1
> xtc_grps                 = EDO  SOL
> energygrps               = EDO  SOL
> nstlist                  = 10
> ns_type                  = grid
> 
rlist                    = 0.8

rlist seems short to me.  Probably ought to use at least 1.0

> coulombtype              = cut-off

Its best to use PME with topolbuild generated topologies.

> rcoulomb                 = 1.4
> rvdw                     = 0.8

rvdw is too small.  Probably ought to be 1.2 to 1.4

> tcoupl                   = Berendsen
> tc-grps                  = EDO  SOL
> tau_t                    = 0.1  0.1
> 
ref_t                    = 300  300
> Pcoupl                   = Berendsen
> tau_p                    = 1.0
> compressibility          = 4.5e-5
> ref_p                    = 1.0
> gen_vel                  = yes
> gen_temp                 = 300
> 
gen_seed                 = 173529
> constraints              = all-bonds
> ========================================
> 
> where "EDO" refers to ethylene glycol.  grompp creates several notes:
> 
> NOTE 1 [file grompp.mdp, line unknown]:
> 
  The Berendsen thermostat does not generate the correct kinetic energy
>   distribution. You might want to consider using the V-rescale thermostat.

Something to think about.

> NOTE 2 [file grompp.mdp, line unknown]:
>   You are using a plain Coulomb cut-off, which might produce artifacts.
> 
  You might want to consider using PME electrostatics.
> 
> NOTE 3 [file grompp.mdp, line unknown]:
>   This run will generate roughly 2500 Mb of data

That is a large amount of data.  Perhaps make nstxout larger.
Also, this looks like a production simulation.  I presume you did
an energy minimization step and a position restrained equilibration
before this, but you did not show us anything about these.  Did they
have any errors?

> I then run the simulation:
> 
> $ mdrun -s ethanediol.tpr -o traj.trr -x traj.xtc -v
> 
> However mdrun produces several errors and warnings:
> 
> t = 0.036 ps: Water molecule starting at atom 2659 can not be settled.
> Check for bad contacts and/or reduce the timestep.
> Wrote pdb files with previous and current coordinates
> 
> t = 0.038 ps: Water molecule starting at atom 2659 can not be settled.
> Check for bad contacts and/or reduce the timestep.
> Wrote pdb files with previous and current coordinates
> 
> t = 0.040 ps: Water molecule starting at atom 2659 can not be settled.
> 
Check for bad contacts and/or reduce the timestep.
> Wrote pdb files with previous and current coordinates
> 
> Step 21  Warning: pressure scaling more than 1%, mu: 1.0372 1.0372 1.0372
> 
> Step 29, time 0.058 (ps)  LINCS WARNING
> 
relative constraint deviation after LINCS:
> rms 2.714709, max 3.541824 (between atoms 2 and 3)
> bonds that rotated more than 30 degrees:
>  atom 1 atom 2  angle  previous, current, constraint length
>       1      4   88.8    0.6606   0.6855      0.1520
>       2      3   90.6    0.4305   0.4542      0.1000
>       1      2   89.2    0.6123   0.4690      0.1435
>       5      6   93.4    0.4193   0.1741      0.1000
>       4      5   89.2    0.6156   0.5033      0.1435
> 
===================================
> 

It is blowing up.  Did you see anything in the energy minimization or the
position restrained equilibration steps that seemed unusual?

> 
> I
am not sure where the errors are occuring (I am also wondering if the
absence of
> explicit non-polar hydrogens in the .gro files are
relevant).  Please adivse.

 
Gromacs force fields are united atom force fields.  As such, explicit non-polar hydrogens
are supposed to be removed.


I hope these few remarks help.

Sincerely,

-- 
Bruce D. Ray, Ph.D.
Associate Scientist
IUPUI
Physics Dept.
402 N. Blackford St.
Indianapolis, IN  46202-3273


      
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090805/10ae28a3/attachment.html>