Hello,<br><br>I obtained the ethylene glycol (1,2-ethanediol) structure from the URL: <a href="http://www.rcsb.org/pdb/files/ligand/EDO_ideal.pdb">http://www.rcsb.org/pdb/files/ligand/EDO_ideal.pdb</a><br><br>I converted the EDO_ideal.pdb file into ethanediol.mol2 using UCSF Chimera (<a href="http://www.cgl.ucsf.edu/chimera/">http://www.cgl.ucsf.edu/chimera/</a>).<br>
<br>I then use topolbuild 1.2 (written by Dr. Bruce D. Ray) to generate the topologies:<br><br>$ .../topolbuild1_2_2/src/topolbuild -n ethanediol -dir .../topolbuild1_2_2/dat/gromacs -ff gmx53a6<br><br>which outputs the files:<br>
<br>ethanediol.gro<br>ethanediol.log<br>ethanediolMOL.mol2<br>ethanediol.top<br>ffethanediol.itp<br>posreethanediol.itp<br><br>I modified the &quot;ffethanediol.itp&quot; file to read:<br><br>=======================================<br>
#define _FF_GROMOS96<br>#define _FF_GROMOS53A6<br>#define _FF_USER<br><br>[ defaults ]<br>;nbfunc     comb-rule      gen-pairs     fudgeLJ      fudgeQQ<br>      1             1             no         1.0          1.0<br><br>
#include &quot;ffG53a6nb.itp&quot;<br>=======================================<br><br>The lack of hydrogen bonds between the solute and solvent molecules was due to the lack of charges in the generated topology file &quot;ethanediol.top&quot;.  So I changed the &quot;atoms&quot; section of the topology file (the original topology file is at the end of this message), and this causes hydrogen bonds between the solute and solvent in numbers comparable to that between the solvent molecules.<br>
<br>=======================================<br> [ atoms ]<br>;  nr    type   resnr   residu   atom   cgnr    charge      mass<br>    1       CH2     1       EDO      C1    1   0.17600  14.02700       ;   0.0000000<br>    2        OA     1       EDO      O1    1   -0.5740  15.99940       ;   0.0000000<br>
    3        H      1       EDO     HO1    1   0.39800   1.00800       ;   0.0000000<br>    4       CH2     1       EDO      C2    2   0.17600  14.02700       ;   0.0000000<br>    5        OA     1       EDO      O2    2   -0.5740  15.99940       ;   0.0000000<br>
    6        H      1       EDO     HO2    2   0.39800   1.00800       ;   0.0000000<br>; total molecule charge =   0.0000000<br>=======================================<br><br>I obtained the charge values from the methanol tutoral included with Gromacs (.../share/gromacs/tutor/methanol/methanol.itp).  I then enlarge the box and solvate the molecule:<br>
