[gmx-users] Fwd: Help with Gromacs 5.1 mdp options for CHARMM27 force field

Justin Lemkul jalemkul at vt.edu
Fri Dec 4 13:34:57 CET 2015



On 12/3/15 6:55 AM, Rakesh Ramachandran wrote:
> Hi Justin,
>
>     Thanks for the reply. I have doubts regarding pdb2gmx options are
> the chargegrp and cmap options required while generating topology file when
> using CHARMM36 force field and do these options have any effect.
>

CHARMM does not use charge groups; every atom is its own group.  Moreover, 
charge groups are ignored when using the Verlet cutoff scheme (which you need to 
be, to use the settings I provided).  The [ cmap ] directives are required if 
you want to actually use CMAP for the backbone, which is an integral part of the 
force field.

>     Moreover I saw in your recent paper "Phosphorylation of PPARγ Affects
> the Collective Motions of the PPARγ-RXRα-DNA Complex", you had used H++ web
> server to assign protonation state. The server generates only PDB file and
> AMBER topology file so how can I use this in Gromacs. I would be grateful

Well, we used an AMBER force field for those simulations, so using an 
AMBER-prepared coordinate file is seamless.  The important part is the set of 
pKa values and subsequent protonation states.  You can set these using pdb2gmx 
interactively (read the help information).

> if you could share the code or script for PCA based energy landscape
> generation shown in your paper as I am also interested in the collective
> dynamics.
>

If you wish to discuss my work specifically, please contact me offline rather 
than posting to the mailing list and boring everyone else :)  The PCA is a 
straightforward application of gmx covar and gmx anaeig.  There's nothing 
special there.

>     Also I want to retain crystallographic water and ions, but I see the box
> size increases and also see that the naming is different HOH for
> crystallographic water and SOL for TIP3P.  Will this have any effect or is
> it always better to remove any crystallographic water or ions before
> simulating.
>

Depends on whether or not those waters are structurally relevant.  Usually there 
is no harm in keeping them, sometimes removing something important in a 
binding/active site is inappropriate.

-Justin

