<html>
<head>
<style>
.hmmessage P
{
margin:0px;
padding:0px
}
body.hmmessage
{
font-size: 10pt;
font-family:Verdana
}
</style>
</head>
<body class='hmmessage'>
Hi,<br><br>Write a (relatively) simple script that reads a pdb file<br>and prints for each set of three C O O lines two additional O lines: O O C O O.<br>Gromacs does not care about atom names, and the O positions will be<br>good enough that the initial constraining in mdrun will correct the mass positions<br>to exactly the right ones.<br>Then you do not need to run pdb2gmx again.<br><br>Berk<br><br>> Date: Wed, 15 Jul 2009 16:44:40 +0100<br>> From: Jennifer.Williams@ed.ac.uk<br>> To: gmx-users@gromacs.org; gmx3@hotmail.com<br>> Subject: RE: [gmx-users] Re: Best way to handle linear rigid molecules.<br>> CC: <br>> <br>> <br>> Hello,<br>> <br>> Thanks for the help so far.<br>> <br>> I?ve understood the concept of introducing the dummy atoms for linear <br>> molecules with advice and by following an example for acetonitrile <br>> which I found in the user forum (Feb 2003).<br>> <br>> I now (hope) that I have an itp file of a have a system of two masses <br>> (D1 and D2) that accurately reproduce the mass and moments of inertia <br>> of CO2 and three 'dummy'atoms that accurately reproduce the normal <br>> atomic interactions of CO2.<br>> <br>> I?ve pasted the file below.<br>> <br>> [ atomtypes ]<br>> ; type mass charge ptype c6 c12<br>> D1 22.004975 0.0000 A 0.0000 0.00000000<br>> D2 22.004975 0.0000 A 0.0000 0.00000000<br>> C_CO 0.0000000 0.5406 V 0.0000 0.00000000<br>> O_CO 0.0000000 -0.2703 V 0.29847 1.10765301<br>> <br>> [ moleculetype ]<br>> ; name nrexcl<br>> CO2 2<br>> <br>> [ atoms ]<br>> ; nr type resnr residu atom cgnr charge mass<br>> 1 D1 1 CO2 DUM 1 0.0000 22.004975<br>> 2 D2 1 CO2 DUM 1 0.0000 22.004975<br>> 3 C_CO 1 CO2 C_CO 1 0.5406 0.0000000<br>> 4 O_CO 1 CO2 O_CO 1 -0.2703 0.0000000<br>> 5 O_CO 1 CO2 O_CO 1 -0.2703 0.0000000<br>> <br>> [ constraints ]<br>> ; ai aj funct c0 c1<br>> 1 2 1 0.2006146447<br>> <br>> [ virtualsites2]<br>> ; ai aj ak funct c0 c1<br>> 3 1 2 1 0.4263443<br>> 4 1 2 1 -0.07365569<br>> 5 1 2 1 0.926344308<br>> <br>> [ exclusions ]<br>> 4 3 5<br>> 5 3 4<br>> 3 5 4<br>> <br>> However, this now means that I need coordinates of the atoms D1 and D2 <br>> in my pdb file!<br>> <br>> I currently have a pdb file (the output of another simulation program) <br>> which contains the coordinates of the C,O and O atoms. i.e<br>> <br>> ATOM 1077 C_CO CO2 1 -21.646 5.405 -4.644 1.00 0.00<br>> ATOM 1078 O_CO CO2 1 -22.709 5.009 -4.226 1.00 0.00<br>> ATOM 1079 O_CO CO2 1 -20.583 5.801 -5.062 1.00 0.00<br>> <br>> I have constructed my .itp and .top files entirely by hand as the <br>> residues are non standard so pdb2gmx doesn?t work for me.<br>> <br>> How do I go about getting a pdb file containing the dummy positions? <br>> My systems will have around 40 CO2 molecules so whatever I use, it has <br>> to be automated. Is there some way of modelling CO2 as a linear rigid <br>> molecule other than using the dummy mass atoms?<br>> <br>> Thanks in advance,<br>> <br>> <br>> Quoting Berk Hess <gmx3@hotmail.com>:<br>> <br>> ><br>> > Hi,<br>> ><br>> > This is quite wrong.<br>> ><br>> > Your normal atoms are virtual sites defined from the two masses<br>> > (you did it the other way around).<br>> > If you put the masses at the oxygen postitions (which gives an incorrect<br>> > moment of intertia, so incorrect dynamics, but correct thermodynamics),<br>> > then your coefficients are 0.5 for C and 0 and 1 for the O's.<br>> ><br>> > Berk<br>> ><br>> >> Date: Tue, 14 Jul 2009 14:39:10 +0100<br>> >> From: Jennifer.Williams@ed.ac.uk<br>> >> To: gmx-users@gromacs.org<br>> >> Subject: [gmx-users] Re: Best way to handle linear rigid molecules.<br>> >><br>> >> Hi gmx users<br>> >><br>> >> I am still having problems modelling small linear molecules. Following<br>> >> advice from the forum (see below), I tried to define CO2 using 2 dummy<br>> >> atoms (with averaged masses) which I placed in the same positions as<br>> >> my oxygen atoms. (This involves having to replicate coordinates for O<br>> >> in my pdb file which will be fiddly when I have hundreds of CO2<br>> >> molecules).<br>> >><br>> >> Here is my .itp file for CO2;<br>> >><br>> >><br>> >> [ atomtypes ]<br>> >> ; type mass charge ptype c6 c12<br>> >> C_CO 0.0000000 0.5406 A 0.0000 0.00000000<br>> >> O_CO 0.0000000 -0.2703 A 0.29847 1.10765301<br>> >> DUM 22.004975 0.0000 V 0.0000 0.00000000<br>> >><br>> >><br>> >> [ moleculetype ]<br>> >> ; name nrexcl<br>> >> CO2 2<br>> >><br>> >> [ atoms ]<br>> >> ; nr type resnr residu atom cgnr charge mass<br>> >> 1 C_CO 1 CO2 C_CO 1 0.5406 0.0000000<br>> >> 2 O_CO 1 CO2 O_CO 1 -0.2703 0.0000000<br>> >> 3 O_CO 1 CO2 O_CO 1 -0.2703 0.0000000<br>> >> 4 DUM 1 CO2 DUM 1 0.0000 22.004975<br>> >> 5 DUM 1 CO2 DUM 1 0.0000 22.004975<br>> >><br>> >><br>> >> [ constraints ]<br>> >> ; ai aj funct c0 c1<br>> >> 2 3 1 0.24176<br>> >><br>> >> [ virtualsites2]<br>> >> ; ai aj ak funct c0 c1<br>> >> 4 2 3 1 0.0000<br>> >> 5 2 3 1 0.24176<br>> >><br>> >> I am not sure what the final terms in the virtual sites section stand<br>> >> for and I have not found a lot of detail in the manual (section 4.7<br>> >> and 5.5.2). I assumed that the first term defines the virtual atom<br>> >> site, the second and third are the positions of the real atoms which<br>> >> the virtual atom is placed in relation to. What does the function ?1?<br>> >> stand for? Here I tried to place the virtual atoms in the same<br>> >> positions as my ?real? oxygens- I am not sure if this will be allowed?<br>> >><br>> >> However, when I run this I get:<br>> >><br>> >><br>> >> ERROR 1 [file MCM_CO2.top, line 8230]:<br>> >> atom C_CO (Res CO2-1) has mass 0<br>> >><br>> >><br>> >> ERROR 2 [file MCM_CO2.top, line 8230]:<br>> >> atom O_CO (Res CO2-1) has mass 0<br>> >><br>> >><br>> >><br>> >> ERROR 3 [file MCM_CO2.top, line 8230]:<br>> >> atom O_CO (Res CO2-1) has mass 0<br>> >><br>> >><br>> >><br>> >> ERROR 4 [file MCM_CO2.top, line 8230]:<br>> >> virtual site DUM (Res CO2-1) has non-zero mass 22.005<br>> >> Check your topology.<br>> >><br>> >><br>> >><br>> >> ERROR 5 [file MCM_CO2.top, line 8230]:<br>> >> virtual site DUM (Res CO2-1) has non-zero mass 22.005<br>> >> Check your topology.<br>> >><br>> >> So the approach of modelling dummy atoms with a mass and the real<br>> >> atoms as mass-less generates errors. How do I solve this? Importantly<br>> >> I need to get the MSD of the CO2 molecules so I am a bit confused with<br>> >> the approach of using virtual atom sites to model CO2. Which groups<br>> >> would I then select for the MSD? The virtual site with the mass or the<br>> >> real atoms site which is mass-less? Is this approach realistic?<br>> >><br>> >> I have also tried to model CO2 2 different approaches<br>> >><br>> >> 1. Using bond constraints and setting the bond angle to 180 and using<br>> >> a large bending constant to keep the molecule effectively rigid.<br>> >><br>> >> Here I got the following error message and the run stopped after 2 <br>> >> timesteps;<br>> >><br>> >> Step 1, time 0.001 (ps) LINCS WARNING<br>> >> relative constraint deviation after LINCS:<br>> >> rms 0.969889, max 10.563440 (between atoms 1072 and 1074)<br>> >> bonds that rotated more than 30 degrees:<br>> >> atom 1 atom 2 angle previous, current, constraint length<br>> >> 1072 1073 89.9 0.1210 0.6876 0.1209<br>> >> 1072 1074 92.4 0.1210 1.3978 0.1209<br>> >> 1075 1076 86.8 0.1209 0.1303 0.1209<br>> >> 1075 1077 86.8 0.1209 0.1307 0.1209<br>> >> Wrote pdb files with previous and current coordinates<br>> >><br>> >><br>> >> 2. Using constraints to keep the C-O bond fixed to 0.12 and another<br>> >> constraint between O and O of 0.24. I hoped this would keep the O-C-O<br>> >> linear.<br>> >><br>> >> Here, the md.log file contained the following;<br>> >><br>> >><br>> >> 6 constraints are involved in constraint triangles,<br>> >> will apply an additional matrix expansion of order 4 for couplings<br>> >> between constraints inside triangles<br>> >><br>> >> The job ran to completion but on viewing the trajectory, the CO2<br>> >> molecules are only approximately linear.<br>> >><br>> >> If anyone has an example file for a triatomic linear molecule like<br>> >> CO2, I?d be very grateful if you could post it so I could use it as a<br>> >> template,<br>> >><br>> >> I am using gromacs version 4.0.5<br>> >><br>> >> Thanks<br>> >><br>> >><br>> >><br>> >><br>> >><br>> >><br>> >><br>> >><br>> >><br>> >><br>> >><br>> >> Hi,<br>> >><br>> >> This question has been asked several times before on this list.<br>> >> The proper way to handle this in Gromacs is to introduce two dummy<br>> >> mass particles<br>> >> in such a way the total mass and the moment of inertia is identical to CO2.<br>> >> Then you can construct the positions of the three mass-less C and O <br>> >> atoms from<br>> >> the positions of the masses using virtual sites (see the manual).<br>> >><br>> >> Berk<br>> >><br>> >><br>> >><br>> >> Quoting Jennifer Williams <Jennifer.Williams@ed.ac.uk>:<br>> >><br>> >> > Hello users,<br>> >> ><br>> >> > I am using gromacs v 4.0.5<br>> >> ><br>> >> > What is the best way to treat linear triatomic molecules such as CO2.<br>> >> > Is there some kind of rigid algorthym as in DL_POLY?<br>> >> ><br>> >> > In the asbsence of a ?rigid? algorthym. I initially intended to ?fix?<br>> >> > the C=O bond length using constraints (as below) and somehow ?fix? the<br>> >> > O=C=O bond angle to 180.<br>> >> ><br>> >> ><br>> >> > [ atomtypes ]<br>> >> > ; type mass charge ptype c6 c12<br>> >> > C_CO2 12.01115 0.5406 A 0.0000 0.00000000<br>> >> > O_CO2 15.9994 -0.2703 A 0.29847 1.10765301<br>> >> ><br>> >> > [ moleculetype ]<br>> >> > ; name nrexcl<br>> >> > CO2 2<br>> >> ><br>> >> > [ atoms ]<br>> >> > ; nr type resnr residu atom cgnr charge mass<br>> >> > 1 C_CO2 1 CO2 C_CO2 1 0.5406 12.01115<br>> >> > 2 O_CO2 1 CO2 O_CO2 1 -0.2703 15.9994<br>> >> > 3 O_CO2 1 CO2 O_CO2 1 -0.2703 15.9994<br>> >> ><br>> >> > [ constraints ]<br>> >> > ; ai aj funct c0 c1<br>> >> > 1 2 1 0.12088<br>> >> > 1 3 1 0.12088<br>> >> ><br>> >> > [ angles ]<br>> >> > ; ai aj ak funct c0 c1<br>> >> > 2 1 3 1 180.00<br>> >> ><br>> >> > First question... How do I go about fixing this angle to 180degrees?<br>> >> > Should I just make c1 extremely large? Is there a more elegant way of<br>> >> > going about this? Someone also mentioned using a virtual site<br>> >> > representation for the C. Any comments on this welcome.<br>> >> ><br>> >> > I read some posts which imply that the angle bending code apparently can't<br>> >> > handle angles of exactly 180 degrees.<br>> >> ><br>> >> > http://www.mail-archive.com/gmx-users@gromacs.org/msg09979.html<br>> >> ><br>> >> > This was a while ago so perhaps this has been dealt with in the current<br>> >> > version of gromacs?<br>> >> ><br>> >> ><br>> >> > Thanks<br>> >><br>> >><br>> >><br>> >> --<br>> >> The University of Edinburgh is a charitable body, registered in<br>> >> Scotland, with registration number SC005336.<br>> >><br>> >><br>> >> _______________________________________________<br>> >> gmx-users mailing list gmx-users@gromacs.org<br>> >> http://lists.gromacs.org/mailman/listinfo/gmx-users<br>> >> Please search the archive at http://www.gromacs.org/search before posting!<br>> >> Please don't post (un)subscribe requests to the list. Use the<br>> >> www interface or send it to gmx-users-request@gromacs.org.<br>> >> Can't post? Read http://www.gromacs.org/mailing_lists/users.php<br>> ><br>> > _________________________________________________________________<br>> > See all the ways you can stay connected to friends and family<br>> > http://www.microsoft.com/windows/windowslive/default.aspx<br>> <br>> <br>> <br>> Dr. Jennifer Williams<br>> Institute for Materials and Processes<br>> School of Engineering<br>> University of Edinburgh<br>> Sanderson Building<br>> The King's Buildings<br>> Mayfield Road<br>> Edinburgh, EH9 3JL, United Kingdom<br>> Phone: ++44 (0)131 650 4 861<br>> <br>> <br>> -- <br>> The University of Edinburgh is a charitable body, registered in<br>> Scotland, with registration number SC005336.<br>> <br>> <br>> _______________________________________________<br>> gmx-users mailing list gmx-users@gromacs.org<br>> http://lists.gromacs.org/mailman/listinfo/gmx-users<br>> Please search the archive at http://www.gromacs.org/search before posting!<br>> Please don't post (un)subscribe requests to the list. Use the <br>> www interface or send it to gmx-users-request@gromacs.org.<br>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php<br><br /><hr />See all the ways you can stay connected <a href='http://www.microsoft.com/windows/windowslive/default.aspx' target='_new'>to friends and family</a></body>
</html>