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

[Justin Lemkul]

On 8/20/18 6:52 AM, Naba wrote:
> 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."

I think this note is a bit imprecise. grompp is guessing that the 
combination of Parrinello-Rahman and restraints implies equilibration, 
and it should really say so, or otherwise do a more rigorous check to 
see if velocities are being generated. In principle, there shouldn't be 
anything wrong with this approach, but Parrinello-Rahman *usually* isn't 
the best choice for equilibration.

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

This only matters for production runs. The Berendsen method should never 
be used for actual data collection. For equilibration, it's perfectly 
fine because you're going to throw this time out, anyway.

-Justin

> 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

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

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

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