[gmx-users] Questions about CO2.

[David van der Spoel]

On 2014-12-25 08:13, sujithkakkat . wrote:
> Dear all,
>
>        I went through the mails about the use of virtual sites for
> representing linear molecules like CO2. I also read that GROMACS cannot
> handle angle bending potentials at 180 degrees in an appropriate manner.
There is an (undocumented) linear angle potential in gromacs. Check 
source code and manual.

>
>         My system consists of EPM2-CO2 dissollved in TIP4- H2O.
>
>          I made topology files for CO2 with virtual sites and also made
> another topology without the use of virtual sites providing the equilibrium
> angle (180) and the bending force constant. Both these are given at the end.
>
>           I am trying to use the EPM2 model for CO2 reported in* J. Phys.
> Chem. 1995,99, 12021-12024* by *Jonathan G. Harris and Kwong H. Yung* under
> the title "Carbon Dioxide's Liquid-Vapor Coexistence Curve and Critical
> Properties As Predicted by a Simple Molecular Model".
>
> I have the following questions;
>
>   (i) I am using gromacs 4.6.5. Is the problem related to treating angle
> bending potentials for linear molecules solved in later versions of gromacs?
>
>   (ii) The authors Harris and Yung  stress that their model is flexible in
> terms of bending in the abstract of their article. So if I build this model
> using virtual sites then I believe I am neglecting the flexibility. Also
> while using virtual sites the model looks more like a dimer with two point
> masses, instead of a a three point EPM2 model.  So can I call this an EPM2
> model any more.? Do you think introducing virtual sites affects the results
> which are expected from a simple looking EPM2 model, with three
> charged-Lennard-Jones point masses with a flexible angle.
>
> I appreciate (and often surprised by) the fact that a simple model like
> EPM2 which can be specified by a total of just nine parameters can
> accurately predict the properties of CO2 at supercritical conditions which
> otherwise may a be a very intricate condensed system to study. This prompts
> me to think that any small change in the model can impart an error that can
> get magnified for a system of thousands of molecules.
>
>   (iii) I found that I am getting highly positive energy values when using
> the topology built using virtual sites. To add to this I am forced to use
> very short time steps (0.2 fs) for simulations, since with 2fs, system was
> unstable which I believe has to do with the positive energy and high force
> values. Whereas, using the topology built in the usual way without virtual
> sites, I get negative energy values and simulations could easily be
> performed using a 2fs time step (system may be more stable) . I checked the
> angles in this case at the final state of the system, and found them lying
> between 175 to 180 degrees,  which I am not sure if it is all right.
>        Is this slight deviation from 180 degree caused by the failure of
> GROMACS to maintain linearity? Or, is it perfectly normal since the model
> itself is flexible.?
>
>
> Please share your views.
>
> Regards,
>
> Sujith.
>
>
>
> *TOPOLOGY USING VIRTUAL SITES:*
>
>
>   [atomtypes]
>
> ; name   mass   charge   ptype    sigma   epsilon
>
>     D      22.0049   0.0000  A         0.0000   0.0000
>    CA     0.0000    0.6512   A         0.2757   0.2339
>    CO     0.0000   -0.3256   A         0.3033   0.6695
>
>   [moleculetype]
>
> ; name  nrexcl
>     CO2    2
>
>   [atoms]
>
> ; nr type resnr residue atom cgnr charge mass
>    1   D      1     CO2     D1      1    0.0000   22.0049
>    2   D      1     CO2     D2      1    0.0000   22.0049
>    3   CA    1     CO2    CA      1    0.6512   0.0000
>    4   CO    1     CO2    OC1    1   -0.3256   0.0000
>    5   CO    1     CO2    OC2    1   -0.3256   0.0000
>
>   [constraints]
> ; i   j   funct   doc
>    1  2   1       0.195948
>
>   [virtual_sites2]
> ; i  j   k funct  a
>   3  1  2   1     0.5
>   4  1  2   1     1.08638006
>   5  2  1   1     1.08638006
>
>
>   *TOPOLOGY WITHOUT VIRTUAL SITES:*
>
>   [ defaults ]
> ; nbfunc   comb-rule   gen-pairs    fudgeLJ   fudgeQQ
>      1           2              no               1.0          1.0
>
>   [atomtypes]
>
> ;name   mass         charge   ptype   sigma     epsilon
>   CA     12.01100    0.6512      A      0.2757    0.2339
>                                         CO     15.99940   -0.3256      A
>    0.3033    0.6695
>
>   [ nonbond_params ]
>
> ; i      j      funct   sigma    epsilon
>
>   CA  CO    1        0.2892   0.3955
>
>   [moleculetype]
>
> ; name nrexcl
>
>     CO2    2
>
>   [atoms]
>
> ; nr  type resnr  residue  atom cgnr   charge  mass
>    1   CA    1       CO2     CA     1      0.6512  12.0110
>    2   CO    1       CO2     OC1   1     -0.3256  15.9994
>    3   CO    1       CO2     OC2   1     -0.3256  15.9994
>
>   [angles]
>
> ; i j k funct ao
> ak
> 3   1  2  1  180  1236
>
>
>   [constraints]
>
> ;  i   j   funct   doc
>
>    1   2   1       0.1149
>    1   3   1       0.1149
>
>   [exclusions]
>
> 1  2  3
> 2  1  3
> 3  1  2
>


-- 
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
spoel at xray.bmc.uu.se    http://folding.bmc.uu.se