[gmx-users] Average pressure in CG simulation of oil-water system at NAPzT ensemble

[Sudip Das]

Dear All,

I am running a coarse-grained martini (elastic network) simulation at
NAPzT ensemble for an oil-water system containing a protein at the
interface using GROMACS-5.1.4. The system has overall 27,000 atoms
including refined polarizable (v2.2refPOL) water molecules (The Journal of
Chemical Physics *146*, 054501 (2017); https://doi.org/10.1063/1.4974833).
Initially, the protein was kept in the water phase, but it got adsorbed at
the interface after around 0.15 microsecond.

Even though I have set the reference pressure to be 1 bar, the average
pressure from the simulations results to be around -34 bar, even after 1.5
microseconds run.

 Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Pressure                   -34.4124      0.018    35.1829  -0.102016  (bar)

I have also run the simulations by reducing the tau_p up to 4.0 ps, but I
don't find any improvement. Following same protocol with isotropic pressure
coupling, average pressure turns out to be around 11 bar.

Are these simulations with average pressure away from reference pressure
assumed to be correct to a reasonable extent?

The corresponding .mdp file is given below:

integrator = md
dt                    = 0.02
nsteps = 110000000

nstcomm         = 1
comm-mode = Linear
comm-grps = Oil   Protein_Water_Ion

nstxout            = 5000               ;[steps] freq to write coordinates
to trajectory
nstvout            = 5000               ;[steps] freq to write velocities
to trajectory
;nstfout             = 5000               ;[steps] freq to write forces to
trajectory
nstlog              = 5000               ;[steps] freq to write energies to
log file
nstenergy        = 5000               ;[steps] freq to write energies to
energy file
energygrps = system            ;Which energy group(s) to write to disk

cutoff-scheme     = Verlet
nstlist                  = 20
ns_type               = grid
pbc                      = xyz
verlet-buffer-tolerance  = 0.005

coulombtype            = PME        ; for normal MARTINI water, use
'reaction-field'
rcoulomb                 = 1.1
epsilon_r                 = 2.5       ; 15 (with normal MARTINI water)
;epsilon_rf               = 0
vdw_type                 = cutoff
vdw-modifier            = Potential-shift-verlet
rvdw                         = 1.1
DispCorr                  = EnerPres   ; account for cut-off vdW scheme

tcoupl            =  v-rescale
tc-grps          =  Oil    Protein_Water_Ion
tau-t              =  1.0   1.0
ref-t               =  300  300

Pcoupl               = parrinello-rahman
Pcoupltype        = semiisotropic
tau_p                 = 12.0                ; PR barostat is more stable
with larger tau-p, DdJ, 20130422
ref_p                  = 1.0   1.0          ; reference pressure, x-y, z
(in bar)
compressibility   = 0      3e-4       ; isothermal compressibility, bar^-1

gen_vel          = no

constraints                = none
constraint_algorithm = Lincs
lincs_order                = 4
lincs_warnangle        = 30


Thanks in advance.

Best regards,
Sudip