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

[Jennifer Williams]

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.