[gmx-users] Unfolding event after restart

[Daniel Luesebrink]

Hello,

In umbrella sampling simulations I have observed unusual unfolding events,
which appear after restarting simulations.

I have simulated a homo-dimer with about 150 residues for each monomer. The
system has been set up with TIP3P explicit solvent and the CHARMM22* force
field.
In the umbrellas I use the pull code (see mdp-file below) to confine the
monomers. I also use the enforced rotation module keep the orientation of
each monomer separately fixed.

The problem arises when I restart simulations. I restart them by using the
check point file and either the original tpr file or a new tpr file for
extending the simulation time (see below). In about 10 of my 200
trajectories the protein had unfolded right after restarting.
The unfolding looks artificial and it happens within a few ps. Before the
protein had been simulated for 20ns without any problems and the unfolding
does not seem to have a physical reason. In some case almost the whole
contour length of the backbone is unfolded. Also the unfolded
configurations sometimes seem to be squeezed into a plane and happen
strangly symmetric with respect to the two monomers. In other cases the
unfolding is not as severe and only local, but still seems artificial from
visual observations of the protein structure.

For restarting I use:
mdrun -deffnm umbrella -pf pullf-umbrella.xvg -px pullx-umbrella.xvg -cpi
umbrella.cpt
(The tpr file is called umbrella.tpr)
or:
mdrun -s umbrella_ext.tpr -deffnm umbrella -pf pullf-umbrella.xvg -px
pullx-umbrella.xvg -cpi umbrella.cpt
The umbrella_ext.tpr is created by:
tpbconv -s umbrella.tpr -extend 10000 -o umbrella_ext.tpr

I was using gromacs-4.6.2 and 4.6.5. My system has about 130k atoms and was
using between 96 and 144 CPU's.

I was wondering if the issue might be related to the enforced rotation
module, since I have only experienced the unfolding events in the case
where I was using this module.

Please find the content of my mdp-file below. I would be happy to know, if
anyone help me to find out what is the issue.

Cheers,
Daniel

MDP file:
title       = Umbrella sampling explicit solvent
; Run parameters
integrator  = md
dt          = 0.002
tinit       = 0
nsteps      = 10000000   ; 20 ns
comm-mode   = None
; Output parameters
nstxout     = 50000     ; every 100 ps
nstvout     = 50000
nstfout     = 50000
nstxtcout   = 5000      ; every 10 ps
nstenergy   = 5000
; Bond parameters
constraint_algorithm    = lincs
constraints             = all-bonds
continuation            = yes
; Single-range cutoff scheme
nstlist     = 5
ns_type     = grid
rlist       = 1.4
rcoulomb    = 1.4
rvdw        = 1.4
; PME electrostatics parameters
coulombtype     = PME
fourierspacing  = 0.12
fourier_nx      = 0
fourier_ny      = 0
fourier_nz      = 0
pme_order       = 4
ewald_rtol      = 1e-5
optimize_fft    = yes
; Berendsen temperature coupling is on in two groups
Tcoupl      = Nose-Hoover
tc_grps     = Protein   Non-Protein
tau_t       = 0.5       0.5
ref_t       = 300       300
; Pressure coupling is on
Pcoupl          = Parrinello-Rahman
pcoupltype      = isotropic
tau_p           = 1.0
compressibility = 4.5e-5
ref_p           = 1.0
refcoord_scaling = com
; Generate velocities is off
gen_vel     = no
; Periodic boundary conditions are on in all directions
pbc     = xyz
; Long-range dispersion correction
DispCorr    = EnerPres
; Pull code
pull            = umbrella
pull_geometry   = position
pull_dim        = Y Y Y         ; umbrella potential in all three dimensions
pull_start      = yes           ; define initial COM distance > 0
pull_ngroups    = 1
pull_group0     = freeze
pull_group1     = pull
pull_vec1       = 0 0 1         ; pull along z
pull_rate1      = 0.0           ; nm per ps
pull_k1         = 10000         ; kJ mol^-1 nm^-2
pull-nstxout    = 10
pull-nstfout    = 10
; ENFORCED ROTATION
; Enforced rotation: No or Yes
rotation                 = Yes
; Output frequency for angle, torque and rotation potential energy for the
whole group
rot-nstrout              = 10
; Output frequency for per-slab data (angles, torques and slab centers)
rot-nstsout              = 10
; Number of rotation groups
rot-ngroups              = 2
; Rotation group name
rot-group0               = freeze_Ca
rot-group1               = pull_Ca
; Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2,
rm2-pf, flex, flex-t, flex2, flex2-t
rot-type0                = rm2-pf
rot-type1                = rm2-pf
; Use mass-weighting of the rotation group positions
rot-massw0               = no ; only Ca, all masses equal
rot-massw1               = no
; Rotation vector, will get normalized
rot-vec0                 = 1 0 0 ; specify arbitrary vector, rot mat is id
regardless, because rot-rate = 0
rot-vec1                 = 1 0 0
; Rotation rate [degree/ps] and force constant [kJ/(mol*nm^2)]
rot-rate0                = 0.0
rot-k0                   = 1000.0
rot-rate1                = 0.0
rot-k1                   = 1000.0
; Fitting method to determine angle of rotation group (rmsd, norm, or
potential)
rot-fit-method0          = rmsd
rot-fit-method1          = rmsd
; Value of additive constant epsilon' [nm^2] for rm2* and flex2* potentials
rot-eps0                 = 0.01
rot-eps1                 = 0.01