[gmx-users] Re: Best way to handle linear rigid molecules.

Berk Hess gmx3 at hotmail.com
Tue Jul 14 15:53:52 CEST 2009


Hi,

This is quite wrong.

Your normal atoms are virtual sites defined from the two masses
(you did it the other way around).
If you put the masses at the oxygen postitions (which gives an incorrect
moment of intertia, so incorrect dynamics, but correct thermodynamics),
then your coefficients are 0.5 for C and 0 and 1 for the O's.

Berk

> Date: Tue, 14 Jul 2009 14:39:10 +0100
> From: Jennifer.Williams at ed.ac.uk
> To: gmx-users at gromacs.org
> Subject: [gmx-users] Re: Best way to handle linear rigid molecules.
> 
> Hi gmx users
> 
> I am still having problems modelling small linear molecules. Following  
> advice from the forum (see below), I tried to define CO2 using 2 dummy  
> atoms (with averaged masses) which I placed in the same positions as  
> my oxygen atoms. (This involves having to replicate coordinates for O  
> in my pdb file which will be fiddly when I have hundreds of CO2  
> molecules).
> 
> Here is my .itp file for CO2;
> 
> 
> [ atomtypes ]
> ;   type      mass    charge    ptype       c6            c12
>     C_CO    0.0000000     0.5406       A	 0.0000      0.00000000
>     O_CO    0.0000000    -0.2703       A	 0.29847     1.10765301
>     DUM    22.004975      0.0000       V	 0.0000      0.00000000
> 
> 
> [ moleculetype ]
> ; name  nrexcl
> CO2         2
> 
> [ atoms ]
> ;   nr  type    resnr   residu  atom    cgnr    charge	mass
> 1       C_CO     1       CO2    C_CO      1        0.5406  0.0000000
> 2       O_CO     1       CO2    O_CO      1       -0.2703  0.0000000
> 3       O_CO     1       CO2    O_CO      1       -0.2703  0.0000000
> 4       DUM      1       CO2    DUM       1        0.0000  22.004975
> 5       DUM      1       CO2    DUM       1        0.0000  22.004975
> 
> 
> [ constraints ]
> ;  ai  aj funct           c0           c1
> 2       3	1	   0.24176
> 
> [ virtualsites2]
> ;  ai    aj    ak       funct   c0      c1
>      4     2     3       1       0.0000
>      5     2     3       1       0.24176
> 
> I am not sure what the final terms in the virtual sites section stand  
> for and I have not found a lot of detail in the manual (section 4.7  
> and 5.5.2). I assumed that the first term defines the virtual atom  
> site, the second and third are the positions of the real atoms which  
> the virtual atom is placed in relation to. What does the function ?1?  
> stand for?  Here I tried to place the virtual atoms in the same  
> positions as my ?real? oxygens- I am not sure if this will be allowed?
> 
> However, when I run this I get:
> 
> 
> ERROR 1 [file MCM_CO2.top, line 8230]:
>    atom C_CO (Res CO2-1) has mass 0
> 
> 
> ERROR 2 [file MCM_CO2.top, line 8230]:
>    atom O_CO (Res CO2-1) has mass 0
> 
> 
> 
> ERROR 3 [file MCM_CO2.top, line 8230]:
>    atom O_CO (Res CO2-1) has mass 0
> 
> 
> 
> ERROR 4 [file MCM_CO2.top, line 8230]:
>    virtual site DUM (Res CO2-1) has non-zero mass 22.005
>         Check your topology.
> 
> 
> 
> ERROR 5 [file MCM_CO2.top, line 8230]:
>    virtual site DUM (Res CO2-1) has non-zero mass 22.005
>         Check your topology.
> 
> So the approach of modelling dummy atoms with a mass and the real  
> atoms as mass-less generates errors. How do I solve this? Importantly  
> I need to get the MSD of the CO2 molecules so I am a bit confused with  
> the approach of using virtual atom sites to model CO2. Which groups  
> would I then select for the MSD? The virtual site with the mass or the  
> real atoms site which is mass-less?  Is this approach realistic?
> 
> I have also tried to model CO2 2 different approaches
> 
> 1.	Using bond constraints and setting the bond angle to 180 and using  
> a large bending constant to keep the molecule effectively rigid.
> 
> Here I got the following error message and the run stopped after 2 timesteps;
> 
> Step 1, time 0.001 (ps)  LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.969889, max 10.563440 (between atoms 1072 and 1074)
> bonds that rotated more than 30 degrees:
>   atom 1 atom 2  angle  previous, current, constraint length
>     1072   1073   89.9    0.1210   0.6876      0.1209
>     1072   1074   92.4    0.1210   1.3978      0.1209
>     1075   1076   86.8    0.1209   0.1303      0.1209
>     1075   1077   86.8    0.1209   0.1307      0.1209
> Wrote pdb files with previous and current coordinates
> 
> 
> 2.	Using constraints to keep the C-O bond fixed to 0.12 and another  
> constraint between O and O of 0.24. I hoped this would keep the O-C-O  
> linear.
> 
> Here, the md.log file contained the following;
> 
> 
> 6 constraints are involved in constraint triangles,
> will apply an additional matrix expansion of order 4 for couplings
> between constraints inside triangles
> 
> The job ran to completion but on viewing the trajectory, the CO2  
> molecules are only approximately linear.
> 
> If anyone has an example file for a triatomic linear molecule like  
> CO2, I?d be very grateful if you could post it so I could use it as a  
> template,
> 
> I am using gromacs version 4.0.5
> 
> Thanks
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> Hi,
> 
> This question has been asked several times before on this list.
> The proper way to handle this in Gromacs is to introduce two dummy  
> mass particles
> in such a way the total mass and the moment of inertia is identical to CO2.
> Then you can construct the positions of the three mass-less C and O atoms from
> the positions of the masses using virtual sites (see the manual).
> 
> Berk
> 
> 
> 
> Quoting Jennifer Williams <Jennifer.Williams at ed.ac.uk>:
> 
> > Hello users,
> >
> > I am using gromacs v 4.0.5
> >
> > What is the best way to treat linear triatomic molecules such as CO2.
> > Is there some kind of rigid algorthym as in DL_POLY?
> >
> > In the asbsence of a ?rigid? algorthym. I initially intended to ?fix?
> > the C=O bond length using constraints (as below) and somehow ?fix? the
> > O=C=O bond angle to 180.
> >
> >
> > [ atomtypes ]
> > ;   type      mass    charge    ptype       c6            c12
> >    C_CO2    12.01115    0.5406      A	0.0000      0.00000000
> >    O_CO2    15.9994    -0.2703      A	0.29847     1.10765301
> >
> > [ moleculetype ]
> > ; name  nrexcl
> > CO2         2
> >
> > [ atoms ]
> > ;   nr  type    resnr   residu  atom    cgnr    charge	mass
> > 1       C_CO2     1       CO2    C_CO2      1        0.5406  12.01115
> > 2       O_CO2     1       CO2    O_CO2      1       -0.2703  15.9994
> > 3       O_CO2     1       CO2    O_CO2      1       -0.2703  15.9994
> >
> > [ constraints ]
> > ;  ai  aj funct           c0           c1
> > 1       2	1	   0.12088
> > 1       3	1	   0.12088
> >
> > [ angles ]
> > ;  ai    aj    ak       funct   c0      c1
> >     2     1     3       1       180.00
> >
> > First question... How do I go about fixing this angle to 180degrees?
> > Should I just make c1 extremely large? Is there a more elegant way of
> > going about this? Someone also mentioned using a virtual site
> > representation for the C. Any comments on this welcome.
> >
> > I read some posts which imply that the angle bending code apparently can't
> > handle angles of exactly 180 degrees.
> >
> > http://www.mail-archive.com/gmx-users@gromacs.org/msg09979.html
> >
> > This was a while ago so perhaps this has been dealt with in the current
> > version of gromacs?
> >
> >
> > Thanks
> 
> 
> 
> -- 
> The University of Edinburgh is a charitable body, registered in
> Scotland, with registration number SC005336.
> 
> 
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/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.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php

_________________________________________________________________
See all the ways you can stay connected to friends and family
http://www.microsoft.com/windows/windowslive/default.aspx
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090714/90d6e68c/attachment.html>


More information about the gromacs.org_gmx-users mailing list