[gmx-users] Re: MD Simulation Errors for Ethylene Glycol-Water System
Justin A. Lemkul
jalemkul at vt.edu
Wed Aug 5 22:03:29 CEST 2009
Nancy wrote:
> Hello,
>
> I have energy minimised my ethylene glycol system. I try to run
> position restrained equilibration with constant volume and temperature
> on it with the following .mdp file:
>
> =============================================
> title = Ethanediol equilibration 1
> define = -DPOSRES
> ; Run parameters
> integrator = md
> nsteps = 5000
> dt = 0.0002
5000 steps * 0.0002 ps/step = 1 ps. I doubt that length of time is sufficient
to equilibrate even the simplest system.
<snip>
> rlist = 1.0
> rcoulomb = 1.0
> rvdw = 1.3
If you're still using Gromos96, I'd suggest you read the original derivation of
the force field to get the correct values for these parameters, most notably rvdw.
<snip>
> I ran the equilibration with the command:
>
> $ grompp -f nvt.mdp -c em.gro -p ethanediol.top -o nvt.tpr
>
> where "em.gro" is the output coordinates of the minimisation, and
> "ethanediol.top" is the topology of the system. I then run the
> equilibration with the following command:
>
You haven't told us much about the energy minimization. What Fmax and Epot did
your system converge to? Do you believe those values to be acceptable?
> $ mdrun -v -deffnm nvt
>
> The system was "blowing up", and mdrun would not complete the
> equilbration. I tried decreasing dt, and although mdrun finishes
> without error, the output trajectory (viewed through ngmx) shows that
> the system has "blown up". Additionally, the structure of the
> ethanediol solute seems to be incorrect throughout the equilbration.
>
Have you corrected for periodicity effects? The trajectory will look scrambled
if any molecules cross periodic boundaries.
Beyond this, I doubt anyone can provide any magical insights into what's going
wrong without a better description of what "incorrect" means in the context of
the structure.
-Justin
> I am unsure what is going wrong, please help.
>
> Thank you.
>
> Nancy
>
>
> On Wed, Aug 5, 2009 at 9:13 AM, Bruce D. Ray <brucedray at yahoo.com
> <mailto:brucedray at yahoo.com>> wrote:
>
> On Tuesday, August 4, 2009, at 6:30:57 PM, Nancy
> <nancy5villa at gmail.com <mailto: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
>
>
>
> _______________________________________________
> gmx-users mailing list gmx-users at gromacs.org
> <mailto: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 the
> www interface or send it to gmx-users-request at gromacs.org
> <mailto:gmx-users-request at gromacs.org>.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> 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 the
> www 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
========================================
More information about the gromacs.org_gmx-users
mailing list