[gmx-users] Virtual sites and diffusion coefficients for CO2

[贺仲金]

Dear all,

    I am using rigid CO2 model with virtual sites. The bond lengthe and LJ parameters are taken from TraPPE force field (Vapor–liquid equilibria of mixtures containing alkanes, carbon dioxide, and nitrogen. JJ Potoff, JI Siepmann. AIChE journal 47 (7), 1676-1682, 2001. 469). I followed Justin A. Lemkul's tutorial to construct the virtual sites (http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/vsites/index.html). I find everything is OK except for the diffusion coeffients for CO2 are strange. The value with virtual sites (144.3087 (+/- 2.4934) 1e-5 cm^2/s) is half of that without virtual sites (282.5691 (+/- 20.2510) 1e-5 cm^2/s) . I just simulated 256 CO2 molecules at 298K and 55.2 bar with GROMACS 5.0.4. The itp files for CO2 with virtual sites is attched below. I found someone else also encounter such problems as shown in https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2015-February/095205.html and http://comments.gmane.org/gmane.science.biology.gromacs.user/23002 , and nobody solve out the origin of this issue.  Could anybody help me to check the itp files I build and find out the origin of this issue. Thanks in advance!



CO2_atomtypes.itp
[ atomtypes ]
;name  bond_type    mass    charge   ptype          sigma      epsilon
OZ        OZ      15.99940   0.0000    A           3.0500e-01  0.656806e0
CZ        CZ      12.01100   0.0000    A           2.8000e-01  0.224478e0
MCO       MCO       0.000    0.000     A              0.000      0.000


CO2_rigid.itp
[ moleculetype ]
; Name           nrexcl
CO2              2


[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB    chargeB
     1         OZ      1    CO2      O1      1   -0.35    0.000000
     2         CZ      1    CO2      C      1    0.70    0.000000
     3         OZ      1    CO2      O2      1   -0.35          0.000000
     4        MCO      1    CO2      M1      1     0.000          22.0049     
     5        MCO      1    CO2      M2      1     0.000          22.0049


[ constraints ]
; Overall structure is O1--M1---C---M2--O2 where "M" indicate mass centers
; There are no bonds in this system
; Instead, we fix the distance between the mass centers such that
; the virtual sites can be reconstructed
    4   5   1   0.197995


[ virtual_sites2 ]
; the M1--O1 distance is 0.1161 - (0.197995/2) = 0.0171025
; the M1--M2 distance is 0.197995
; therefore, the fraction of the distance O1M2 along the M1--M2 length is (0.197995+0.0171025)/0.197995 = 1.0863784
; thus placing the virtual O1 sites beyond the M1--M2 distance 
; site  ai(from)  aj  funct   a
     1   4   5      1   1.0863784  ; relative to mass center 4, extends beyond mass center 5 
     2   4   5      1   0.5000      ; right in the middle
     3   5   4      1   1.0863784   ; as in the case of site 1


[ exclusions ]
1 2 3
2 1 3
3 1 2




itp files without Virtual sites


CO2_atomtypes.itp
[ atomtypes ]
;name  bond_type    mass    charge   ptype          sigma      epsilon
OZ        OZ      16.0000   0.0000    A           3.0500e-01  0.656806e0
CZ        CZ      12.0000   0.0000    A           2.8000e-01  0.224478e0


CO2_rigid.itp
[ moleculetype ]
; Name           nrexcl
CO2              2


[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB    chargeB
     1         OZ      1    CO2      O1      1   -0.35    16.000000
     2         CZ      1    CO2      C      1    0.70    12.000000
     3         OZ      1    CO2      O2      1   -0.35          16.000000


[ Constraints ]
;  ai    aj funct  b0          
    1     2     1  0.1161      
    2     3     1  0.1161       
    1     3     1  0.2322




Best wishes!


Zhongjin HE