[gmx-users] Cyclic Peptide in Cyclohexane - LINCS warnings

[Billy Williams-Noonan]

Hi GROMACS Users,

Am running a simulation of a cyclic peptide with some unusual residues in a
box of cyclohexane.

The parameters for this peptide are sources from the ATB, as are the
parameters for the cyclohexane solvent.  Am using the Gromos 54a7
parameters with GROMACS-2016.

My simulation keeps crashing after repeated LINCS warnings and after
re-running the simulations with just a box of cyclohexane, as well as
re-running with a simulation of the cyclic peptide in vacuum, I've found
that it is the cyclohexane, specifically, that is the problem.  By looking
at the step files, you can see that the atoms of each solvent molecule come
into close-contact with each other, creating a high amount of force to be
experienced by the relevant atoms, and thereby causing the relevant
molecules to 'blow up'.

I was wondering why this would be the case.

Here are my *.mdp *details.

title           = test_cylohexane
; Run parameters
integrator      = md            ; leap-frog integrator
nsteps          = 50000000              ;  = 2 * 50000 ps = 100 ns
dt                  = 0.002             ; 1 fs
; Output control
nstxout         = 5000          ; save coordinates every 1.0 ps
nstvout         = 5000          ; save velocities every 1.0 ps
nstenergy       = 5000          ; save energies every 1.0 ps
nstlog          = 5000          ; update log file every 1.0 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
cutoff-scheme   = Verlet
ns_type             = grid             ; search neighboring grid cells
nstlist             = 40                 ; 80 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         = half_A half_B         ; two coupling groups - more
accurate
tau_t           = 0.1     0.1           ; time constant, in ps
ref_t           = 300     300           ; 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                   = 2.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

Here are is the topology for cyclohexane.

[ moleculetype ]
; Name   nrexcl
_VMQ     3
[ atoms ]
;  nr  type  resnr  resid  atom  cgnr  charge    mass    total_charge
    1  CH2r    1    _VMQ     C1    1    0.000  14.0270
    2  CH2r    1    _VMQ     C6    1    0.000  14.0270     ;  0.000
    3  CH2r    1    _VMQ     C2    2    0.000  14.0270
    4  CH2r    1    _VMQ     C3    2    0.000  14.0270     ;  0.000
    5  CH2r    1    _VMQ     C4    3    0.000  14.0270
    6  CH2r    1    _VMQ     C5    3    0.000  14.0270     ;  0.000
; total charge of the molecule:   0.000
[ bonds ]
;  ai   aj  funct   c0         c1
    1    2    2   0.1530   7.1500e+06
    1    3    2   0.1530   7.1500e+06
    2    6    2   0.1530   7.1500e+06
    3    4    2   0.1530   7.1500e+06
    4    5    2   0.1530   7.1500e+06
    5    6    2   0.1530   7.1500e+06
[ pairs ]
;  ai   aj  funct  ;  all 1-4 pairs but the ones excluded in GROMOS itp
    1    5    1
    2    4    1
    3    6    1
[ angles ]
;  ai   aj   ak  funct   angle     fc
    2    1    3    2    111.00   530.00
    1    2    6    2    111.00   530.00
    1    3    4    2    111.00   530.00
    3    4    5    2    111.00   530.00
    4    5    6    2    111.00   530.00
    2    6    5    2    111.00   530.00
[ dihedrals ]
; GROMOS improper dihedrals
;  ai   aj   ak   al  funct   angle     fc
[ dihedrals ]
;  ai   aj   ak   al  funct    ph0      cp     mult
    3    1    2    6    1      0.00     5.92    3
    2    1    3    4    1      0.00     5.92    3
    1    2    6    5    1      0.00     5.92    3
    1    3    4    5    1      0.00     5.92    3
    3    4    5    6    1      0.00     5.92    3
    4    5    6    2    1      0.00     5.92    3
[ exclusions ]
;  ai   aj  funct  ;  GROMOS 1-4 exclusions


The periodic box dimensions were set up so that the density of the system
was 0.779 g/mL, which is the density of cyclohexane from experiment.

I should probably add that there was one simulation where I could not see
two cyclohexane molecules obviously 'touching' (this was run in triplicate
for 100 ns) - but two molecules of cyclohexane 'blew up' regardless.
Further, reducing the time-step to 1 fs, or using a Berendsen barostat
eliminates the LINCS warnings an prevents crashing, but I'm unsure if I'm
actually fixing the problem by doing that.  Further, my understanding is
that Berendsen gives the wrong ensemble pressure, and I'd rather avoid
using it, just as I'd like my simulations to run more quickly and have a 2
fs time-step.

Is there anything wrong with what I'm doing or anything I should add to the
*.mdp* file to make it work?

Cheers,

Billy




-- 
Billy Noonan*    |    *PhD Student    *|*    Bsci ( *Adv* ), IA Hon

*LinkedIn Profile
<http://www.linkedin.com/profile/preview?locale=en_US&trk=prof-0-sb-preview-primary-button>
**|*   +61420 382 557

Monash Institute for Pharmaceutical Sciences ( *MIPS* )
Royal Parade, Parkville, 3052