[gmx-users] Problem during equilibration using semi-isotropic pressure coupling

[ARNAB MUKHERJEE]

Hi everyone,

I have a system of infinite CG (Martini) DNA (length of the DNA same as
length of the periodic box) protein system. I am trying to equilibrate this
system to 1 bar pressure. Since my DNA is same as the size of the box, I
don't want the box to fluctuate along Z direction. Hence I am using semi
isotropic pressure coupling, with the box fluctuating only along XY.
I have built 2 systems. One where the DNA is completely frozen, and the 2nd
one only DNA end points are frozen. For both the systems I ran
equilibration NPT run for over 20 ns, and still the average pressure
doesn't show 1 bar, which is the target pressure. For the only end points
frozen DNA protein system, the average pressure shows 1.5 bar, and for the
completely frozen DNA system, the average pressure shows around 5.5 bar,
after running for over 20 ns for both the systems.
Here is how my .mdp file looks :

 title       = NVT equilibration with position restraint on all solute
(topology modified)
; Run parameters
integrator  = sd        ; leap-frog integrator
tinit       = 0
init-step   = 16000000
nsteps      = 12000000
dt          = 0.001     ; 1 fs
; Output control
nstxout     = 0 ; save coordinates every 10 ps
nstvout     = 0 ; save velocities every 10 ps
nstcalcenergy   = 50
nstenergy   = 1000  ; save energies every 1 ps
nstxtcout       = 2500
;nstxout-compressed  = 5000   ; save compressed coordinates every 1.0 ps
                             ; nstxout-compressed replaces nstxtcout
;compressed-x-grps  = System  ; replaces xtc-grps
nstlog      = 1000  ; update log file every 1 ps
; Bond parameters
continuation    = no   ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = none          ; all bonds (even heavy atom-H bonds)
constrained
;lincs_iter = 2                 ; accuracy of LINCS
lincs_order = 4                 ; also related to accuracy
epsilon_r       = 15
; Neighborsearching
cutoff-scheme   = Verlet
ns_type     = grid      ; search neighboring grid cels
nstlist     = 10            ; 20 fs
rvdw_switch     = 1.0
rlist       = 1.2       ; short-range neighborlist cutoff (in nm)
rcoulomb    = 1.2       ; short-range electrostatic cutoff (in nm)
rvdw        = 1.2       ; short-range van der Waals cutoff (in nm)
vdwtype         = Cut-off       ; Twin range cut-offs rvdw >= rlist
;vdw-modifier    = Force-switch
;Electrostatics
coulombtype = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order   = 4         ; cubic interpolation
fourierspacing  = 0.12      ; grid spacing for FFT
; Temperature coupling is on
tcoupl          = v-rescale
tc_grps         = System
tau_t           = 2.0
ld-seed  =  -1

;energygrps = DNA W_ION_Protein
;energygrp-excl = DNA DNA
freezegrps = DNA
freezedim = Y N Y

; Pressure coupling is off
;pcoupl     = no        ; no pressure coupling in NVT
Pcoupl     = parrinello-rahman
Pcoupltype  = semiisotropic
;tau_p       = 5.0
tau_p       = 2.5
compressibility = 3e-4 0
ref_p       = 1.0 1.0
; Periodic boundary conditions
pbc     = xyz       ; 3-D PBC
; Dispersion correctiion
DispCorr    = no        ; account for cut-off vdW scheme
; Velocity generation
gen_vel     = no        ; assign velocities from Maxwell distribution
gen_temp    = 300       ; temperature for Maxwell distribution
gen_seed    = -1        ; generate a random seed
; COM motion removal
; These options remove motion of the protein/bilayer relative to the
solvent/ions
nstcomm     = 50
comm-mode   = Linear
comm-grps       = System
;
;refcoord_scaling = com
refcoord_scaling = all

As you see I used Parrinello Rahman Pressure coupling. Initially I had used
tau_p = 5.0 ps, and then I reduced it to 2.5 ps but there is still not much
change. I also did equilibration run for Finite DNA system with isotropic
pressure coupling and using tau_p = 5.0 ps, I got an average pressure of 1
bar just after 5ns NPT run. While calculating the average pressure for both
these systems, the RMSD shows around 10 bar.

For the infinite DNA (completely frozen DNA system) I have tried using
Berensden thermostat with tau_p = 1ps, and after 4 ns run, I have average
pressure around 7 bar. I am continuing this run to check if it shows better
results.

I do not understand what is problem occuring while equlibrating for the
infinite DNA systems using semi-isotropic pressure coupling? Why it cannot
reach the target pressure?

I would highly appreciate any help.

Thank you very much,

Regards,

Arnab