[gmx-users] Combining position restraints with Parrinello-Rahman pressure coupling leads to instabilities

[Naba]

Dear Gromacs users and developers,

I am using Gromacs 2018.2.
Following the membrane protein simulation tutorial, I am planning to run
long simulations of a tetramer that needs larger lipid bilayer than 128
lipids as described in the tutorial. So, I have replicated the
pre-equilibrated POPC bilayer of 128 lipids from lipidbook (
https://lipidbook.bioch.ox.ac.uk/lipid/) to get a larger bilayer of 512
lipids using genconf. Moreover, I have extended gromos54a7_lipid.ff
forcefield to Berger lipid parameters for long membrane protein
simulations. I have set 10 ns for NVT and 20 ns for NPT equilibration.
NVT simulation ran successfully, but when I reached to NPT equilibration,
grompp showed me a note like the following:
"You are combining position restraints with Parrinello-Rahman pressure
  coupling, which can lead to instabilities. If you really want to combine
  position restraints with pressure coupling, we suggest to use Berendsen
  pressure coupling instead."

When I used Berendsen pressure coupling, grompp terminated with a warning
as:
"Using Berendsen pressure coupling invalidates the true ensemble for the
  thermostat"

What to do with this? Can I proceed with the note using Parrinello-Rahman
pressure coupling? Please help.

Following is the mdp parameters I have used:

title       = NPT Equilibration for KALP15-DPPC
define      = -DPOSRES  ; position restrain the protein
; Run parameters
integrator  = md        ; leap-frog integrator
nsteps      = 10000000    ; 2 * 10000000 = 20000 ps (20 ns)
dt          = 0.002     ; 2 fs
cutoff-scheme = verlet
; Output control
nstxout     = 2500       ; save coordinates every 5 ps
nstvout     = 2500       ; save velocities every 5 ps
nstenergy   = 2500       ; save energies every 5 ps
nstlog      = 2500       ; update log file every 5 ps
; Bond parameters
continuation            = yes       ; Restarting after NVT
constraint_algorithm    = lincs     ; holonomic constraints
constraints             = all-bonds ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Neighborsearching
ns_type     = grid      ; search neighboring grid cels
nstlist     = 5         ; 10 fs
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)
; Electrostatics
coulombtype     = PME       ; Particle Mesh Ewald for long-range
electrostatics
pme_order       = 4         ; cubic interpolation
fourierspacing  = 0.16      ; grid spacing for FFT
; Temperature coupling is on
tcoupl      = Nose-Hoover                       ; More accurate thermostat
tc-grps     = Protein_POPC  Water_and_ions      ; two coupling groups -
more accurate
tau_t       = 0.5           0.5                 ; time constant, in ps
ref_t       = 300           300                 ; reference temperature,
one for each group, in K
; Pressure coupling is on
pcoupl      = Berendsen     ; Pressure coupling on in NPT
pcoupltype  = semiisotropic         ; uniform scaling of x-y box vectors,
independent z
tau_p       = 5.0                   ; time constant, in ps
ref_p       = 1.0   1.0             ; reference pressure, x-y, z (in bar)
compressibility = 4.5e-5    4.5e-5  ; isothermal compressibility, bar^-1
; Periodic boundary conditions
pbc         = xyz       ; 3-D PBC
; Dispersion correction
DispCorr    = EnerPres  ; account for cut-off vdW scheme
; Velocity generation
gen_vel     = no        ; Velocity generation is off
; COM motion removal
; These options remove motion of the protein/bilayer relative to the
solvent/ions
nstcomm         = 1
comm-mode       = Linear
comm-grps       = Protein_POPC Water_and_ions
; Scale COM of reference coordinates
refcoord_scaling = com
nstcalcenergy = 1
nhchainlength = 1


Thank in advance.

Regards,

Nabajyoti Goswami

Research Associate
Bioinformatics Infrastructure Facility
Department of Animal Biotechnology
College of Veterinary Science
Khanapara,Guwahati 781022
Assam, India