[gmx-users] Unraveling Helix: Simulation Data not in Agreement with Experiment

[Matthew Fisher]

Dear all,


I'm simulating P450 proteins currently and having issues with my simulation data not agreeing with experimentally obtained data. A helix, located immediately above the heme on the distal side, is either unravelling or distorting significantly during production MD runs, whereas crystal structures of identical/similar mutants indicate no such distortion from the WT.


Briefly, I'm using GROMACS 2018.1 and the AMBER99SB-ILDN force field. I'm doing a 2 stage EM: Steepest Descent with a criterion of 500 kJ / mol nm followed by conjugate gradient (100 kJ / mol nm) and then position restrained NVT equilibriation with the v-rescale thermostat (200ps)  followed by position restrained NPT with Berendsen barostat (I'm then changing to Parrinello-Rahman for the production runs.)


I've tried altering time constants to no avail and was wondering if anyone had any ideas what might be causing this / places to focus my attention?



I've included my .mdp parameters below for the production run.


Many thanks,

Matthew Fisher


; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 50000000  ; 2 * 50000000 fs = 100000 ps (100 ns)
dt                      = 0.002     ; 2 fs
; Output control
nstxout                 = 5000         ; save coordinates every 10.0 ps
nstvout                 = 5000         ; save velocities every 10.0 ps
nstfout                 = 5000         ; save forces every 10.0 ps
nstenergy               = 5000         ; save energies every 10.0 ps
nstlog                  = 5000         ; update log file every 10.0 ps
nstxout-compressed      = 5000         ; save compressed coordinates every 10.0 ps
compressed-x-grps       = System    ; save the whole system
; Bond parameters
continuation            = yes       ; Restarting after NPT
constraint_algorithm    = lincs     ; holonomic constraints
constraints             = all-bonds ; all bonds are constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Neighborsearching
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; 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                  = V-rescale             ; modified Berendsen thermostat
tc-grps                 = Protein Non-Protein   ; two coupling groups - more accurate
tau_t                   = 0.1     0.1           ; time constant, in ps
ref_t                   = 298     298           ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl                  = Parrinello-Rahman     ; Pressure coupling on in NPT
pcoupltype              = isotropic             ; uniform scaling of box vectors
tau_p                   = 5.0                   ; time constant, in ps
ref_p                   = 1.0                   ; reference pressure, in bar
compressibility         = 4.5e-5                ; isothermal compressibility of water, 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