<br>$ editconf -f ethanediol.gro -o ethanediol_box.gro -box 2.05 2.05 2.05<br><br>$ genbox -cp ethanediol_box.gro -cs spc216.gro -o ethanediol_solv.gro -p ethanediol.top -box 2.05 2.05 2.05<br><br>The solvated system consists of one ethylene glycol and 232 water molecules.  I then configure the minimisation:<br>
<br>$ grompp -f em.mdp -c ethanediol_solv.gro -p ethanediol.top -o em.tpr<br><br>with the &quot;em.mdp&quot; file:<br><br>=======================================<br>constraints         =  none<br>integrator          =  steep<br>
nsteps              =  10000<br>;<br>;    Energy minimizing stuff<br>;<br><br>emtol               =  2.0<br>emstep              =  0.01<br><br>nstlist             =  2<br>coulombtype         =  PME<br>nstcomm             =  2<br>
<br>ns_type             =  grid<br>;rlist               =  1.0<br>;rcoulomb            =  1.0<br>;rvdw                =  1.0<br>nstxout             =  2<br><br>pbc                 =  xyz<br>pme_order           =  4<br>=======================================<br>
<br>I then run the minimisation:<br><br>$ mdrun -v -deffnm em<br><br>The last few lines of output are:<br><br>=======================================<br>Step= 8069, Dmax= 1.5e-04 nm, Epot= -1.22107e+04 Fmax= 8.11110e+01, atom= 5<br>
Step= 8070, Dmax= 1.8e-04 nm, Epot= -1.22107e+04 Fmax= 5.63915e+01, atom= 5<br>Step= 8071, Dmax= 2.2e-04 nm, Epot= -1.22108e+04 Fmax= 1.20155e+02, atom= 5<br>Step= 8072, Dmax= 2.6e-04 nm, Epot= -1.22108e+04 Fmax= 7.97799e+01, atom= 5<br>
Step= 8074, Dmax= 1.6e-04 nm, Epot= -1.22108e+04 Fmax= 4.58228e+01, atom= 5<br>Step= 8078, Dmax= 2.3e-05 nm, Epot= -1.22108e+04 Fmax= 2.84317e+01, atom= 571<br>Step= 8082, Dmax= 3.5e-06 nm, Epot= -1.22108e+04 Fmax= 2.59027e+01, atom= 571<br>
Step= 8085, Dmax= 1.1e-06 nm, Epot= -1.22108e+04 Fmax= 2.51778e+01, atom= 5<br>Stepsize too small, or no change in energy.<br>Converged to machine precision,<br>but not to the requested precision Fmax &lt; 2<br><br>Double precision normally gives you higher accuracy.<br>
You might need to increase your constraint accuracy, or turn<br>off constraints alltogether (set constraints = none in mdp file)<br><br>writing lowest energy coordinates.<br><br>Steepest Descents converged to machine precision in 8086 steps,<br>
but did not reach the requested Fmax &lt; 2.<br>Potential Energy  = -1.2210801e+04<br>Maximum force     =  2.5902737e+01 on atom 5<br>Norm of force     =  2.9901226e+00<br>=======================================<br><br>I viewed the minimisation trajectory using VMD, and there are hydrogen bonds between the solute and solvent.<br>
<br>Then I proceed to equilibration:<br><br>$ grompp -f nvt.mdp -c em.gro -p ethanediol.top -o nvt.tpr<br><br>where the &quot;nvt.mdp&quot; file is:<br><br>=======================================<br>title        = Ethanediol NVT equilibration <br>
<br>integrator    = md<br>nsteps        = 50000<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><br>coulombtype    = PME<br>pme_order    = 4<br>fourierspacing    = 0.16<br><br>tcoupl        = V-rescale<br>
tc-grps        = System<br>tau_t        = 0.1<br>ref_t        = 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>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?<br>
<br>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?<br><br>Thank you.<br><br>Nancy<br><br><br><br><br>===========ethanediol.top===========<br>
;<br>; Topology from .mol2 file<br>; topolbuild<br>;<br>; The force field files to be included<br>#include &quot;ffethanediol.itp&quot;<br><br> [ moleculetype ]<br>; name  nrexcl<br>EDO_ideal.pdb   3<br><br> [ atoms ]<br>
;  nr    type   resnr   residu   atom   cgnr    charge      mass<br>    1       CH2     1       EDO      C1    1   0.00000  14.02700       ;   0.0000000<br>    2        OA     1       EDO      O1    1   0.00000  15.99940       ;   0.0000000<br>
    3        H      1       EDO     HO1    1   0.00000   1.00800       ;   0.0000000<br>    4       CH2     1       EDO      C2    2   0.00000  14.02700       ;   0.0000000<br>    5        OA     1       EDO      O2    2   0.00000  15.99940       ;   0.0000000<br>
    6        H      1       EDO     HO2    2   0.00000   1.00800       ;   0.0000000<br>; total molecule charge =   0.0000000<br><br> [ bonds ]<br>;   ai  aj   funct      b0          kb<br>       1     4   2     0.15200     5430000.       ;    C1-    C2<br>
       4     5   2     0.14350     6100000.       ;    C2-    O2<br>       5     6   2     0.10000    15700000.       ;    O2-   HO2<br>       1     2   2     0.14350     6100000.       ;    C1-    O1<br>       2     3   2     0.10000    15700000.       ;    O1-   HO1<br>
<br> [ pairs ]<br>       2       5  1       ;    O1-    O2<br>       1       6  1       ;    C1-   HO2<br>       4       3  1       ;    C2-   HO1<br><br>[ angles ]<br>; ai  aj  ak  funct      th0         cth<br>     2     1     4   2     109.500    320.0000     ;    O1-    C1-    C2<br>
     1     4     5   2     109.500    320.0000     ;    C1-    C2-    O2<br>     4     5     6   2     109.500    450.0000     ;    C2-    O2-   HO2<br>     1     2     3   2     109.500    450.0000     ;    C1-    O1-   HO1<br>
<br>[ dihedrals ]<br>; ai  aj   ak  al  funct    phi0          cp      mult<br>     2     1     4     5   1       0.000       2.530         3  ; dih     O1-    C1-    C2-    O2<br>     2     1     4     5   1       0.000       5.350         1  ; dih     O1-    C1-    C2-    O2<br>
     1     4     5     6   1       0.000       1.260         3  ; dih     C1-    C2-    O2-   HO2<br>     4     1     2     3   1       0.000       1.260         3  ; dih     C2-    C1-    O1-   HO1<br><br>; Include Position restraint file<br>
