[gmx-users] Problems with CO2-water simulation

[sujithkakkat .]

Dear all,

       I am trying to simulate a water-carbondioxide system at 250K and
30bar. The system consists of 85 CO2 and 1204 water molecules in a 3x3x4 nm
box.

       I used EPM2 and TIP4P potentials for CO2 and water following  *Energy
Environ. Sci., 2012, 5, 7033-7041* where the authors used the same.\

       I am facing problems with the energy potential energy values which
are very high and positive ( of the order 10^7).

       I tried the same simulation with the OPLSAA topology of CO2 which I
obtained from Dr. Justin Lemkul's tutorial on virtual site with TIP4P
water. The problem still remains, potential energy values being highly
positive.

       I am totally stuck at this point and is unable to find how to tackle
the problem.

       Firstly, I am wondering whether the positive potential energy is
something to be worried about. If so, it would be great if you could
provide me some direction in which I should think, or suggest me something
to read about on this issue.

       Given below is the parameter file I used (I have redueced the step
size due to a avoid persistent blowing up problem);

 title                   =  co2--in-water

; Run parameters
integrator            = md
nsteps                = 2500000
dt                       = 0.0002

; Output control
nstxout               = 5000
nstvout               = 5000
nstenergy            = 5000
nstlog                  = 5000
nstxtcout             = 5000

; Bond parameters
continuation               = no
constraint_algorithm    = lincs
constraints                  = all-bonds
lincs_iter                     = 2
lincs_order                   = 4

; Neighborsearching
cutoff-scheme              = Verlet
ns_type                        = grid
nstlist                           = 10
rcoulomb                       = 1.2
rvdw                              = 1.2

; Electrostatics
coulombtype                  = PME
pme_order                      = 4
fourierspacing                 = 0.16

; Temperature coupling is on
tcoupl                             = Nose-Hoover
tc-grps                            = CO2   SOL
tau_t                               = 0.4   0.4
ref_t                                = 250   250

; Pressure coupling is on
pcoupl                              = Parrinello-Rahman
pcoupltype                        = isotropic
tau_p                                = 2.0
ref_p                                 = 30.0
compressibility                   = 4.5e-5
refcoord_scaling                 = com

; Periodic boundary conditions
pbc                                    = xyz

; Dispersion correction
DispCorr                             = EnerPres

; Velocity generation
gen_vel                               = no
gen_temp                            = 250
gen_seed                            = -1



The EPM2 topology file which I made is given below;

[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



Regards,

Sujith.