<html><head><style type="text/css"><!-- DIV {margin:0px;} --></style></head><body><div style="font-family:times new roman,new york,times,serif;font-size:10pt"><span style="font-family: times new roman,new york,times,serif;">On </span><font style="font-family: times new roman,new york,times,serif;" face="Tahoma" size="2">Tuesday, August 4, 2009, at 6:30:57 PM</font><span style="font-family: times new roman,new york,times,serif;">, </span><font style="font-family: times new roman,new york,times,serif;" face="Tahoma" size="2">Nancy <nancy5villa@gmail.com> wrote:</font><br>> I am trying to simulate a simple system of ethylene
glycol (ethane-1,2-diol) solvated in a water box.<br>> I converted the pdb
file to a mol2 file and used topolbuild 1.2.1 to generate topology
files:<br>> <br>> $ topolbuild -n ethanediol -dir .../topolbuild1_2_1/dat/gromacs -ff gmx53a6<br>
<br>I prefer to use an absolute reference to the parameters directory, but<br>I also have my molecules files in directories outside of the topolbuild<br>directory.<br><br>> I the used editconf to enlarge the box. I then solvated the molecule using the following command:<br><br>I presume that you checked the topology produced to be sure that parameters<br>were found for all atoms, bonds, angles, dihedrals, and impropers in the molecule<br>before proceeding. This is always a necessary step with any automatic topology<br>generation program because the force field chosen might not have parameters<br>for everything in your molecule. (In this case, I believe that ethylene glycol<br>should not present any problems, but always checking the topology generated<br>is a good habit to develop.)<br><br>> $ genbox -cp ethanediol_box.gro -cs spc216.gro -o ethanediol_solv.gro -box 3 3 3 -p ethanediol.top<br>> <br>> where "ethanediol_box.gro"
is the structure of the molecule. I used grompp to generate the mdrun input file:<br>
> <br>> $ grompp -f grompp.mdp -c ethanediol_solv.gro -p ethanediol.top -o ethanediol.tpr<br>> <br>> grompp.mdp is the following:<br>> <br>> ========================================<br>> title = Ethanediol<br>> cpp = /lib/cpp<br>> include = -I../top<br>> define = <br>> integrator = md<br>>
dt = 0.002<br>> nsteps = 50000<br>> nstxout = 1<br>> nstvout = 5<br>>
nstlog = 5<br>> nstenergy = 10<br>> nstxtcout = 1<br>> xtc_grps = EDO SOL<br>> energygrps = EDO SOL<br>> nstlist = 10<br>> ns_type = grid<br>>
rlist = 0.8<br><br>rlist seems short to me. Probably ought to use at least 1.0<br><br>> coulombtype = cut-off<br><br>Its best to use PME with topolbuild generated topologies.<br><br>> rcoulomb = 1.4<br>> rvdw = 0.8<br><br>rvdw is too small. Probably ought to be 1.2 to 1.4<br><br>> tcoupl = Berendsen<br>> tc-grps =
EDO SOL<br>> tau_t = 0.1 0.1<br>>
ref_t = 300 300<br>> Pcoupl = Berendsen<br>> tau_p = 1.0<br>> compressibility = 4.5e-5<br>> ref_p = 1.0<br>> gen_vel = yes<br>> gen_temp = 300<br>>
gen_seed = 173529<br>> constraints = all-bonds<br>> ========================================<br>> <br>> where "EDO" refers to ethylene glycol. grompp creates several notes:<br>> <br>> NOTE 1 [file grompp.mdp, line unknown]:<br>>
The Berendsen thermostat does not generate the correct kinetic energy<br>> distribution. You might want to consider using the V-rescale thermostat.<br><br>Something to think about.<br><br>> NOTE 2 [file grompp.mdp, line unknown]:<br>> You are using a plain Coulomb cut-off, which might produce artifacts.<br>>
You might want to consider using PME electrostatics.<br>> <br>> NOTE 3 [file grompp.mdp, line unknown]:<br>> This run will generate roughly 2500 Mb of data<br><br>That is a large amount of data. Perhaps make nstxout larger.<br>Also, this looks like a production simulation. I presume you did<br>an energy minimization step and a position restrained equilibration<br>before this, but you did not show us anything about these. Did they<br>have any errors?<br><br>> I then run the simulation:<br>> <br>> $ mdrun -s ethanediol.tpr -o traj.trr -x traj.xtc -v<br>
> <br>> However mdrun produces several errors and warnings:<br>> <br>> t = 0.036 ps: Water molecule starting at atom 2659 can not be settled.<br>> Check for bad contacts and/or reduce the timestep.<br>> Wrote pdb files with previous and current coordinates<br>
> <br>> t = 0.038 ps: Water molecule starting at atom 2659 can not be settled.<br>> Check for bad contacts and/or reduce the timestep.<br>> Wrote pdb files with previous and current coordinates<br>> <br>> t = 0.040 ps: Water molecule starting at atom 2659 can not be settled.<br>>
Check for bad contacts and/or reduce the timestep.<br>> Wrote pdb files with previous and current coordinates<br>> <br>> Step 21 Warning: pressure scaling more than 1%, mu: 1.0372 1.0372 1.0372<br>> <br>> Step 29, time 0.058 (ps) LINCS WARNING<br>>
relative constraint deviation after LINCS:<br>> rms 2.714709, max 3.541824 (between atoms 2 and 3)<br>> bonds that rotated more than 30 degrees:<br>> atom 1 atom 2 angle previous, current, constraint length<br>> 1 4 88.8 0.6606 0.6855 0.1520<br>
> 2 3 90.6 0.4305 0.4542 0.1000<br>> 1 2 89.2 0.6123 0.4690 0.1435<br>> 5 6 93.4 0.4193 0.1741 0.1000<br>> 4 5 89.2 0.6156 0.5033 0.1435<br>>
===================================<br>> <br><br>It is blowing up. Did you see anything in the energy minimization or the<br>position restrained equilibration steps that seemed unusual?<br><br>> <br>> I
am not sure where the errors are occuring (I am also wondering if the
absence of<br>> explicit non-polar hydrogens in the .gro files are
relevant). Please adivse.<br><div> <br>Gromacs force fields are united atom force fields. As such, explicit non-polar hydrogens<br>are supposed to be removed.<br><br><br>I hope these few remarks help.<br><br>Sincerely,<br><br></div>-- <br>Bruce D. Ray, Ph.D.<br>Associate Scientist<br>IUPUI<br>Physics Dept.<br>402 N. Blackford St.<br>Indianapolis, IN 46202-3273<br><br></div><br>
</body></html>