[gmx-users] rigid tetrahedral molecule

Justin A. Lemkul jalemkul at vt.edu
Thu Apr 28 12:40:35 CEST 2011 


Sanku M wrote:
> Thanks Mark for your comment. But, as far as degrees of freedom is 
> concerned, if there is a tetrahedral molecule as I shown below with atom 
> #1 ( B) being center of the tetrahedron , if we had defined the 
> molecules in terms of bonds and angles ( in stead of constraints ), we 
> would have 4 bonds ( 1-2,1-3,1-4,1-5 ) and 6 angles ( 
> 2-1-3,2-1-4,2-1-5,3-1-4,3-1-5 , 4-1-5 )at 109.45 degrees  , thus total ( 
> 4+6)=10 descriptors .
> In that case, the simulation goes fine . But, then wonder, why, for 
> constraints, we can put only 9 descriptors .  
> Also, I have tried all possible combination to get 9 constraints, but , 
> in each case, gromacs crashes. However, defining 4 bonds and 6 angles 
> using high force-constant lets the simulation go fine.  So, on a second 
> thought, I wonder whether defining 9 bond-constraints to get a stable 
> tetrahedral molecule is at all possible. Do I need to go for some sort 
> of compromise,  where 4 bonds ( 1-2,1-3,1-4,1-5 )  are defined as 
> constraints and 6 angles with high force-constant are defined to 
> maintain tetrahedral nature ? 
> I am not sure what should be the best possible 9 constraints .

For an example, see the CHCL3 definition in the Gromos96 53A6 force field.

-Justin