>
> On 1 December 2015 at 18:14, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>>
>> On 11/30/15 3:17 PM, Rakesh Ramachandran wrote:
>>
>>> Hi Justin,
>>>
>>>       Thanks for the reply. I have few more doubts. I had mentioned
>>> CHARMM27
>>> just for consistency with Gromacs naming. In the website it mentions that
>>> those settings are for CHARMM36 and I read at lot of places that this
>>> force
>>> field is good mainly for membrane proteins, so is it ok for me to use it
>>>
>>
>> I wouldn't limit the applicability of CHARMM36 to simply membrane
>> proteins. It contains updated CMAP parameters that are superior to the old
>> C22/CMAP force field and thus C36 is recommended in general over C22/CMAP.
>>
>> with CHARMM22/CMAP. Moreover I am using Verlet cutoff with GPU, so are
>>> these options still applicable or I have to use vdw-modifier:
>>> potential-shift-verlet.
>>>
>>>
>> The settings I linked have been thoroughly tested and validated. I do not
>> know the accuracy of any other settings.
>>
>> I also noticed that the Dispcorr is set to 'no', is this ok. Also I
>>>
>>
>> With the CHARMM force field, dispersion correction is only applied in the
>> case of lipid monolayers, per the link I sent.
>>
>> prepared my input files with CHARMM-GUI and I see tau_t - 1.0 and tau_p 5.0
>>> whereas your tutorial has different values. I am confused as to how to
>>> determine the optimal set of parameters and are all these parameters
>>> optimal for CHARMM36 only.
>>>
>>>
>> There will be minor differences in ensembles, but average properties
>> should agree.  See, for instance
>> http://pubs.acs.org/doi/abs/10.1021/acs.jctc.5b00935.
>>
>> -Justin
>>
>>
>> I would be grateful if you could let me know what kind of tests or values
>>> one needs to check to know if the parameters are appropriate for running
>>> md
>>> with a protein in a particular force field or point me to the relevant
>>> literature.
>>>
>>> On 30 November 2015 at 06:10, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>
>>>
>>>>
>>>> On 11/29/15 7:35 PM, Rakesh Ramachandran wrote:
>>>>
>>>> Dear all,
>>>>>
>>>>>         I am using Gromacs 5.1 with CHARMM27 force field for protein
>>>>> simulation and using the following mdp options. I am really confused
>>>>> whether to use PME-switch option or only PME and what is the basic
>>>>> difference. Moreover for CHARMM I see switching needs to be performed,
>>>>> but
>>>>> with Verlet cutoff should I use Potential-switch, Force-switch
>>>>> or potential-shift-verlet. Also let me know if any other options need to
>>>>> be
>>>>> changed.
>>>>>
>>>>>
>>>>> ; NEIGHBORSEARCHING PARAMETERS
>>>>> ; nblist update frequency
>>>>> nstlist = 40
>>>>>
>>>>> ; ns algorithm (simple or grid)
>>>>> ns_type = grid ; search neighboring grid cells
>>>>>
>>>>> ; Periodic boundary conditions: xyz, no, xy
>>>>> pbc                 = xyz ; 3-D PBC
>>>>>
>>>>> ; nblist cut-off
>>>>> ; NBOND CUTNB  (see notes on ELEC below)
>>>>> ;rlist               = 1.4 ; Cut-off for making neighbor list (short
>>>>> range
>>>>> forces). This is ignored in GPU
>>>>>
>>>>> ; OPTIONS FOR ELECTROSTATICS AND VDW
>>>>> ; Method for doing electrostatics
>>>>> ; From the CHARMM docs (ewald.doc):
>>>>> ; NBOND EWALD PMEWald KAPPa 0.34 ORDEr 6 CTOFNB 12.0 CUTNB 14.0
>>>>> coulombtype         = PME-switch ; Treatment of long range electrostatic
>>>>> interactions
>>>>> rcoulomb             = 1.2 ; long range electrostatic cut-off
>>>>>
>>>>> ; Relative dielectric constant for the medium and the reaction field
>>>>> epsilon_r           = 1
>>>>> epsilon_rf           = 1
>>>>>
>>>>> ; Method for doing Van der Waals
>>>>> ; NBOND VATOM VSWI CTONNB 10.0 CTOFNB 12.0 CUTNB 14.0
>>>>> cutoff-scheme     = Verlet
>>>>> vdw-type                = Cut-off
>>>>> vdw-modifier = Potential-switch
>>>>>
>>>>> ; cut-off lengths
>>>>> rvdw-switch         = 1.0
>>>>> rvdw                 = 1.2
>>>>>
>>>>> ; Apply long range dispersion corrections for Energy and Pressure
>>>>> ; NBOND LRC
>>>>> DispCorr             = EnerPres ; account for cut-off vdW scheme
>>>>>
>>>>> ; Seperate tables between energy group pairs
>>>>> energygrp_table         =
>>>>>
>>>>> ; Spacing for the PME/PPPM FFT grid
>>>>> ; CHARMM: EWALD recommended spacing: 0.8 A - 1.2 A and 6th Order spline
>>>>> fourierspacing       = 0.12
>>>>>
>>>>> ; EWALD/PME/PPPM parameters
>>>>> ; (possibly increase pme_order to 6 to match the CHARMM recommendation)
>>>>> pme_order           = 4
>>>>> ewald_rtol           = 1e-05
>>>>> ewald_geometry       = 3d
>>>>> epsilon_surface     = 0
>>>>>
>>>>> ; OPTIONS FOR WEAK COUPLING ALGORITHMS
>>>>> ; Temperature coupling
>>>>> Tcoupl               = V-rescale ; modified Berendsen thermostat
>>>>> tau_t               = 0.1 0.1 ; time constant, in ps
>>>>> tc-grps             = Protein non-Protein ; two coupling groups - more
>>>>> accurate
>>>>> ref_t               = 300 300 ; reference temperature, one for each
>>>>> group,
>>>>> in K
>>>>>
>>>>> ; Pressure coupling
>>>>> Pcoupl               = Parrinello-Rahman ; Pressure coupling on in NPT
>>>>> Pcoupltype           = isotropic ; uniform scaling of box vectors
>>>>>
>>>>> ; Time constant (ps), compressibility (1/bar) and reference P (bar)
>>>>> tau_p               = 2.5 ; time constant, in ps
>>>>> compressibility     = 4.5e-5 ; isothermal compressibility of water,
>>>>> bar^-1
>>>>> ref_p               = 1.0 ; reference pressure, in bar
>>>>>
>>>>> ; OPTIONS FOR BONDS
>>>>> ; CHARMM uses SHAKE with tol 1e-6
>>>>> constraints   = h-bonds   ; Constrain hydrogen bonds
>>>>> constraint_algorithm = LINCS     ; Type of constraint algorithm
>>>>> continuation   = yes       ; Do not constrain the start configuration
>>>>> (yes/no)
>>>>> lincs_iter = 1 ; accuracy of LINCS
>>>>> shake_tol             = 0.0001   ; Relative tolerance of shake
>>>>> lincs_order         = 4         ; Highest order in the expansion of the
>>>>> constraint coupling matrix
>>>>> lincs_warnangle       = 30       ; Rotate over more degrees than
>>>>>
>>>>> ; Velocity generation
>>>>> ;
>>>>> gen_vel = no ; Velocity generation is off
>>>>>
>>>>> ; the output
>>>>> ;
>>>>> nstxout       = 2500             ; Frequency to write coordinates to
>>>>> output
>>>>> trajectory file, save coordinates every 5 ps
>>>>> nstvout       = 2500             ; Frequency to write velocities to
>>>>> output
>>>>> trajectory file, save velocities every 5 ps
>>>>>
>>>>> ; Output frequency for energies to log file and energy file
>>>>> nstlog       = 2500             ; Frequency to write energies to log
>>>>> file,
>>>>> update log file every 5 ps
>>>>> nstenergy   = 2500             ; Frequency to write energies to energy
>>>>> file, save energies every 5 ps
>>>>>
>>>>> ; Output frequency and precision for xtc file
>>>>> nstxout-compressed     = 2500               ; Frequency to write
>>>>> coordinates to xtc trajectory, xtc compressed trajectory every 5 ps
>>>>> compressed-x-grps     = System           ; Group(s) to write to xtc
>>>>> trajectory
>>>>> energygrps   = System ; Group(s) to write to energy file
>>>>>
>>>>> comm_mode               = Linear              ; remove center of mass
>>>>> translation
>>>>> nstcomm                 = 1000                ; [steps] frequency of
>>>>> mass
>>>>> motion removal
>>>>> comm_grps               = System    ; group(s) for center of mass motion
>>>>> removal
>>>>>
>>>>>
>>>>> http://www.gromacs.org/Documentation/Terminology/Force_Fields/CHARMM
>>>>
>>>> Settings there apply, as well.  Note that there is no such thing as a
>>>> CHARMM27 protein force field.  What you are using is CHARMM22/CMAP.  Due
>>>> to
>>>> unfortunate file naming, the misnomer gets perpetuated.
>>>>
>>>> -Justin
>>>>
>>>> --
>>>> ==================================================
>>>>
>>>> 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
>>>>
>>>> ==================================================
>>>> --
>>>> Gromacs Users mailing list
>>>>
>>>> * Please search the archive at
>>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>>>> posting!
>>>>
>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>>
>>>> * For (un)subscribe requests visit
>>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>>> send a mail to gmx-users-request at gromacs.org.
>>>>
>>>>
>> --
>> ==================================================
>>
>> 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
>>
>> ==================================================
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-request at gromacs.org.
>>

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

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

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


More information about the gromacs.org_gmx-users mailing list