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

[Justin Lemkul]

On 10/9/18 6:02 AM, ARNAB MUKHERJEE wrote:
> 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.

The average needs to be interpreted in light of the fluctuations, which 
for pressure are usually massive.

http://www.gromacs.org/Documentation/Terminology/Pressure

> 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

Freezing in the absence of proper energygrp-excl leads to spurious 
contributions to the virial, as stated in the manual. Freezing is 
totally artificial, anyway, so you should carefully evaluate whether 
this is an appropriate approach.

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Assistant Professor
Virginia Tech Department of Biochemistry

303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com

==================================================