>                       1  opls_1014   1    BF4    B       1      0.8276
>>>      2  opls_1015   1    BF4    F1      1      -0.4569
>>>      3  opls_1015   1    BF4    F2      1      -0.4569
>>>      4  opls_1015   1    BF4    F3      1      -0.4569
>>>      5  opls_1015   1    BF4    F4      1      -0.4569
> 
> ------------------------------------------------------------------------
> *From:* Mark Abraham <Mark.Abraham at anu.edu.au>
> *To:* Discussion list for GROMACS users <gmx-users at gromacs.org>
> *Sent:* Thu, April 28, 2011 1:00:42 AM
> *Subject:* Re: [gmx-users] rigid tetrahedral molecule
> 
> On 4/28/2011 3:49 PM, Sanku M wrote:
>>  I went through the LINCS manual .  But, I am still struggling with 
>> coming up with the idea of putting correct constraint to maintain the 
>> rigidity of tetrahedral molecule .  I seem to understand  from your 
>> suggestion that the tetrahedral can be seen as a combination of  4 
>>  coupled triangles.( or am I still wrong about it ?)
> 
> I think it's seven coupled triangles, but that's not a relevant way to 
> think about it. You have 5 atoms, so 3N-6 means 9 degrees of freedom, so 
> you need 9 independent descriptors of relative atomic positions, so 9 
> constraints. You were trying 10 and 5.
> 
> Because those constraints form triangles, you may need to take care with 
> LINCS to get a stable simulation. Read up on the details here, I don't 
> remember them.
> 
>> In that case, am I supposed to use multiple settle to keep the 
>> molecule in a tetrahedral fashion ? I am sorry but if you can explain 
>> it in bit more details, I might get the point.
> 
> Forget about SETTLE - it was just an example to illustrate that this is 
> not easy to do right. There's a specialised algorithm for rigid water 
> (with three coupled bond constraints) because it is fiddly to get such 
> things right (and fast).
> 
> Mark
> 
>> ------------------------------------------------------------------------
>> *From:* Mark Abraham <Mark.Abraham at anu.edu.au>
>> *To:* Discussion list for GROMACS users <gmx-users at gromacs.org>
>> *Sent:* Wed, April 27, 2011 11:14:34 PM
>> *Subject:* Re: [gmx-users] rigid tetrahedral molecule
>>
>> On 4/28/2011 1:54 PM, Sanku M wrote:
>>> Hi,
>>>   I tried to keep the geometry of the BF4 fixed by using constraints 
>>> using lincs. But , unfortunately, my simulation is crashing 
>>> immediately and if I try minimization with only 2 molecules, it 
>>> provides a lot of LINCS warning and generate a lot of step*.pdb file 
>>> . If I try to visualize the minimized snapshot in VMD, it looks like 
>>> all the distances I tried to constrain decreased drastically. 
>>> Finally, trying MD run with this "minimized" configuration results in 
>>> crashing due to bad contacts.
>>>
>>>  I am sure I am doing something wrong and it might be that my itp 
>>> file is wrong . So any help will be highly appreciated.
>>> Here is the details of what I did.
>>>
>>>       The geometry of the molecule is tetrahedral with B at the 
>>> center and 4 F atoms is surrounding it in a tetrahedral manner.
>>> I first generated a itp file for BF4 which is shown below: I first 
>>> got the LJ parameters and charges for B and F atom and put them in 
>>> ffoplsnb.itp file as new atom types opls_1014 and opls_1015 
>>> . Initially I tried to put contsraint along all bonds ( i.e among F 
>>> atoms as well ). But, grompp provides warning that number of 
>>> constraint is more than number of degrees of freedom. So, I reduced 
>>> number of constraints by only putting constraint among B and F. But, 
>>> it did not work either.
>>
>> Sure, you need as many constraints as 
>> http://en.wikipedia.org/wiki/Degrees_of_freedom_%28physics_and_chemistry%29
>>
>> You should also do your homework about using LINCS and coupled 
>> triangles of constraints, as I suggested last time.
>>
>> Mark
>>
>>> Here is  the .itp file I wrote for rigid BF4 . It will be great if 
>>> someone can point me what I am doing wrong.
>>>
>>>   [ moleculetype ]
>>> ; molname       nrexcl
>>> BF4             3
>>>
>>> [ atoms ]
>>> #ifdef _FF_OPLS
>>>      1  opls_1014   1    BF4    B       1      0.8276
>>>      2  opls_1015   1    BF4    F1      1      -0.4569
>>>      3  opls_1015   1    BF4    F2      1      -0.4569
>>>      4  opls_1015   1    BF4    F3      1      -0.4569
>>>      5  opls_1015   1    BF4    F4      1      -0.4569
>>> #endif
>>> [ constraints ]
>>>   1  2   1   0.146
>>>   1  3   1   0.146
>>>   1  4   1   0.146
>>>   1  5   1   0.146
>>>   ;  2  3   1   0.238
>>> ;  2  4   1   0.238
>>> ;  2  5   1   0.238
>>> ;  3  4   1   0.238
>>> ;  3  5   1   0.238
>>> ;  4  5   1   0.238
>>>
>>> [ exclusions ]
>>> 1       2       3    4  5
>>> 2       1       3    4  5
>>> 3       1       2    4  5
>>> 4       1       2    3  5
>>> 5       1       2    3  4
>>>
>>>
>>> ------------------------------------------------------------------------
>>> *From:* Mark Abraham <Mark.Abraham at anu.edu.au>
>>> *To:* Discussion list for GROMACS users <gmx-users at gromacs.org>
>>> *Sent:* Wed, April 27, 2011 8:39:23 PM
>>> *Subject:* Re: [gmx-users] rigid tetrahedral molecule
>>>
>>> On 4/28/2011 11:25 AM, Justin A. Lemkul wrote:
>>> >
>>> >
>>> > Sanku M wrote:
>>> >> Hi,
>>> >>  I am interested in simulating a anionic molecule BF4(-)  ( Boron 
>>> tetrafluoride).  In the paper which developed the parameters for this 
>>> molecule, it is mentioned that it has been used as 'rigid' molecule 
>>> i.e the molecule only has non-bonding interaction but there was no 
>>> intramolecular motion as the geometry was fixed.
>>> >>  I am trying to simulate this molecule in gromacs treating it as 
>>> rigid.  But, I was looking for best way to 'rigidify' this molecule.
>>> >>
>>> >> I was wondering whether using LINCS to constrain all B-F and F-F 
>>> bonds will be good enough .  Or, Should I use virtual sites ? If I 
>>> really need to use virtual site, will it be something like TIP5P 
>>> water model ?
>>> >>  Can someone suggest the best wayout ?
>>> >>
>>> >
>>> > Constraints should do the trick, but probably the best approach is 
>>> to simply contact the authors who developed the model and ask how 
>>> they did it.  Then you know you're exactly reproducing what they did.
>>>
>>> Yep.
>>>
>>> Be aware that the coupled constraints make life tricky, and you 
>>> should read up in the manual and literature for how best to use 
>>> P-LINCS in such cases. Algorithms like SETTLE for rigid water exist 
>>> for a reason...
>>>
>>> Mark
>>>
>>>
>>>
>>> -- gmx-users mailing list    gmx-users at gromacs.org 
>>> <mailto:gmx-users at gromacs.org>
>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>> Please search the archive at 
>>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>> Please don't post (un)subscribe requests to the list. Use the www 
>>> interface or send it to gmx-users-request at gromacs.org 
>>> <mailto:gmx-users-request at gromacs.org>.
>>> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================

More information about the gromacs.org_gmx-users mailing list