[gmx-users] Fatal Equilibration Errors
Justin A. Lemkul
jalemkul at vt.edu
Sun Aug 9 23:36:17 CEST 2009
Nancy wrote:
<snip>
Sounds like things are going fine.
> I viewed the trajectory of the equilibration, and it looks reasonable.
> The O-C-C-O conformation is mostly trans (anti) during the entire NVT
> equilibration, and the OH hydrogens form numerous hydrogen bonds with
> the solvent (it seems there are no intramolecular hydrogen bonds).
> Therefore, the hydrogen bonds do form, given the charges adopted from
> the methanol topology (O: -0.574, H: +0.398, CH2: +0.176). Do these
> charge values seem reasonable?
>
"Reasonable" and "valid" are two separate ideas that may or may not be true
simultaneously. If you look, for example, at the charges on a serine side chain
in 53a6, all of your charges are very different. Since methanol was used to
parameterize serine, I don't know where the charges in the tutorial file came
from. The tutorial files claim to use 43a1, in any case, so be careful of
blindly copying parameters from some other location.
To convince a skeptical audience (i.e., reviewers) that your parameters are
valid, at least use the parameters from the functional groups present in the
force field parameter set you claim to be using. Beyond that, you should really
be validating the parameters based on the original force field methodology,
i.e., calculations of thermodynamic observables, described in full in the 2004
JCC paper dedicated to the 53a5 and 53a6 parameter sets.
> The ethylene glycol lacks explicit non-polar hydrogens; is that normal
> for this type of system, force field (based on GROMOS96 53a6) and SPC/E
> water model?
>
Since G53a6 is a UA force field, by definition, yes. If this comes as a
surprise, you should certainly do a little background reading about this type of
force field and why you want to use it. Choosing a force field for simulation
is perhaps the most important part of running your experiments. If you don't
know the ins and outs (benefits, limitations, accuracy) of a particular
parameter set, you're asking for trouble if a reviewer asks you why you felt it
was best :)
-Justin
> Thank you.
>
> Nancy
>
>
>
>
> ===========ethanediol.top===========
> ;
> ; Topology from .mol2 file
> ; topolbuild
> ;
> ; The force field files to be included
> #include "ffethanediol.itp"
>
> [ moleculetype ]
> ; name nrexcl
> EDO_ideal.pdb 3
>
> [ atoms ]
> ; nr type resnr residu atom cgnr charge mass
> 1 CH2 1 EDO C1 1 0.00000 14.02700
> ; 0.0000000
> 2 OA 1 EDO O1 1 0.00000 15.99940
> ; 0.0000000
> 3 H 1 EDO HO1 1 0.00000 1.00800
> ; 0.0000000
> 4 CH2 1 EDO C2 2 0.00000 14.02700
> ; 0.0000000
> 5 OA 1 EDO O2 2 0.00000 15.99940
> ; 0.0000000
> 6 H 1 EDO HO2 2 0.00000 1.00800
> ; 0.0000000
> ; total molecule charge = 0.0000000
>
> [ bonds ]
> ; ai aj funct b0 kb
> 1 4 2 0.15200 5430000. ; C1- C2
> 4 5 2 0.14350 6100000. ; C2- O2
> 5 6 2 0.10000 15700000. ; O2- HO2
> 1 2 2 0.14350 6100000. ; C1- O1
> 2 3 2 0.10000 15700000. ; O1- HO1
>
> [ pairs ]
> 2 5 1 ; O1- O2
> 1 6 1 ; C1- HO2
> 4 3 1 ; C2- HO1
>
> [ angles ]
> ; ai aj ak funct th0 cth
> 2 1 4 2 109.500 320.0000 ; O1- C1- C2
> 1 4 5 2 109.500 320.0000 ; C1- C2- O2
> 4 5 6 2 109.500 450.0000 ; C2- O2- HO2
> 1 2 3 2 109.500 450.0000 ; C1- O1- HO1
>
> [ dihedrals ]
> ; ai aj ak al funct phi0 cp mult
> 2 1 4 5 1 0.000 2.530 3 ;
> dih O1- C1- C2- O2
> 2 1 4 5 1 0.000 5.350 1 ;
> dih O1- C1- C2- O2
> 1 4 5 6 1 0.000 1.260 3 ;
> dih C1- C2- O2- HO2
> 4 1 2 3 1 0.000 1.260 3 ;
> dih C2- C1- O1- HO1
>
> ; Include Position restraint file
> ; WARNING: Position restraints and distance restraints ought not be done
> together
> #ifdef POSRES
> #include "posreethanediol.itp"
> #endif
>
> ; Include water topology
> #include "spce.itp"
>
> #ifdef POSRES_WATER
> ; Position restraint for each water oxygen
> [ position_restraints ]
> ; i funct fcx fcy fcz
> 1 1 1000 1000 1000
> #endif
>
> ; Include generic topology for ions
> #include "ions.itp"
>
> [ system ]
> ; title from mol2 input
> EDO_ideal.pdb
>
> [ molecules ]
> ; molecule name nr.
> EDO_ideal.pdb 1
>
> ===========ethanediol.top===========
>
>
>
>
>
>
> On Sat, Aug 8, 2009 at 7:46 PM, Nagy, Peter I. <PNAGY at utnet.utoledo.edu
> <mailto:PNAGY at utnet.utoledo.edu>> wrote:
>
>
> It is extremely surprising that you have not found solute-solvent
> hydrogen bonds in aqueous solution. What is the O-C-C-O
> conformation? If it is trans, then what do the OH hydrogens do?
> If O-C-C-O is gauche, then existence of two internal hydrogen bonds
> is possible, but it is
> far from the quantum mechanically most stable form in the gas phase.
> I think that at least one solute-solvent hydrogen bond should have
> been seen. There are papers in the literature which do show
> solute-solvent hydrogen bonds in water. The exciting question is:
> there is one
> H-O-C-C trans moiety in the gas phase, and the other H is H-bonded
> to this oxygen, when the O-C-C-O skeleton is gauche. Then the
> question of solvation is, whether this intramolecular hydrogen bond
> is maintained or disrupts. In aqueous solution it may, whereas in
> chloroform or CCl4 both OH may be internally H-bonded. Experiments
> found 12% O-C-C-O trans and 88% O-C-C-O gauche conformers in aqueous
> solution. Modeling should find what are the positions of the hydroxy
> hydrogens.
> If you really do not find solute-solvent hydrogen bonds, your
> ethylene-glycol charges may be
> questionable. What charge values did you use? A possibility is that
> the glycol OH charges cannot compete with those of the water model
> (SPC, TIPxP or what) in forming intermolecular hydrogen bonds.
>
> Peter Nagy
> The University of Toledo,
> Toledo, OH 43606
>
>
> -----Original Message-----
> From: gmx-users-bounces at gromacs.org
> <mailto:gmx-users-bounces at gromacs.org> on behalf of Nancy
> Sent: Sat 8/8/2009 3:33 PM
> To: Discussion list for GROMACS users
> Subject: [gmx-users] Fatal Equilibration Errors
>
> Hello,
>
> I have successfully minimised and equilibrated ethylene glycol in a
> water
> box. I have noticed that there seem to be no hydrogen bonds between the
> solute and solvent, but there are hydrogen bonds forming and breaking
> between solvent molecules. Is this a normal behavior during
> minimsation and
> equilibration?
>
> I am now trying to run MD on glycerol. I start with a
> "glycerol.mol2" file
> which contains the structure of glycerol. I am using the following
> commands
> to setup and run minimisation and equilibration on glycerol:
>
> $ .../topolbuild1_2_1/src/topolbuild -n glycerol -dir
> .../topolbuild1_2_1/dat/gromacs -ff gmx53a6
>
> I modified the "defaults" part of the "ffglycerol.itp" (generated by
> topolbuild) file to read:
>
> ===========================
> [ defaults ]
> ;nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> 1 1 no 1.0 1.0
> ===========================
>
> $ editconf -f glycerol.gro -o glycerol_box.gro -box 2.65 2.65 2.65
> $ genbox -cp glycerol_box.gro -cs spc216.gro -o glycerol_solv.gro -p
> glycerol.top -box 2.65 2.65 2.65
> $ grompp -f em.mdp -c glycerol_solv.gro -p glycerol.top -o em.tpr
>
> my "em.mdp" file is as follows:
>
> ===========================
> constraints = none
> integrator = steep
> nsteps = 10000
>
> emtol = 10.0
> emstep = 0.01
>
> nstlist = 2
> coulombtype = PME
> nstcomm = 2
>
> ns_type = grid
> rlist = 1.0
> rcoulomb = 1.0
> rvdw = 1.3
> nstxout = 2
>
> pbc = xyz
> pme_order = 4
> ===========================
>
> $ mdrun -v -deffnm em
>
> the last few lines of output are:
>
> ===============================================================
> Step= 6990, Dmax= 4.7e-04 nm, Epot= -3.34590e+04 Fmax= 9.13695e+01,
> atom= 4
> Step= 6992, Dmax= 2.8e-04 nm, Epot= -3.34590e+04 Fmax= 1.40130e+02,
> atom= 5
> Step= 6993, Dmax= 3.4e-04 nm, Epot= -3.34591e+04 Fmax= 8.57813e+01,
> atom= 4
> Step= 6994, Dmax= 4.0e-04 nm, Epot= -3.34591e+04 Fmax= 2.41887e+02,
> atom= 5
> Step= 6995, Dmax= 4.8e-04 nm, Epot= -3.34592e+04 Fmax= 9.36159e+01,
> atom= 4
> Step= 6999, Dmax= 7.3e-05 nm, Epot= -3.34592e+04 Fmax= 4.95606e+01,
> atom= 4
> Step= 7006, Dmax= 1.4e-06 nm, Epot= -3.34592e+04 Fmax= 4.89221e+01,
> atom= 4
> Stepsize too small, or no change in energy.
> Converged to machine precision,
> but not to the requested precision Fmax < 10
>
> Double precision normally gives you higher accuracy.
> You might need to increase your constraint accuracy, or turn
> off constraints alltogether (set constraints = none in mdp file)
>
> writing lowest energy coordinates.
>
> Steepest Descents converged to machine precision in 7007 steps,
> but did not reach the requested Fmax < 10.
> Potential Energy = -3.3459230e+04
> Maximum force = 4.9560604e+01 on atom 4
> Norm of force = 3.4399371e+00
> ===============================================================
>
> As I believe these forces to be acceptable, I proceed to equilibration:
>
> $ grompp -f nvt.mdp -c em.gro -p glycerol.top -o nvt.tpr
>
> where my "nvt.mdp" file is:
>
> ===========================
> title = Glycerol NVT equilibration
> define = -DPOSRES
>
> integrator = md
> nsteps = 5000
> dt = 0.002
>
> nstxout = 10
> nstvout = 10
> nstenergy = 10
> nstlog = 10
>
> continuation = no
> constraint_algorithm = lincs
> constraints = all-angles
> lincs_iter = 1
> lincs_order = 4
>
> ns_type = grid
> nstlist = 5
> rlist = 1.0
> rcoulomb = 1.0
> rvdw = 1.3
>
> coulombtype = PME
> pme_order = 4
> fourierspacing = 0.16
>
> tcoupl = V-rescale
> tc-grps = GOL SOL
> tau_t = 0.1 0.1
> ref_t = 300 300
>
> pcoupl = no
>
> pbc = xyz
>
> DispCorr = EnerPres
>
> gen_vel = yes
> gen_temp = 300
> gen_seed = -1
> ===========================
>
> where "GOL" refers to GlycerOL.
>
> $ mdrun -v -deffnm nvt
>
> When I execute the above command, I get numerous errors of the following
> type:
>
> ===============================================================
> Step 0, time 0 (ps) LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.326337, max 0.500490 (between atoms 1 and 3)
> bonds that rotated more than 30 degrees:
> atom 1 atom 2 angle previous, current, constraint length
> starting mdrun 'glycerol.pdb in water'
> 5000 steps, 10.0 ps.
>
> Step 0, time 0 (ps) LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.663850, max 1.227203 (between atoms 1 and 7)
> bonds that rotated more than 30 degrees:
> atom 1 atom 2 angle previous, current, constraint length
> 1 3 172.3 0.2091 0.1644 0.2004
> 7 9 164.1 0.2082 0.1625 0.2004
> 2 3 124.5 0.1022 0.0787 0.1000
> 1 2 148.3 0.1477 0.0882 0.1435
> 5 6 30.4 0.0992 0.1265 0.1000
> 8 9 134.1 0.1014 0.0660 0.1000
> 7 8 138.2 0.1476 0.0991 0.1435
>
> Warning: 1-4 interaction between 2 and 7 at distance 6.002 which is
> larger
> than the 1-4 table size 2.300 nm
> These are ignored for the rest of the simulation
> This usually means your system is exploding,
> if not, you should increase table-extension in your mdp file
> or with user tables increase the table size
>
> t = 0.004 ps: Water molecule starting at atom 1234 can not be settled.
> Check for bad contacts and/or reduce the timestep.
> Wrote pdb files with previous and current coordinates
>
> There were 12 inconsistent shifts. Check your topology
> Segmentation fault
> ===============================================================
>
> I don't know what is happening wrong. Please advise.
>
> Thank you.
>
> Nancy
>
>
> _______________________________________________
> 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