[gmx-users] Fwd: Help with Gromacs 5.1 mdp options for CHARMM27 force field
Rakesh Ramachandran
korenpol3 at gmail.com
Sun Dec 6 14:58:22 CET 2015
Hi Justin,
Thanks for the reply. Sorry if I am asking too many questions. I wanted
to clarify some more doubts just to make sure I am performing the
simulation the right way.
1) What is the ideal Tcoupl option for CHARMM36 Nose-hoover or V-rescale.
Moreover I get the following note and warning when I use Nose-Hoover
thermostat.
NOTE 1 [file npt.mdp]:
leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1
WARNING 1 [file npt.mdp]:
For proper integration of the Nose-Hoover thermostat, tau-t (1) should be
at least 20 times larger than nsttcouple*dt (0.08)
So is V-rescale more optimal for Verlet cutoff and GPU based Gromacs.
2) Moreover I see in some publications that 0.9 nm to 1.0 nm for minimum
protein-edge distance is being used for Gromacs. Even in your tutorial you
mention using 1.0 nm. Taking into account minimum image convention should
it not be greater than 1.2 nm or is it highly dependent on the protein you
are simulating or rcouloumb cutoff or is dependent on Debye length.
I would be grateful if you could clear my doubts on this aspect or point me
to the relevant literature. I am planning to use 1.25 nm for a salt conc of
300 mM MgCl2 with epsilon_r = 1 or should I use higher cutoff such as 1.4
nm. Moreover is the dielectric constant of the medium 82 for TIP3P water
model at 300K (ref http://www1.lsbu.ac.uk/water/water_models.html),
3) In spite of running the NPT equilibration for 5 ns I am getting these
values and the RMSD values are also too high when the pressure parameter is
measured. In your tutorial you had suggested to run longer NPT
equilibration if the density value does not converge to around 1000 kg/m^3
and pressure value around 1 bar. I found the same issue when I ran with
CHARMM27 force field and V-rescale thermostat. But on plotting I find the
results are similar to your plots in the tutorial but for 5 ns in both the
cases (CHARMM 27 and 36).
for 1.25 nm protein-edge distance
Pressure 1.52207 0.52 146.108 1.62992 (bar)
Density 1056.64 0.029 2.7515 -0.0336808
(kg/m^3)
for 1.4 nm protein-edge distance
Pressure 0.321022 0.94 137.028 2.4858 (bar)
Density 1054.45 0.053 2.67943 0.264365
(kg/m^3)
Regards
Rakesh
On 4 December 2015 at 18:04, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>
> 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
>
> ==================================================
> --
> 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.
>
More information about the gromacs.org_gmx-users
mailing list