[gmx-users] Questions about CO2.

[sujithkakkat .]

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.

       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