[gmx-users] Simulating Multiple Solute Particles

[Magnus Lundborg]

On 11/10/2014 01:40 PM, Justin Lemkul wrote:
>
>
> 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.
I actually think Gromacs >= 4.6.4 (or something) can handle triple bonds 
without virtual sites, but I haven't tested it that much.

Cheers,

Magnus
>
> -Justin
>
>> 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
>>
>