; WARNING: Position restraints and distance restraints ought not be done together<br>#ifdef POSRES<br>#include &quot;posreethanediol.itp&quot;<br>#endif<br><br>; Include water topology<br>#include &quot;spce.itp&quot;<br>
<br>#ifdef POSRES_WATER<br>; Position restraint for each water oxygen<br>[ position_restraints ]<br>;  i funct       fcx        fcy        fcz<br>   1    1       1000       1000       1000<br>#endif<br><br>; Include generic topology for ions<br>
#include &quot;ions.itp&quot;<br><br> [ system ]<br>; title from mol2 input<br>EDO_ideal.pdb<br><br> [ molecules ]<br>; molecule name    nr.<br>EDO_ideal.pdb           1<br><br>===========ethanediol.top===========<br><br>
<br><br><br><br><br><div class="gmail_quote">On Sat, Aug 8, 2009 at 7:46 PM, Nagy, Peter I. <span dir="ltr">&lt;<a href="mailto:PNAGY@utnet.utoledo.edu">PNAGY@utnet.utoledo.edu</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">



<p><font size="2">   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?<br>
If O-C-C-O is gauche, then existence of two internal hydrogen bonds is possible, but it is<br>
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<br>

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.<br>

   If you really do not find solute-solvent hydrogen bonds, your ethylene-glycol charges may be<br>
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.<br>
Peter Nagy<br>
The University of Toledo,<br>
Toledo, OH 43606<div><div></div><div class="h5"><br>
-----Original Message-----<br>
From: <a href="mailto:gmx-users-bounces@gromacs.org" target="_blank">gmx-users-bounces@gromacs.org</a> on behalf of Nancy<br>
Sent: Sat 8/8/2009 3:33 PM<br>
To: Discussion list for GROMACS users<br>
Subject: [gmx-users] Fatal Equilibration Errors<br>
I have successfully minimised and equilibrated ethylene glycol in a water<br>
box.  I have noticed that there seem to be no hydrogen bonds between the<br>
solute and solvent, but there are hydrogen bonds forming and breaking<br>
between solvent molecules.  Is this a normal behavior during minimsation and<br>
I am now trying to run MD on glycerol.  I start with a &quot;glycerol.mol2&quot; file<br>
which contains the structure of glycerol.  I am using the following commands<br>
to setup and run minimisation and equilibration on glycerol:<br>
$ .../topolbuild1_2_1/src/topolbuild -n glycerol -dir<br>
.../topolbuild1_2_1/dat/gromacs -ff gmx53a6<br>
I modified the &quot;defaults&quot; part of the &quot;ffglycerol.itp&quot; (generated by<br>
topolbuild) file to read:<br>
[ defaults ]<br>
;nbfunc     comb-rule      gen-pairs     fudgeLJ      fudgeQQ<br>
      1             1             no         1.0          1.0<br>
