Hello,<br><br>I am trying to equilibrate and run MD on a system of solvated ethylene glycol (ethanediol). However I am running into numerous problems.<br><br>First, I try to minimise the system. I start with an ethanediol.mol2 file which contains the structure of the molecule.<br>
<br>$ .../topolbuild1_2_1/src/topolbuild -n ethanediol -dir .../topolbuild1_2_1/dat/gromacs -ff gmx53a6<br><br>the above outputs the following files:<br><br>ethanediol.gro<br>ethanediol.log<br>ethanediolMOL.mol2<br>ethanediol.top<br>
ffethanediol.itp<br>posreethanediol.itp<br><br>Note: the ethanediol.log file contains a section that has several lines with asterisks:<br><br>============================<br> Angles Force Field Results<br>Angle Atoms force angle method measured<br>
1 H12- C1- H11 ****** ****** ****** 109.403<br> 2 C2- C1- H11 ****** ****** ****** 109.456<br> 3 O1- C1- H11 ****** ****** ****** 109.438<br> 4 C2- C1- H12 ****** ****** ****** 109.537<br>
5 O1- C1- H12 ****** ****** ****** 109.468<br> 6 O1- C1- C2 320.0 109.50 1 109.526<br> 7 C1- C2- O2 320.0 109.50 1 109.526<br> 8 C1- C2- H22 ****** ****** ****** 109.453<br>
9 C1- C2- H21 ****** ****** ****** 109.537<br> 10 H22- C2- O2 ****** ****** ****** 109.439<br> 11 H21- C2- O2 ****** ****** ****** 109.468<br> 12 C2- O2- HO2 450.0 109.50 1 106.864<br>
13 H21- C2- H22 ****** ****** ****** 109.404<br> 14 C1- O1- HO1 450.0 109.50 1 106.864<br>============================<br><br>$ editconf -f ethanediol.gro -o ethanediol_box.gro -c -d 1.0 -bt cubic<br>
<br>$ genbox -cp ethanediol_box.gro -cs spc216.gro -o ethanediol_solv.gro -p ethanediol.top<br><br>$ grompp -f minim.mdp -c ethanediol_solv.gro -p ethanediol.top -o em.tpr<br><br>the "minim.mdp" file is:<br><br>
============================<br>define = -DFLEXIBLE<br>integrator = steep<br>emtol = 10.0<br>emstep = 0.01<br>nsteps = 2000<br><br>nstlist = 1<br>ns_type = grid<br>rlist = 1.0<br>
coulombtype = PME<br>rcoulomb_switch = 1.0<br>rvdw_switch = 1.3<br>pbc = xyz<br><br>pme_order = 4<br>constraints = none<br><br>nstxout = 1<br>nstvout = 1<br>nstenergy = 1<br>nstlog = 1<br>
nstcomm = 1<br>Tcoupl = no<br>Pcoupl = no<br>gen_vel = no<br>============================<br><br>$ mdrun -v -deffnm em<br><br>the last few lines of the output of the minimisation are:<br>
<br>============================<br>Step= 1990, Dmax= 1.1e-05 nm, Epot= -3.38051e+04 Fmax= 4.80304e+04, atom= 559<br>Step= 1991, Dmax= 1.3e-05 nm, Epot= -3.38053e+04 Fmax= 7.73723e+04, atom= 795<br>Step= 1992, Dmax= 1.6e-05 nm, Epot= -3.38059e+04 Fmax= 7.31046e+04, atom= 559<br>
Step= 1994, Dmax= 9.5e-06 nm, Epot= -3.38069e+04 Fmax= 2.67253e+04, atom= 739<br>Step= 1995, Dmax= 1.1e-05 nm, Epot= -3.38076e+04 Fmax= 4.96799e+04, atom= 559<br>Step= 1996, Dmax= 1.4e-05 nm, Epot= -3.38078e+04 Fmax= 8.04458e+04, atom= 795<br>
Step= 1997, Dmax= 1.6e-05 nm, Epot= -3.38084e+04 Fmax= 7.56146e+04, atom= 559<br>Step= 1999, Dmax= 9.9e-06 nm, Epot= -3.38094e+04 Fmax= 2.68226e+04, atom= 739<br>Step= 2000, Dmax= 1.2e-05 nm, Epot= -3.38102e+04 Fmax= 5.49575e+04, atom= 559<br>
<br>writing lowest energy coordinates.<br><br>Steepest Descents did not converge to Fmax < 10 in 2001 steps.<br>Potential Energy = -3.3810168e+04<br>Maximum force = 5.4957539e+04 on atom 559<br>Norm of force = 2.7119446e+03<br>
============================<br><br>I have already tried increasing the number of steps, but I found no way of lowering the forces. I have found that by limiting the number of solvent molecules to ~5, I can manage to acheive lower energies, however, the amount of time it takes for a larger number of water molecules increases exponentially.<br>
<br>So, I am trying to proceed to equilibration:<br><br>$ grompp -f nvt.mdp -c em.gro -p ethanediol.top -o nvt.tpr<br><br>"nvt.mdp" is:<br><br>============================<br>title = Ethanediol NVT equilibration <br>
define = -DFLEXIBLE<br><br>integrator = md<br>nsteps = 2000<br>dt = 0.002<br><br>nstxout = 10<br>nstvout = 10<br>nstenergy = 10<br>nstlog = 10<br><br>continuation = no<br>
constraint_algorithm = lincs<br>constraints = all-angles<br>lincs_iter = 1<br>lincs_order = 4<br><br>ns_type = grid<br>nstlist = 5<br>rlist = 1.0<br>rcoulomb = 1.0<br>rvdw = 1.2<br>
<br>coulombtype = PME<br>pme_order = 4<br>fourierspacing = 0.16<br><br>tcoupl = V-rescale<br>tc-grps = EDO SOL<br>tau_t = 0.1 0.1<br>ref_t = 300 300<br><br>pcoupl = no<br><br>pbc = xyz<br>
<br>DispCorr = EnerPres<br><br>gen_vel = yes<br>gen_temp = 300<br>gen_seed = -1<br>============================<br><br>Unfortunately, I receive numerous errors such as this one:<br><br>============================<br>
Step 3, time 0.006 (ps) LINCS WARNING<br>relative constraint deviation after LINCS:<br>rms 0.067249, max 2.300835 (between atoms 691 and 692)<br>bonds that rotated more than 30 degrees:<br> atom 1 atom 2 angle previous, current, constraint length<br>
406 407 33.5 0.1005 0.0956 0.1000<br> 692 693 91.2 0.1631 0.3470 0.1633<br> 691 692 88.9 0.0999 0.3301 0.1000<br> 691 693 86.0 0.0999 0.0881 0.1000<br>
794 795 35.4 0.1707 0.1587 0.1633<br> 793 795 59.5 0.1063 0.0995 0.1000<br>Wrote pdb files with previous and current coordinates<br>Segmentation fault<br>============================<br>
<br>I believe that equilibration is failing because of the large forces on the molecules during minimisation (although I am not sure, as I was able to run a simulation without proper minimisation before).<br><br>If at all possible, please advise on how to reduce the energies and forces during minimisation.<br>
<br><br>Thank you.<br><br>Nancy<br><br><br><br><br><br>