[gmx-users] using charmm36 force field with gromacs

[Justin Lemkul]

On 11/3/15 4:32 PM, Yoav Atsmon-Raz wrote:
> Hi there,
>
> I'm a noob with Gromacs 4.6.7 which is interested in doing all-atom simulations with the charmm36 force field for a protein-membrane system while applying the Verlet cutoff  scheme. Basically I have two questions regarding this matter:
> 1) Does using the CHARMM36 ff necessitates any special modifications to the mdp files (minimization,equilibration,prod. run)?
> 2) Besides the "cutoff-scheme=verlet" line in the mdp files, do I need to add anything else?
>

Several.  See below.  There used to be a page on the wiki with all this 
information but it somehow got deleted.  I'm trying to get it restored.

> This is my prod. run mdp file:
>
> ; Run parameters
> integrator    = md        ; leap-frog integrator
> nsteps        = 50        ; 2 * 500000 = 1000 ps (1 ns)
> dt            = 0.002        ; 2 fs
> ; Output control
> nstxout        = 10        ; save coordinates every 2 ps
> nstvout        = 10        ; save velocities every 2 ps
> nstxtcout    = 10        ; xtc compressed trajectory output every 2 ps
> nstenergy    = 10    ; save energies every 2 ps
> nstlog        = 10        ; update log file every 2 ps
> ; Bond parameters
> continuation    = yes            ; Restarting after NPT
> constraint_algorithm = lincs    ; holonomic constraints
> constraints    = all-bonds            ; all bonds (even heavy atom-H bonds) constrained

Use constraints = h-bonds

You'll get incorrect balance of 1-4 interactions if you constrain all bonds with 
CHARMM.

> 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)
> cutoff-scheme   = Verlet
> ; Electrostatics
> coulombtype    = PME        ; Particle Mesh Ewald for long-range electrostatics
> pme_order    = 4            ; cubic interpolation
> fourierspacing    = 0.16        ; grid spacing for FFT

The proper nonbonded settings are:

cutoff-scheme           = Verlet
vdwtype                 = cutoff
vdw-modifier            = force-switch  ; same as vfswitch
rlist                   = 1.2           ; will be tuned automatically
rvdw                    = 1.2
rvdw-switch             = 1.0
coulombtype             = PME
rcoulomb                = 1.2
pme_order               = 4
ewald_rtol              = 1e-5
optimize_fft            = yes
fourierspacing          = 0.12

Note that these are essential for getting membrane properties right.

> ; Temperature coupling is on
> tcoupl        = Nose-Hoover            ; More accurate thermostat
> tc-grps        = POPC_Protein Water            ; two coupling groups - more accurate
> tau_t        = 0.5    0.5                ; time constant, in ps
> ref_t        = 323     323                ; reference temperature, one for each group, in K
> ; Pressure coupling is on
> pcoupl        = Parrinello-Rahman        ; Pressure coupling on in NPT
> pcoupltype    = semiisotropic            ; uniform scaling of x-y box vectors, independent z
> tau_p        = 2.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

Do not apply dispersion correction for bilayers.  Do apply it for monolayers.

-Justin

> ; 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       = POPC_Protein Water
>
>
> Thanks in advance to anyone answering this!
>

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

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

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