[gmx-users] Bond lengths

[Justin Lemkul]

On 1/4/17 2:28 AM, suniba shuaib wrote:
> Dear Users
> I have performed a 200 ns MD simulation of a protein and a small molecule
> inhibitor of 30 atoms, using GROMOS 43a1 ff and Gromacs v 5.0. Taking
> reference from a number of papers and similar systems, I have used LINCS
> algorithm to constrain all bond lengths.
> I have reported interactions like pi-pi and hydrogen bonding between
> protein and ligand plus some other analysis and  recently i submitted my
> article. Reviewers have asked few question one of which is:
>
> *Q: Why you have constrained all bond lengths? A good force field should
> permit free variation of the bond lengths and only high frequency motions
> (C-H) should be constrained.*
>
> As i visited mailing list, I found an answer by Justin Lemkul:
>
> *"The effects of H-bond vibrations aren't of interest for the problems
> at hand; if you want to model bond formation, you can't accept a model
> that has no "electron vibrations," and *that* drives your choice of
> time step."*
>
> Another reply by Mark Abraham:
> "*Force field parameters are (unfortunately) co-optimized with and
> dependent on*
>
> *implementation issues like size of time step and lengths of cut-offs. The
> use of constraints=all-bonds couples constraints across a whole polymer; if
> that's a protein simulated using a spatial domain decomposition, those
> atoms are widely distributed on the machine and the communication required
> for the iterative constraint solver can become rate limiting at extreme
> scale. Reducing the constraint iteration to a single node is unattractive
> if the rest of the machine lies idle - which is part of why we'd like to
> improve the ability to deploy task parallelism in GROMACS at scale. The
> main alternative is using shorter time steps (maybe no constraints at all!)
> - if they only way you can get 10x greater sampling rate for a single
> simulation is by using 100x the hardware at one tenth of the time step,
> then that could be a serious option.
>
> I think that feeling that "constrained H bonds" is natural and "constrained
> heavy-atom bonds" is artificial is wrongly conflating experience with
> correctness."*
>
>
> After reading these answers, I am little confused how should i respond
> to the raised question. In my case, I believe, it was not necessary to
> keep the system unconstrained. Plus, from Mark's answer, I concluded
> that it is better to keep all bonds constrained? Please clear my
> confusion.
>

In both the 43A1 and the 53A5/53A6 papers, the authors state that all bonds were 
constrained to their equilibrium lengths.  Thus, rigid bonds are assumed to be 
part of the inherent parametrization of the GROMOS force field and thus your 
usage is correct.  This is not true for OPLS, AMBER, and CHARMM, which only fix 
bonds to H atoms.

-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

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