[gmx-users] Re: MD Simulation Errors for Ethylene Glycol-Water System
Nancy
nancy5villa at gmail.com
Wed Aug 5 21:06:37 CEST 2009
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
; Output control
nstxout = 10
nstvout = 10
nstenergy = 10
nstlog = 10
; Bond parameters
continuation = no
constraint_algorithm = lincs
constraints = all-angles
lincs_iter = 1
lincs_order = 4
; Neighborsearching
ns_type = grid
nstlist = 5
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.3
; Electrostatics
coulombtype = PME
pme_order = 4
fourierspacing = 0.16
; Temperature coupling is on
tcoupl = V-rescale
tc-grps = EDO SOL
tau_t = 0.1 0.1
ref_t = 300 300
; Pressure coupling is off
pcoupl = no
; Periodic boundary conditions
pbc = xyz
; Dispersion correction
DispCorr = EnerPres
; Velocity generation
gen_vel = yes
gen_temp = 300
gen_seed = -1
=============================================
and the following position restraint file:
=============================================
[ position_restraints ]
; atom type fx fy fz
1 1 1000 1000 1000
2 1 1000 1000 1000
4 1 1000 1000 1000
5 1 1000 1000 1000
=============================================
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:
$ 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.
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> wrote:
> 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
>
>
>
> _______________________________________________
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090805/4523fe41/attachment.html>
More information about the gromacs.org_gmx-users
mailing list