$ editconf -f glycerol.gro -o glycerol_box.gro -box 2.65 2.65 2.65<br>
$ genbox -cp glycerol_box.gro -cs spc216.gro -o glycerol_solv.gro -p<br>
glycerol.top -box 2.65 2.65 2.65<br>
$ grompp -f em.mdp -c glycerol_solv.gro -p glycerol.top -o em.tpr<br>
my &quot;em.mdp&quot; file is as follows:<br>
constraints         =  none<br>
integrator          =  steep<br>
nsteps              =  10000<br>
emtol               =  10.0<br>
emstep              =  0.01<br>
nstlist             =  2<br>
coulombtype         =  PME<br>
nstcomm             =  2<br>
ns_type             =  grid<br>
rlist               =  1.0<br>
rcoulomb            =  1.0<br>
rvdw                =  1.3<br>
nstxout             =  2<br>
pbc                 =  xyz<br>
pme_order           =  4<br>
$ mdrun -v -deffnm em<br>
the last few lines of output are:<br>
Step= 6990, Dmax= 4.7e-04 nm, Epot= -3.34590e+04 Fmax= 9.13695e+01, atom= 4<br>
Step= 6992, Dmax= 2.8e-04 nm, Epot= -3.34590e+04 Fmax= 1.40130e+02, atom= 5<br>
Step= 6993, Dmax= 3.4e-04 nm, Epot= -3.34591e+04 Fmax= 8.57813e+01, atom= 4<br>
Step= 6994, Dmax= 4.0e-04 nm, Epot= -3.34591e+04 Fmax= 2.41887e+02, atom= 5<br>
Step= 6995, Dmax= 4.8e-04 nm, Epot= -3.34592e+04 Fmax= 9.36159e+01, atom= 4<br>
Step= 6999, Dmax= 7.3e-05 nm, Epot= -3.34592e+04 Fmax= 4.95606e+01, atom= 4<br>
Step= 7006, Dmax= 1.4e-06 nm, Epot= -3.34592e+04 Fmax= 4.89221e+01, atom= 4<br>
Stepsize too small, or no change in energy.<br>
Converged to machine precision,<br>
but not to the requested precision Fmax &lt; 10<br>
Double precision normally gives you higher accuracy.<br>
You might need to increase your constraint accuracy, or turn<br>
off constraints alltogether (set constraints = none in mdp file)<br>
writing lowest energy coordinates.<br>
Steepest Descents converged to machine precision in 7007 steps,<br>
but did not reach the requested Fmax &lt; 10.<br>
Potential Energy  = -3.3459230e+04<br>
Maximum force     =  4.9560604e+01 on atom 4<br>
Norm of force     =  3.4399371e+00<br>
As I believe these forces to be acceptable, I proceed to equilibration:<br>
$ grompp -f nvt.mdp -c em.gro -p glycerol.top -o nvt.tpr<br>
where my &quot;nvt.mdp&quot; file is:<br>
title        = Glycerol NVT equilibration<br>
define        = -DPOSRES<br>
integrator    = md<br>
nsteps        = 5000<br>
dt        = 0.002<br>
nstxout        = 10<br>
nstvout        = 10<br>
nstenergy    = 10<br>
nstlog        = 10<br>
continuation    = no<br>
constraint_algorithm = lincs<br>
constraints    = all-angles<br>
lincs_iter    = 1<br>
lincs_order    = 4<br>
ns_type        = grid<br>
nstlist        = 5<br>
rlist        = 1.0<br>
rcoulomb    = 1.0<br>
rvdw        = 1.3<br>
coulombtype    = PME<br>
pme_order    = 4<br>
fourierspacing    = 0.16<br>
tcoupl        = V-rescale<br>
tc-grps        = GOL SOL<br>
tau_t        = 0.1 0.1<br>
ref_t        = 300 300<br>
pcoupl        = no<br>
pbc        = xyz<br>
DispCorr    = EnerPres<br>
gen_vel        = yes<br>
gen_temp    = 300<br>
gen_seed    = -1<br>
where &quot;GOL&quot; refers to GlycerOL.<br>
$ mdrun -v -deffnm nvt<br>
When I execute the above command, I get numerous errors of the following<br>
Step 0, time 0 (ps)  LINCS WARNING<br>
relative constraint deviation after LINCS:<br>
rms 0.326337, max 0.500490 (between atoms 1 and 3)<br>
bonds that rotated more than 30 degrees:<br>
 atom 1 atom 2  angle  previous, current, constraint length<br>
starting mdrun &#39;glycerol.pdb in water&#39;<br>
5000 steps,     10.0 ps.<br>
Step 0, time 0 (ps)  LINCS WARNING<br>
relative constraint deviation after LINCS:<br>
rms 0.663850, max 1.227203 (between atoms 1 and 7)<br>
bonds that rotated more than 30 degrees:<br>
 atom 1 atom 2  angle  previous, current, constraint length<br>
      1      3  172.3    0.2091   0.1644      0.2004<br>
      7      9  164.1    0.2082   0.1625      0.2004<br>
      2      3  124.5    0.1022   0.0787      0.1000<br>
      1      2  148.3    0.1477   0.0882      0.1435<br>
      5      6   30.4    0.0992   0.1265      0.1000<br>
      8      9  134.1    0.1014   0.0660      0.1000<br>
      7      8  138.2    0.1476   0.0991      0.1435<br>
Warning: 1-4 interaction between 2 and 7 at distance 6.002 which is larger<br>
than the 1-4 table size 2.300 nm<br>
These are ignored for the rest of the simulation<br>
This usually means your system is exploding,<br>
if not, you should increase table-extension in your mdp file<br>
or with user tables increase the table size<br>
t = 0.004 ps: Water molecule starting at atom 1234 can not be settled.<br>
Check for bad contacts and/or reduce the timestep.<br>
Wrote pdb files with previous and current coordinates<br>
There were 12 inconsistent shifts. Check your topology<br>
Segmentation fault<br>
I don&#39;t know what is happening wrong.  Please advise.<br>
Thank you.<br>

gmx-users mailing list    <a href="mailto:gmx-users@gromacs.org">gmx-users@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-users" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-users</a><br>
Please search the archive at <a href="http://www.gromacs.org/search" target="_blank">http://www.gromacs.org/search</a> before posting!<br>
Please don&#39;t post (un)subscribe requests to the list. Use the<br>
www interface or send it to <a href="mailto:gmx-users-request@gromacs.org">gmx-users-request@gromacs.org</a>.<br>
Can&#39;t post? Read <a href="http://www.gromacs.org/mailing_lists/users.php" target="_blank">http://www.gromacs.org/mailing_lists/users.php</a><br></blockquote></div><br>