[gmx-users] All-atom vs United-atom force fields and hydrogen bonds

[Jason Hill]

Thank you Justin for the detailed and informative response. I’ve now realized that I was using inappropriate parameters for the GROMOS force field and that if I used the recommended cutoffs rvdw =! rc and I could not use the Vartlet scheme. Switching to CHARMM36 and using the recommended cutoffs puts me back into GPU powered simulations so I’ll move forward with that. Although my output is down to 20 ns/day from 50 ns/day I suppose with all atoms and a more complex water model (tip3p vs spc) that cannot be avoided.

If anyone is interested, the log file is here: https://drive.google.com/file/d/0BwAaTxAET7c5NE9Dejh2S080WXM/view?usp=sharing <https://drive.google.com/file/d/0BwAaTxAET7c5NE9Dejh2S080WXM/view?usp=sharing>

-Jason

> On 12/30/14 9:03 AM, Jason Hill wrote:
>> Hello list,
>> 
>> I?m working with some allozymes that differ by only a couple of amino acids
>> but that I have observed to have empirically very different binding
>> affinities and thermal stabilities. I?m investigating the structural basis of
>> the functional differences so I crystallized one form and am using that
>> structure as the template for a molecular dynamic approach. Among results
>> that I am interested in are hydrogen bonding patterns. I?ve been reading the
>> primary literature for the force fields and I?m sorry to say I am still
>> unclear on what I think are a couple of basic issues, ones that I hope you
>> can help me figure out.
>> 
>> 1) When a united-atom force field is used suck as GROMOS96 54a7 are the
>> hydrogens of ALL carbons in the system subsumed into the carbon and if so are
>> their position on the carbon kept track of?
>> 
> 
> Only nonpolar hydrogens are affected.  In 54A7, aromatic H are explicit, but 
> anything that is CH, CH2, or CH3 is united.  There is no information about these 
> H stored or tracked in any way; they simply don't exist.  The CHn atom types 
> have their nonbonded parameters derived such that they include the effect of the 
> attached hydrogens.
> 
>> 2) When selecting a force field is information lost in a united-atom force
>> field like GROMOS96 54a7 that renders downstream analysis by tools like
>> g_hbond less meaningful?
> 
> No.  Hydrogen bonding (at least within g_hbond, because atom names are 
> hard-coded) involves only polar groups.
> 
>> 
>> 3) One of my favorite up to date tutorials
>> (http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/index.html <http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/index.html>
>> <http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/index.html <http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/index.html>>)
>> uses the OPLS-AA/L force field but a lot of the literature seems to indicate
>> that the newer force fields perform much better (blog review:
>> http://md.chem.rug.nl/cgmartini/index.php/blog/265-comparingforcefields <http://md.chem.rug.nl/cgmartini/index.php/blog/265-comparingforcefields>
>> <http://md.chem.rug.nl/cgmartini/index.php/blog/265-comparingforcefields <http://md.chem.rug.nl/cgmartini/index.php/blog/265-comparingforcefields>>)!
>> Is the all atom parameterization so much more important than factors like
>> agreement with NMR structure for down stream analysis for g_hbond or other
>> things that OPLS-AA/L would still be a better choice?
>> 
> 
> The use of OPLS-AA in my tutorial is for convenience for myself; I reasoned I 
> would get fewer questions like "where did all of my hydrogens go?" so I went 
> with an all-atom force field.  I leave the tutorial as-is because (1) I have a 
> data set that I know is highly reproducible and (2) I really don't feel like 
> rewriting it.  I'm fairly sure I added a disclaimer a while back that indicated 
> that users should always look into suitable force fields, rather than thinking 
> the tutorial is simply a recipe of "this is how you do a protein MD simulation." 
>  It's not.
> 
>> 4) For a protein in water simulation in which you were investigating the
>> movement and stability of the protein, what force field would you use, and
>> why?
>> 
> 
> Everyone will have their own reasons here, and this is a never-ending battle in 
> the literature, depending on what metrics you're using.  I believe strongly in 
> CHARMM36, because I see the results it produces every day in our group.  But 
> then too, I believe in polarizable force fields over additive models, anyway :)
> 
> Most force fields will perform reasonably well for a fully folded protein. 
> Gromos96 53A6 has some issues with helices, but 54A7 seems to have largely fixed 
> these.  AMBER, CHARMM, and OPLS-AA variants all seem to do well, in general.
> 
> -Justin