[gmx-users] Re: Best way to handle linear rigid molecules.
Jennifer Williams
Jennifer.Williams at ed.ac.uk
Wed Jul 15 17:44:40 CEST 2009
Hello,
Thanks for the help so far.
I?ve understood the concept of introducing the dummy atoms for linear
molecules with advice and by following an example for acetonitrile
which I found in the user forum (Feb 2003).
I now (hope) that I have an itp file of a have a system of two masses
(D1 and D2) that accurately reproduce the mass and moments of inertia
of CO2 and three 'dummy'atoms that accurately reproduce the normal
atomic interactions of CO2.
I?ve pasted the file below.
[ atomtypes ]
; type mass charge ptype c6 c12
D1 22.004975 0.0000 A 0.0000 0.00000000
D2 22.004975 0.0000 A 0.0000 0.00000000
C_CO 0.0000000 0.5406 V 0.0000 0.00000000
O_CO 0.0000000 -0.2703 V 0.29847 1.10765301
[ moleculetype ]
; name nrexcl
CO2 2
[ atoms ]
; nr type resnr residu atom cgnr charge mass
1 D1 1 CO2 DUM 1 0.0000 22.004975
2 D2 1 CO2 DUM 1 0.0000 22.004975
3 C_CO 1 CO2 C_CO 1 0.5406 0.0000000
4 O_CO 1 CO2 O_CO 1 -0.2703 0.0000000
5 O_CO 1 CO2 O_CO 1 -0.2703 0.0000000
[ constraints ]
; ai aj funct c0 c1
1 2 1 0.2006146447
[ virtualsites2]
; ai aj ak funct c0 c1
3 1 2 1 0.4263443
4 1 2 1 -0.07365569
5 1 2 1 0.926344308
[ exclusions ]
4 3 5
5 3 4
3 5 4
However, this now means that I need coordinates of the atoms D1 and D2
in my pdb file!
I currently have a pdb file (the output of another simulation program)
which contains the coordinates of the C,O and O atoms. i.e
ATOM 1077 C_CO CO2 1 -21.646 5.405 -4.644 1.00 0.00
ATOM 1078 O_CO CO2 1 -22.709 5.009 -4.226 1.00 0.00
ATOM 1079 O_CO CO2 1 -20.583 5.801 -5.062 1.00 0.00
I have constructed my .itp and .top files entirely by hand as the
residues are non standard so pdb2gmx doesn?t work for me.
How do I go about getting a pdb file containing the dummy positions?
My systems will have around 40 CO2 molecules so whatever I use, it has
to be automated. Is there some way of modelling CO2 as a linear rigid
molecule other than using the dummy mass atoms?
Thanks in advance,
Quoting Berk Hess <gmx3 at hotmail.com>:
>
> 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
Dr. Jennifer Williams
Institute for Materials and Processes
School of Engineering
University of Edinburgh
Sanderson Building
The King's Buildings
Mayfield Road
Edinburgh, EH9 3JL, United Kingdom
Phone: ++44 (0)131 650 4 861
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the gromacs.org_gmx-users
mailing list