[gmx-users] Unexpected lincs errors for OH groups in protein

[Erik Marklund]

Dear GMX-users,

I am initiating a series of simulations of a dimeric protein in water, i.e. something that should be fairly trivial. I use gromacs 4.6.5, vsites for the hydrogens, amber99sb-ildn for the protein, and tip4p-Ew for the water. With a time step of 0.004 I keep getting occasional lincs errors for OH groups that swing around to violently, to my surprise. There are no clashes in the structure that I can see and the OH groups are mostly at the protein-water interface and as such quite unconstrained by the surroundings. I am using a similar setup to what I have successfully been using with gromacs 4.5.X on numerous occasions, but non-bonded mdp options have changed so I suspect that's where the problem is rooted, but so far I have been unable to overcome these errors. EM proceeded just fine, a 100 ps NVT run with position restraints and a 0.002 fs time step, and a 1 ns NVT run with time step 0.002 fs preceded this NVT run  Is there anything obviously wrong with the following combination of mdp parameters?

integrator               = md
; Start time and timestep in ps
tinit                    = 0
dt                       = 0.004
nsteps                   = 250000
…

cutoff-scheme            = group
; nblist update frequency
nstlist                  = 10
…

rlist                    = 1.0
; long-range cut-off for switched potentials
rlistlong                = -1
nstcalclr                = -1

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype              = PME
coulomb-modifier         = Potential-shift-Verlet ; <- moot
rcoulomb-switch          = 0
rcoulomb                 = 1.0
; Relative dielectric constant for the medium and the reaction field
epsilon-r                = 1
epsilon-rf               = 0
; Method for doing Van der Waals
vdw-type                 = Cut-off
vdw-modifier             = Potential-shift
; cut-off lengths       
rvdw-switch              = 0
rvdw                     = 0.8
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                 = EnerPres
…

; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling  
tcoupl                   = v-rescale
nsttcouple               = -1
nh-chain-length          = 10
print-nose-hoover-chain-variables = no
; Groups to couple separately
tc-grps                  = Protein      Water_and_ions
; Time constant (ps) and reference temperature (K)
tau-t                    = 0.2  0.2
ref-t                    = 300  300
; pressure coupling     
Pcoupl                   = Berendsen
pcoupltype               = Isotropic
nstpcouple               = -1
…

; GENERATE VELOCITIES FOR STARTUP RUN
gen-vel                  = no
gen-temp                 = 300
gen-seed                 = 173529

; OPTIONS FOR BONDS    
constraints              = all-bonds
; Type of constraint algorithm
constraint-algorithm     = Lincs
; Do not constrain the start configuration
continuation             = no
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR                = no
; Relative tolerance of shake
shake-tol                = 0.0001
; Highest order in the expansion of the constraint coupling matrix
lincs-order              = 6
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter               = 1
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle          = 30 
...


I tried various alterations to this setup, e.g. rvdw = 1.0, but still no success. Any help coming to grips with the new non-bonded stuff is much welcome.

Kind regards,

Erik Marklund, PhD
Postdoctoral Research Associate

Department of Chemistry
Physical & Theoretical Chemistry Laboratory
University of Oxford
South Parks Road
Oxford
OX1 3QZ