[gmx-users] Problems with CO2-water simulation

[Tsjerk Wassenaar]

Hi Sujith,

There is nothing special about positive potential energy. In this case it
explains why soft drinks gives bubbles: the CO2 wants to get out. In
addition, if it could, it would drive the system to HCO3-/H3O+, which is an
energetically more favorable way to store CO2 in water.

Cheers,

Tsjerk
On Dec 8, 2014 6:50 AM, "sujithkakkat ." <sujithks58 at gmail.com> wrote:

> 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.
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>