[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