[gmx-users] Simulating Multiple Solute Particles

Nathan K Houtz nhoutz at purdue.edu
Wed Nov 12 08:58:56 CET 2014

Dr. Lemkul,

Thank you again - your tutorials are very helpful for me. I want to get through your Virtual Sites tutorial but I got stuck near the beginning. Gromacs thinks the minim.mdp file is for an earlier version but the introduction page said it's for 5.0 and later. I'm using gromacs 5.0.2, so I thought I should be fine. I copied this command:

> gmx grompp -f minim.mdp -c box.pdb -p topol.top -o min.tpr

after having made my box, of course, and got the following error:

>Command line:
>  gmx grompp -f minim.mdp -c box.pdb -p topol.top -o min.tpr
>NOTE 1 [file minim.mdp, line 19]:
>  minim.mdp did not specify a value for the .mdp option "cutoff-scheme".
>  Probably it was first intended for use with GROMACS before 4.6. In 4.6,
>  the Verlet scheme was introduced, but the group scheme was still the
>  default. The default is now the Verlet scheme, so you will observe
>  different behaviour.
>Ignoring obsolete mdp entry 'title'
>ERROR 1 [file minim.mdp]:
>  With Verlet lists only full pbc or pbc=xy with walls is supported
>ERROR 2 [file minim.mdp]:
>  With Verlet lists nstlist should be larger than 0
>NOTE 2 [file minim.mdp]:
>  With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note
>  that with the Verlet scheme, nstlist has no effect on the accuracy of
>  your simulation.

I tried to figure out how I could modify it myself but I don't really know what I'm doing. The sample mdp file (http://manual.gromacs.org/online/mdp.html) is much longer than the one you linked (http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/vsites/Files/minim.mdp) and I'm not sure which fields are necessary. Is the file outdated, or have I made a mistake?


----- Original Message -----
From: "Justin Lemkul" <jalemkul at vt.edu>
To: gmx-users at gromacs.org
Sent: Monday, November 10, 2014 7:40:58 AM
Subject: Re: [gmx-users] Simulating Multiple Solute Particles

On 11/10/14 3:07 AM, Nathan K Houtz wrote:
> Thank you again, Dr. Lemkul! I went with the first option because my advisor wants me to use tip3p water. Unfortunately, I immediately got myself stuck on the next step, but this time I think the solution may not be as trivial. The command to output em.tpr seems to work without any errors. The command,

Well, strictly speaking, either option works because the parameters of TIP3P are 
not force field-dependent, just translated between different conventions for 
nonbonded parameters, but that's just a pedantic issue.

>> gmx mdrun -deffnm em
> gives me a lot of errors. (It stops after 1000). Here's a small sample:
>> Steepest Descents:
>>    Tolerance (Fmax)   =  1.00000e+03
>>    Number of steps    =        50000
>> Step 1, time 0.001 (ps)  LINCS WARNING
>> relative constraint deviation after LINCS:
>> rms 0.373893, max 0.619462 (between atoms 1747 and 1748)
>> bonds that rotated more than 30 degrees:
>> atom 1 atom 2  angle  previous, current, constraint length
>>       4      5   74.8    0.1130   0.0501      0.1204
>>      11     12   75.0    0.1132   0.0510      0.1204
>>      18     19   74.9    0.1129   0.0497      0.1204
>>      25     26   75.6    0.1126   0.0488      0.1204
>>      32     33   74.9    0.1132   0.0513      0.1204
>>      39     40   75.7    0.1129   0.0503      0.1204
> etc. There are 7 atoms (the methyl group is a united atom, otherwise there would be 10) in each TTA molecule, so you'll notice that it's complaining about the bond between the 4th and 5th atoms on every one of them. My best guess is that my topology file has insufficient or incorrect constraints, allowing this particular bond to rotate. The entire molecule is supposed to be rigid. I did constrain all the distances, but not the angles or dihedrals. One of the example topology files I used did the same: constrained distances but not angles or dihedrals. I was skeptical because it seems to allow some degrees of freedom but I assumed gromacs somehow took care of it. Even in this instance, it doesn't make sense to me that only one of the bonds is problematic if none of them are constrained. So, I'm unsure of what exactly is being left out. The bond in question is between the double bonded C=O. My molecule looks like this:
> H              O
>    \          //
> H - C ≡ C - C
>    /          \
> H             O-H
> My topology file is available from a previous message in this thread (minus the last force field correction), linked here from the archive: https://www.mail-archive.com/gromacs.org_gmx-users@maillist.sys.kth.se/msg07776.html
> If anyone knows directly how to fix the error, I would greatly appreciate it, but helping me understand the error would also be very useful: Does "relative constraint deviation" mean that a constraint has been broken, or that there is no constraint where there should be one? Also, when it's saying the bonds are rotating more than 30 degrees, is that the bond angle changing or does it mean that a dihedral is changing by that much?

Your system is http://www.gromacs.org/Documentation/Terminology/Blowing_Up 
because the topology is probably unstable.  Triple bonds need to be handled 
using virtual sites.  See previous discussions about linear construction 
(examples are acetonitrile and others) and my tutorial on doing just this.


> Thanks again in advance - I really appreciate the help!
> N.H.
> ----- Original Message -----
> From: "Justin Lemkul" <jalemkul at vt.edu>
> To: gmx-users at gromacs.org
> Sent: Tuesday, November 4, 2014 7:45:39 AM
> Subject: Re: [gmx-users] Simulating Multiple Solute Particles
> On 11/4/14 12:45 AM, Nathan K Houtz wrote:
>> Thanks for your reply. However, I'm still confused. I thought that the command:
>> #include "oplsaa.ff/tip3p.itp"
>> is a call to the opls-aa force field. If this is not the correct way to include the force field parameters, how should I do that?
> That line does not call a force field, it calls a topology for TIP3P, which
> makes use of OPLS-AA parameters (well, more correctly TIP3P parameters as
> translated into OPLS-AA atom types and functional form).  Your options are
> either (1) #include "oplsaa.ff/forcefield.itp" at the start of the topology and
> remove the [defaults] directive in your TETR topology or (2) translate TIP3P
> into atom types that are self-contained in your topology.
> -Justin


Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441

Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.

More information about the gromacs.org_gmx-users mailing list