force field research paper

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

• View all journals
• Explore content
• About the journal
• Publish with us
• Sign up for alerts
• Published: 29 March 2021

Martini 3: a general purpose force field for coarse-grained molecular dynamics

• Paulo C. T. Souza   ORCID: orcid.org/0000-0003-0660-1301 1 , 2 ,
• Riccardo Alessandri   ORCID: orcid.org/0000-0003-1948-5311 1 ,
• Jonathan Barnoud 1 , 3 ,
• Sebastian Thallmair   ORCID: orcid.org/0000-0002-3396-5840 1 , 4 ,
• Ignacio Faustino 1 ,
• Fabian Grünewald   ORCID: orcid.org/0000-0001-6979-1363 1 ,
• Ilias Patmanidis 1 ,
• Haleh Abdizadeh 1 ,
• Bart M. H. Bruininks   ORCID: orcid.org/0000-0001-5136-0864 1 ,
• Tsjerk A. Wassenaar 1 ,
• Peter C. Kroon   ORCID: orcid.org/0000-0001-9273-4850 1 ,
• Josef Melcr   ORCID: orcid.org/0000-0003-4729-3990 1 ,
• Vincent Nieto 2 ,
• Valentina Corradi   ORCID: orcid.org/0000-0001-6546-9780 5 ,
• Hanif M. Khan 5 , 6 ,
• Jan Domański 7 , 8 ,
• Matti Javanainen   ORCID: orcid.org/0000-0003-4858-364X 9 , 10 ,
• Hector Martinez-Seara   ORCID: orcid.org/0000-0001-9716-1713 9 ,
• Nathalie Reuter   ORCID: orcid.org/0000-0002-3649-7675 6 ,
• Robert B. Best 8 ,
• Ilpo Vattulainen   ORCID: orcid.org/0000-0001-7408-3214 10 , 11 ,
• Luca Monticelli   ORCID: orcid.org/0000-0002-6352-4595 2 ,
• Xavier Periole 12 ,
• D. Peter Tieleman   ORCID: orcid.org/0000-0001-5507-0688 5 ,
• Alex H. de Vries 1 &
• Siewert J. Marrink 1  

Nature Methods volume  18 ,  pages 382–388 ( 2021 ) Cite this article

38k Accesses

578 Citations

148 Altmetric

Metrics details

• Computational biology and bioinformatics
• Computational biophysics
• Computational models

The coarse-grained Martini force field is widely used in biomolecular simulations. Here we present the refined model, Martini 3 ( http://cgmartini.nl ), with an improved interaction balance, new bead types and expanded ability to include specific interactions representing, for example, hydrogen bonding and electronic polarizability. The updated model allows more accurate predictions of molecular packing and interactions in general, which is exemplified with a vast and diverse set of applications, ranging from oil/water partitioning and miscibility data to complex molecular systems, involving protein–protein and protein–lipid interactions and material science applications as ionic liquids and aedamers.

This is a preview of subscription content, access via your institution

Access options

Access Nature and 54 other Nature Portfolio journals

Get Nature+, our best-value online-access subscription

24,99 € / 30 days

cancel any time

Subscribe to this journal

Receive 12 print issues and online access

251,40 € per year

only 20,95 € per issue

Buy this article

• Purchase on SpringerLink
• Instant access to full article PDF

Prices may be subject to local taxes which are calculated during checkout

Similar content being viewed by others

The HADDOCK2.4 web server for integrative modeling of biomolecular complexes

Physics-driven coarse-grained model for biomolecular phase separation with near-quantitative accuracy

Macromolecular modeling and design in Rosetta: recent methods and frameworks

Data availability.

Force-field parameters and procedures (for example, tutorials) are publicly available at http://cgmartini.nl . Simulation data (for example, trajectories) supporting the results of this paper are available from the corresponding authors upon reasonable request.

Code availability

Martinize2 (for which the manuscript is in preparation) and Martinate codes used in this work are publicly available at https://github.com/marrink-lab/ . For more detailed information, see Supplementary Codes .

Bottaro, S. & Lindorff-Larsen, K. Biophysical experiments and biomolecular simulations: a perfect match? Science 361 , 355–360 (2018).

Article   CAS   PubMed   Google Scholar  

Ingólfsson, H. I. et al. The power of coarse graining in biomolecular simulations. Wiley Interdiscip. Rev. Comput. Mol. Sci. 4 , 225–248 (2014).

Article   PubMed   Google Scholar  

Marrink, S. J., De Vries, A. H. & Mark, A. E. Coarse grained model for semiquantitative lipid simulations. J. Phys. Chem. B 108 , 750–760 (2004).

Article   CAS   Google Scholar  

Marrink, S. J., Risselada, H. J., Yefimov, S., Tieleman, D. P. & de Vries, A. H. The MARTINI force field: coarse grained model for biomolecular simulations. J. Phys. Chem. B 111 , 7812–7824 (2007).

Uusitalo, J. J., Ingólfsson, H. I., Akhshi, P., Tieleman, D. P. & Marrink, S. J. Martini coarse-grained force field: extension to DNA. J. Chem. Theory Comput. 11 , 3932–3945 (2015).

Abellón-Ruiz, J. et al. Structural basis for maintenance of bacterial outer membrane lipid asymmetry. Nat. Microbiol. 2 , 1616–1623 (2017).

Yen, H. Y. et al. PtdIns(4,5)P2 stabilizes active states of GPCRs and enhances selectivity of G-protein coupling. Nature 559 , 423–427 (2018).

Article   CAS   PubMed   PubMed Central   Google Scholar  

Van Eerden, F. J., Melo, M. N., Frederix, P. W. J. M., Periole, X. & Marrink, S. J. Exchange pathways of plastoquinone and plastoquinol in the photosystem II complex. Nat. Commun. 8 , 15214 (2017).

Article   PubMed   PubMed Central   Google Scholar  

Vögele, M., Köfinger, J. & Hummer, G. Hydrodynamics of diffusion in lipid membrane simulations. Phys. Rev. Lett. 120 , 268104 (2018).

Agostino, M. D., Risselada, H. J., Lürick, A., Ungermann, C. & Mayer, A. A tethering complex drives the terminal stage of SNARE-dependent membrane fusion. Nature 551 , 634–638 (2017).

Jeena, M. T. et al. Mitochondria localization induced self-assembly of peptide amphiphiles for cellular dysfunction. Nat. Commun. 8 , 26 (2017).

Jiang, Z. et al. Subnanometre ligand-shell asymmetry leads to Janus-like nanoparticle membranes. Nat. Mater. 14 , 912–917 (2015).

Maingi, V. et al. Stability and dynamics of membrane-spanning DNA nanopores. Nat. Commun. 8 , 14784 (2017).

Frederix, P. W. J. M. et al. Exploring the sequence space for (tri-)peptide self-assembly to design and discover new hydrogels. Nat. Chem. 7 , 30–37 (2015).

Bochicchio, D., Salvalaglio, M. & Pavan, G. M. Into the dynamics of a supramolecular polymer at submolecular resolution. Nat. Commun. 8 , 147 (2017).

Stark, A. C., Andrews, C. T. & Elcock, A. H. Toward optimized potential functions for protein-protein interactions in aqueous solutions: osmotic second virial coefficient calculations using the MARTINI coarse-grained force field. J. Chem. Theory Comput. 9 , 4176–4185 (2013).

Javanainen, M., Martinez-Seara, H. & Vattulainen, I. Excessive aggregation of membrane proteins in the Martini model. PLoS ONE 12 , e0187936 (2017).

Schmalhorst, P. S., Deluweit, F., Scherrers, R., Heisenberg, C.-P. & Sikora, M. Overcoming the limitations of the MARTINI force field in simulations of polysaccharides. J. Chem. Theory Comput. 13 , 5039–5053 (2017).

Alessandri, R. et al. Pitfalls of the Martini model. J. Chem. Theory Comput. 15 , 5448–5460 (2019).

Uusitalo, J. J., Ingólfsson, H. I., Marrink, S. J. & Faustino, I. Martini coarse-grained force field: extension to RNA. Biophys. J. 113 , 246–256 (2017).

Ben-Naim, A. Molecular Theory of Solutions (Oxford Univ. Press, 2006).

Ploetz, E. A., Bentenitis, N. & Smith, P. E. Kirkwood–Buff integrals for ideal solutions. J. Chem. Phys. 132 , 164501 (2010).

Zych, A. J. & Iverson, B. L. Synthesis and conformational characterization of tethered, self-complexing 1,5-dialkoxynaphthalene/1,4,5,8-naphthalenetetracarboxylic diimide systems. J. Am. Chem. Soc. 122 , 8898–8909 (2000).

Gabriel, G. J. & Iverson, B. L. Aromatic oligomers that form hetero duplexes in aqueous solution. J. Am. Chem. Soc. 124 , 15174–15175 (2002).

Liu, W. et al. Structural basis for allosteric regulation of GPCRs by sodium ions. Science 337 , 232–236 (2012).

Gao, Z. G. & Ijzerman, A. P. Allosteric modulation of A(2A) adenosine receptors by amiloride analogues and sodium ions. Biochem. Pharmacol. 60 , 669–676 (2000).

Okur, H. I. et al. Beyond the Hofmeister series: ion-specific effects on proteins and their biological functions. J. Phys. Chem. B 121 , 1997–2014 (2017).

Dupont, D., Depuydt, D. & Binnemans, K. Overview of the effect of salts on biphasic ionic liquid/water solvent extraction systems: anion exchange, mutual solubility, and thermomorphic properties. J. Phys. Chem. B 119 , 6747–6757 (2015).

Naert, P., Rabaey, K. & Stevens, C. V. Ionic liquid ion exchange: exclusion from strong interactions condemns cations to the most weakly interacting anions and dictates reaction equilibrium. Green Chem. 20 , 4277–4286 (2018).

Khan, H. M. et al. Capturing choline-aromatics cation–π interactions in the MARTINI force field. J. Chem. Theory Comput. 16 , 2550–2560 (2020).

Tanaka, K., Caaveiro, J. M. M., Morante, K., González-Manãs, J. M. & Tsumoto, K. Structural basis for self-assembly of a cytolytic pore lined by protein and lipid. Nat. Commun. 6 , 6337 (2015).

Huang, G., Willems, K., Soskine, M., Wloka, C. & Maglia, G. Electro-osmotic capture and ionic discrimination of peptide and protein biomarkers with FraC nanopores. Nat. Commun. 8 , 935 (2017).

Alessandri, R., Uusitalo, J. J., De Vries, A. H., Havenith, R. W. A. & Marrink, S. J. Bulk heterojunction morphologies with atomistic resolution from coarse-grain solvent evaporation simulations. J. Am. Chem. Soc. 139 , 3697–3705 (2017).

Chiu, M. Y., Jeng, U. S., Su, C. H., Liang, K. S. & Wei, K. H. Simultaneous use of small- and wide-angle X-ray techniques to analyze nanometerscale phase separation in polymer heterojunction solar cells. Adv. Mater. 20 , 2573–2578 (2008).

Petrov, D. & Zagrovic, B. Are current atomistic force fields accurate enough to study proteins in crowded environments? PLoS Comput. Biol. 10 , e1003638 (2014).

Højgaard, C. et al. A soluble, folded protein without charged amino acid residues. Biochemistry 55 , 3949–3956 (2016).

Ruckenstein, E. & Shulgin, I. L. Effect of salts and organic additives on the solubility of proteins in aqueous solutions. Adv. Colloid Interface Sci. 123–126 , 97–103 (2006).

Zhou, F. X., Cocco, M. J., Russ, W. P., Brunger, A. T. & Engelman, D. M. Interhelical hydrogen bonding drives strong interactions in membrane proteins. Nat. Struct. Biol. 7 , 154–160 (2000).

Zhou, F. X., Merianos, H. J., Brunger, A. T. & Engelman, D. M. Polar residues drive association of polyleucine transmembrane helices. Proc. Natl Acad. Sci. USA 98 , 2250–2255 (2001).

Grau, B. et al. The role of hydrophobic matching on transmembrane helix packing in cells. Cell Stress 1 , 90–106 (2017).

Chen, L., Merzlyakov, M., Cohen, T., Shai, Y. & Hristova, K. Energetics of ErbB1 transmembrane domain dimerization in lipid bilayers. Biophys. J. 96 , 4622–4630 (2009).

Artemenko, E. O., Egorova, N. S., Arseniev, A. S. & Feofanov, A. V. Transmembrane domain of EphA1 receptor forms dimers in membrane-like environment. Biochim. Biophys. Acta 1778 , 2361–2367 (2008).

Sarabipour, S. & Hristova, K. Glycophorin A transmembrane domain dimerization in plasma membrane vesicles derived from CHO, HEK 293T, and A431 cells. Biochim. Biophys. Acta - Biomembr. 1828 , 1829–1833 (2013).

Chen, L., Novicky, L., Merzlyakov, M., Hristov, T. & Hristova, K. Measuring the energetics of membrane protein dimerization in mammalian membranes. J. Am. Chem. Soc. 132 , 3628–3635 (2010).

Nash, A., Notman, R. & Dixon, A. M. De novo design of transmembrane helix–helix interactions and measurement of stability in a biological membrane. Biochim. Biophys. Acta - Biomembr. 1848 , 1248–1257 (2015).

Finger, C. et al. The stability of transmembrane helix interactions measured in a biological membrane. J. Mol. Biol. 358 , 1221–1228 (2006).

Hong, H., Blois, T. M., Cao, Z. & Bowie, J. U. Method to measure strong protein–protein interactions in lipid bilayers using a steric trap. Proc. Natl Acad. Sci. USA 107 , 19802–19807 (2010).

Sparr, E. et al. Self-association of transmembrane α-helices in model membranes: importance of helix orientation and role of hydrophobic mismatch. J. Biol. Chem. 280 , 39324–39331 (2005).

MacKenzie, K. R., Prestegard, J. H. & Engelman, D. M. Transmembrane helix dimer: structure and implications. Science 276 , 131–133 (1997).

Trenker, R., Call, M. E. & Call, M. J. Crystal structure of the glycophorin A transmembrane dimer in lipidic cubic phase. J. Am. Chem. Soc. 137 , 15676–15679 (2015).

Domański, J., Sansom, M. S. P., Stansfeld, P. J. & Best, R. B. Balancing force field protein–lipid interactions to capture transmembrane helix–helix association. J. Chem. Theory Comput. 14 , 1706–1715 (2018).

Souza, P. C. T., Thallmair, S., Marrink, S. J. & Mera-Adasme, R. An allosteric pathway in copper, zinc superoxide dismutase unravels the molecular mechanism of the G93A amyotrophic lateral sclerosis-linked mutation. J. Phys. Chem. Lett. 10 , 7740–7744 (2019).

Brini, E. et al. Systematic coarse-graining methods for soft matter simulations-a review. Soft Matter 9 , 2108–2119 (2013).

Foley, T. T., Shell, M. S. & Noid, W. G. The impact of resolution upon entropy and information in coarse-grained models. J. Chem. Phys. 143 , 243104 (2015).

Wagner, J. W., Dama, J. F., Durumeric, A. E. P. & Voth, G. A. On the representability problem and the physical meaning of coarse-grained models. J. Chem. Phys. 145 , 044108 (2016).

Wörner, S. J., Bereau, T., Kremer, K. & Rudzinski, J. F. Direct route to reproducing pair distribution functions with coarse-grained models via transformed atomistic cross correlations. J. Chem. Phys. 151 , 244110 (2019).

Noid, W. G., Chu, J. W., Ayton, G. S. & Voth, G. A. Multiscale coarse-graining and structural correlations: connections to liquid-state theory. J. Phys. Chem. B. 111 , 4116–4127 (2007).

Wu, Z., Cui, Q. & Yethiraj, A. Driving force for the association of hydrophobic peptides: the importance of electrostatic interactions in coarse-grained water models. J. Phys. Chem. Lett. 2 , 1794–1798 (2011).

Jin, J., Yu, A. & Voth, G. A. Temperature and phase transferable bottom-up coarse-grained models. J. Chem. Theory Comput. 16 , 6823–6842 (2020).

Yesylevskyy, S. O., Schäfer, L. V., Sengupta, D. & Marrink, S. J. Polarizable water model for the coarse-grained MARTINI force field. PLoS Comput. Biol. 6 , e1000810 (2010).

Michalowsky, J., Schäfer, L. V., Holm, C. & Smiatek, J. A refined polarizable water model for the coarse-grained MARTINI force field with long-range electrostatic interactions. J. Chem. Phys. 146 , 054501 (2017).

Marrink, S. J. & Tieleman, D. P. Perspective on the Martini model. Chem. Soc. Rev. 42 , 6801–22 (2013).

Bruininks, B. M. H., Souza, P. C. T. & Marrink, S. J. in Biomolecular Simulations: Methods and Protocols (eds Bonomi, M. & Camilloni, C.) 105–127 (Humana Press Inc., 2019).

Liu, J. et al. Enhancing molecular n-type doping of donor-acceptor copolymers by tailoring side chains. Adv. Mater. 30 , 1704630 (2018).

Article   Google Scholar  

Vazquez-Salazar, L. I., Selle, M., de Vries, A., Marrink, S. J. & Souza, P. C. T. Martini coarse-grained models of imidazolium-based ionic liquids: from nanostructural organization to liquid–liquid extraction. Green Chem. 22 , 7376–7386 (2020).

Souza, P. C. T. et al. Protein–ligand binding with the coarse-grained Martini model. Nat. Commun. 11 , 3714 (2020).

López, C. A. et al. Martini coarse-grained force field: extension to carbohydrates. J. Chem. Theory Comput. 5 , 3195–3210 (2009).

Monticelli, L. et al. The MARTINI coarse-grained force field: extension to proteins. J. Chem. Theory Comput. 4 , 819–834 (2008).

Grunewald, F., Rossi, G., de Vries, A. H., Marrink, S. J. & Monticelli, L. Transferable MARTINI model of poly(ethylene oxide). J. Phys. Chem. B 122 , 7436–7449 (2018).

de Jong, D. H. et al. Improved parameters for the martini coarse-grained protein force field. J. Chem. Theory Comput. 9 , 687–97 (2013).

Herzog, F. A., Braun, L., Schoen, I. & Vogel, V. Improved side chain dynamics in MARTINI simulations of protein–lipid interfaces. J. Chem. Theory Comput. 12 , 2446–2458 (2016).

Poma, A. B., Cieplak, M. & Theodorakis, P. E. Combining the MARTINI and structure-based coarse-grained approaches for the molecular dynamics studies of conformational transitions in proteins. J. Chem. Theory Comput. 13 , 1366–1374 (2017).

Periole, X., Cavalli, M., Marrink, S.-J. & Ceruso, M. A. Combining an elastic network with a coarse-grained molecular force field: structure, dynamics, and intermolecular recognition. J. Chem. Theory Comput. 5 , 2531–2543 (2009).

Wassenaar, T. A., Ingólfsson, H. I., Böckmann, R. A., Tieleman, D. P. & Marrink, S. J. Computational lipidomics with insane: a versatile tool for generating Custom membranes for molecular simulations. J. Chem. Theory Comput. 11 , 2144–2155 (2015).

Melo, M. N., Ingólfsson, H. I. & Marrink, S. J. Parameters for Martini sterols and hopanoids based on a virtual-site description. J. Chem. Phys. 143 , 243152 (2015).

López, C. A., Sovova, Z., van Eerden, F. J., de Vries, A. H. & Marrink, S. J. Martini force field parameters for glycolipids. J. Chem. Theory Comput. 9 , 1694–1708 (2013).

Carpenter, T. S. et al. Capturing phase behavior of ternary lipid mixtures with a refined Martini coarse-grained force field. J. Chem. Theory Comput. 14 , 6050–6062 (2018).

de Jong, D. H., Baoukina, S., Ingólfsson, H. I. & Marrink, S. J. Martini straight: boosting performance using a shorter cutoff and GPUs. Comput. Phys. Commun. 199 , 1–7 (2016).

Hockney, R. W., Goel, S. P. & Eastwood, J. W. Quiet high-resolution computer models of a plasma. J. Comput. Phys. 14 , 148–158 (1974).

Páll, S. & Hess, B. A flexible algorithm for calculating pair interactions on SIMD architectures. Comput. Phys. Commun. 184 , 2641–2650 (2013).

Verlet, L. Computer ‘experiments’ on classical fluids. I. Thermodynamical properties of Lennard–Jones molecules. Phys. Rev. 159 , 98–103 (1967).

Tironi, I. G., Sperb, R., Smith, P. E. & Van Gunsteren, W. F. A generalized reaction field method for molecular dynamics simulations. J. Chem. Phys. 102 , 5451–5459 (1995).

Essmann, U. et al. A smooth particle mesh Ewald method. J. Chem. Phys. 103 , 8577–8593 (1995).

Bussi, G., Donadio, D. & Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 126 , 014101 (2007).

Parrinello, M. & Rahman, A. Polymorphic transitions in single crystals: a new molecular dynamics method. J. Appl. Phys. 52 , 7182–7190 (1981).

Abraham, M. J. et al. GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 1–2 , 19–25 (2015).

Van Der Spoel, D. et al. GROMACS: fast, flexible, and free. J. Comput. Chem. 26 , 1701–1718 (2005).

Wassenaar, T. A., Ingólfsson, H. I., Prieß, M., Marrink, S. J. & Schäfer, L. V. Mixing MARTINI: electrostatic coupling in hybrid atomistic-coarse-grained biomolecular simulations. J. Phys. Chem. B. 117 , 3516–3530 (2013).

Wassenaar, T. A. et al. High-throughput simulations of dimer and trimer assembly of membrane proteins. The DAFT approach. J. Chem. Theory Comput. 11 , 2278–91 (2015).

Humphrey, W., Dalke, A. & Schulten, K. VMD—visual molecular dynamics. J. Molec. Graph. 14 , 33–38 (1996).

Gowers, R. J. et al. MDAnalysis: a Python package for the rapid analysis of molecular dynamics simulations. in Proc. 15th Python Sci. Conference 98–105 (2016).

Download references

Acknowledgements

We thank all members of the S.J.M. group and also external users for testing Martini 3 in its open-beta version. In particular, we thank C. F. E. Schroer, P. W. J. M. Frederix, W. Pezeshkian, M. N. Melo, H. I. Ingólfsson, M. Tsanai, M. König, P. A. Vainikka, T. Zijp, L. Gaifas, J. H. van der Woude, M. Espinoza Cangahuala, M. Scharte, J. Cruiming, L. M. van der Sleen, V. Verduijn, A. H. Beck Frederiksen, B. Schiøtt, M. Sikora, P. Schmalhorst, R. A. Moreira, A. B. Poma, K. Pluhackova, C. Arnarez, C. A. López, E. Jefferys and M. S. P. Sansom for their preliminary tests with a lot of different systems including aedamers, sugars, amino acids, deep eutectic solvents, lipids, peptides and proteins. We also thank the Center for Information Technology of the University of Groningen for providing access to the Peregrine high-performance computing cluster. We acknowledge the National Computing Facilities Foundation of The Netherlands Organization for Scientific Research (NWO), CSC–IT Center for Science Ltd (Espoo, Finland) and CINES (France) for providing computing time. Work in the S.J.M. group was supported by an European Research Council advanced grant no. ‘COMP-MICR-CROW-MEM’. R.A. thanks the NWO (Graduate Programme Advanced Materials, no. 022.005.006) for financial support. L.M. acknowledges the Institut National de la Santé et de la Recherche Medicale and the Agence Nationale de la Recherche for funding (grant no. ANR-17-CE11-0003) and GENCI-CINES for computing time (grant no. A0060710138). S.T. acknowledges the support from the European Commission via a Marie Skłodowska-Curie Actions individual fellowship (MicroMod-PSII, grant agreement no. 748895). M.J. thanks the Emil Aaltonen foundation for financial support. I.V. thanks the Academy of Finland (Center of Excellence program (grant no. 307415)), Sigrid Juselius Foundation, the Helsinki Institute of Life Science fellow program and the HFSP (research grant no. RGP0059/2019). R.B.B. and J.D. were supported by the intramural research program of the NIDDK, NIH. Their work used the computational resources of the NIH HPC Biowulf cluster ( http://hpc.nih.gov ). H.M.-S. acknowledges the Czech Science Foundation (grant no. 19-19561S). J.B. acknowledges funding from the TOP grant from S.J.M. (NWO) and the EPSRC program grant no. EP/P021123/1. Work in D.P.T.’s group is supported by the Natural Sciences and Engineering Research Council (Canada) and Compute Canada, funded by the Canada Foundation for Innovation. D.P.T. acknowledges further support from the Canada Research Chairs program. N.R. acknowledges funding from the Norwegian Research Council (FRIMEDBIO nos. 251247 and 288008) and computational resources from UNINETT SIGMA2 AS (grant no. NN4700K). H.M.K. acknowledges funding from the University of Calgary through the ‘Eyes High Postdoctoral Fellowship’ program.

Author information

Authors and affiliations.

Groningen Biomolecular Sciences and Biotechnology Institute and Zernike Institute for Advanced Material, University of Groningen, Groningen, the Netherlands

Paulo C. T. Souza, Riccardo Alessandri, Jonathan Barnoud, Sebastian Thallmair, Ignacio Faustino, Fabian Grünewald, Ilias Patmanidis, Haleh Abdizadeh, Bart M. H. Bruininks, Tsjerk A. Wassenaar, Peter C. Kroon, Josef Melcr, Alex H. de Vries & Siewert J. Marrink

Molecular Microbiology and Structural Biochemistry, UMR 5086 CNRS and University of Lyon, Lyon, France

Paulo C. T. Souza, Vincent Nieto & Luca Monticelli

Intangible Realities Laboratory, University of Bristol, School of Chemistry, Bristol, UK

Jonathan Barnoud

Frankfurt Institute for Advanced Studies, Frankfurt am Main, Germany

Sebastian Thallmair

Centre for Molecular Simulation and Department of Biological Sciences, University of Calgary, Calgary, Alberta, Canada

Valentina Corradi, Hanif M. Khan & D. Peter Tieleman

Department of Chemistry and Computational Biology Unit, University of Bergen, Bergen, Norway

Hanif M. Khan & Nathalie Reuter

Department of Biochemistry, University of Oxford, Oxford, UK

Jan Domański

Laboratory of Chemical Physics, National Institute of Diabetes and Digestive and Kidney Diseases, National Institutes of Health, Bethesda, MD, USA

Jan Domański & Robert B. Best

Institute of Organic Chemistry and Biochemistry, Czech Academy of Sciences, Prague, Czech Republic

Matti Javanainen & Hector Martinez-Seara

Computational Physics Laboratory, Tampere University, Tampere, Finland

Matti Javanainen & Ilpo Vattulainen

Department of Physics, University of Helsinki, Helsinki, Finland

Ilpo Vattulainen

Department of Chemistry, Aarhus University, Aarhus C, Denmark

Xavier Periole

You can also search for this author in PubMed   Google Scholar

Contributions

P.C.T.S. and S.J.M. conceived the project with suggestions from R.A., A.H.V., J.B. and S.T. P.C.T.S. generated and optimized all bead parameters. P.C.T.S., R.A. and J.B. generated the topology and bonded parameters of all CG models with suggestions from S.T. and I.F. P.C.T.S., R.A., A.H.V. and F.G. performed the simulations and analysis involving transfer free energies, solvent and polymer properties. P.C.T.S., S.T., J.B. and J.M. performed the simulations and analysis involving lipid bilayers. P.C.T.S., I.F. and R.A. performed the simulations and analysis involving nucleobases. P.C.T.S., I.P. and A.H.V. generated the models and performed the simulations and analysis involving aedamers. P.C.T.S. and F.G. generated the models and performed the simulations and analysis involving ionic liquids and ionic water solutions. R.A. generated the models and performed the simulations and analysis involving bulk heterojunctions, with suggestions from L.M. regarding the fullerene model. P.C.T.S., J.B., H.A., R.A., B.M.H.B., S.T., J.M., V.N., X.P., M.J., H.M.K., J.D., V.C. and H.M.-S. performed the simulations and analysis involving amino acids, peptides and proteins. J.B., T.A.W., P.C.K. and S.T. developed some tools and scripts used to generate the CG models and to run the molecular dynamics simulations. L.M., R.B.B., P.T., N.R., I.V., A.H.V. and S.J.M. provided guidance and supervision in the studies performed by their respective group members and collaborators. P.C.T.S. and S.J.M. wrote the main manuscript, with contributions from all the authors. P.C.T.S. prepared the figures with contributions from R.A., B.M.H.B., H.M.K. and A.H.V. P.C.T.S. wrote the Methods with contributions from all the authors. P.C.T.S. wrote the Supplementary Information , with contributions from all the authors. All the authors revised and approved the final version of the manuscript, Methods and Supplementary Information .

Corresponding authors

Correspondence to Paulo C. T. Souza or Siewert J. Marrink .

Ethics declarations

Competing interests.

The authors declare no competing interests.

Additional information

Peer review information Nature Methods thanks the anonymous reviewers for their contribution for the peer review of this work. Arunima Singh was the primary editor on this article and managed its editorial process and peer review in collaboration with the rest of the editorial team.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Supplementary information.

Supplementary Notes, Results and Codes.

Reporting Summary

Supplementary table 1.

Comparison between simulation and experimental results

Rights and permissions

Reprints and permissions

About this article

Cite this article.

Souza, P.C.T., Alessandri, R., Barnoud, J. et al. Martini 3: a general purpose force field for coarse-grained molecular dynamics. Nat Methods 18 , 382–388 (2021). https://doi.org/10.1038/s41592-021-01098-3

Download citation

Received : 16 June 2020

Accepted : 22 February 2021

Published : 29 March 2021

Issue Date : April 2021

DOI : https://doi.org/10.1038/s41592-021-01098-3

Share this article

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

This article is cited by

The molecular picture of the local environment in a stable model coacervate.

• Atanu Baksi
• Hasan Zerze
• Gül H. Zerze

Communications Chemistry (2024)

Inhibition and transport mechanisms of the ABC transporter hMRP5

• Chenyang Xue
• Zhongmin Liu

Nature Communications (2024)

Vesicle protrusion induced by antimicrobial peptides suggests common carpet mechanism for short antimicrobial peptides

• Danilo K. Matsubara
• Iolanda M. Cuccovia

Scientific Reports (2024)

Disordered regions in the IRE1α ER lumenal domain mediate its stress-induced clustering

• Paulina Kettel
• Laura Marosits
• G Elif Karagöz

The EMBO Journal (2024)

Structure of the MlaC-MlaD complex reveals molecular basis of periplasmic phospholipid transport

• Peter Wotherspoon
• Hannah Johnston
• Timothy J. Knowles

Quick links

• Explore articles by subject
• Guide to authors
• Editorial policies

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Academia.edu no longer supports Internet Explorer.

To browse Academia.edu and the wider internet faster and more securely, please take a few seconds to  upgrade your browser .

•  We're Hiring!
•  Help Center

Force Field Analysis

• Most Cited Papers
• Most Downloaded Papers
• Newest Papers
• Last »
• Neoliberalism and community organizing Follow Following
• Radical Community Organizing Follow Following
• Teaching Macro Practice Follow Following
• Community Organizing Follow Following
• Catia Follow Following
• Decision and Risk Analysis Follow Following
• COLLISION Follow Following
• Vibrational Analysis Follow Following
• 8051 Microcontroller Follow Following
• Ansys Follow Following

Enter the email address you signed up with and we'll email you a reset link.

• Academia.edu Journals
•   We're Hiring!
•   Help Center
• Find new research papers in:
• Health Sciences
• Earth Sciences
• Cognitive Science
• Mathematics
• Computer Science
• Academia ©2024

• Previous Article
• Next Article

I. INTRODUCTION

Ii. non-reactive potentials, a. lennard-jones (lj) potentials, b. class i force fields, iii. reactive potentials, a. embedded atom potentials [embedded-atom method (eam), modified embedded-atom method (meam)], b. bond-order potentials [bop, rebo, adaptive intermolecular rebo (airebo), mod-lj airebo, airebo-m], c. bond-order potentials with screening (rebo-s), d. charge-optimized many body (comb), iv. machine learning (ml), acknowledgments, review of force fields and intermolecular potentials used in atomistic computational materials research.

• Split-Screen
• Article contents
• Figures & tables
• Supplementary Data
• Peer Review
• Open the PDF for in another window
• Reprints and Permissions
• Cite Icon Cite
• Search Site

Judith A. Harrison , J. David Schall , Sabina Maskey , Paul T. Mikulski , M. Todd Knippenberg , Brian H. Morrow; Review of force fields and intermolecular potentials used in atomistic computational materials research. Appl. Phys. Rev. 1 September 2018; 5 (3): 031104. https://doi.org/10.1063/1.5020808

Download citation file:

• Ris (Zotero)
• Reference Manager

Molecular simulation is a powerful computational tool for a broad range of applications including the examination of materials properties and accelerating drug discovery. At the heart of molecular simulation is the analytic potential energy function. These functions span the range of complexity from very simple functions used to model generic phenomena to complex functions designed to model chemical reactions. The complexity of the mathematical function impacts the computational speed and is typically linked to the accuracy of the results obtained from simulations that utilize the function. One approach to improving accuracy is to simply add more parameters and additional complexity to the analytic function. This approach is typically used in non-reactive force fields where the functional form is not derived from quantum mechanical principles. The form of other types of potentials, such as the bond-order potentials, is based on quantum mechanics and has led to varying levels of accuracy and transferability. When selecting a potential energy function for use in molecular simulations, the accuracy, transferability, and computational speed must all be considered. In this focused review, some of the more commonly used potential energy functions for molecular simulations are reviewed with an eye toward presenting their general forms, strengths, and weaknesses.

Rapid advancements in computing power for the modeling of materials at the atomic scale have vastly expanded the array of phenomena that can be examined using computer simulations. Simulations are used in physics, chemistry, astronomy, material science, and biology to lend insight into physical and chemical processes as well as to provide predictions about properties and behavior. If a desired property or dynamic phenomenon is to be studied, before undertaking any simulations a method must be selected that has the ability to “capture” the behavior of interest and decisions need to be made about the desired level of accuracy of the results. When making these decisions, trade-offs between system size, i.e., the number of atoms, and the accuracy of the method are inevitable. For example, the highest level of accuracy can be obtained using quantum mechanical (QM) methods but because electronic degrees of freedom have computational cost that scales roughly as the cube of the number of these degrees of freedom, system sizes are limited to tens, hundreds, or perhaps thousands of atoms, depending on the level of theory. 1–5  

Despite their computational cost, quantum mechanical methods have been widely used for carefully selected problems with great success. For example, the temperature at which carbon nanotubes have a carrier density similar to an undoped metal was predicted using density functional theory (DFT). 6 Subsequently, these nanotubes were used in the construction of field effect transistors, providing experimental evidence for this prediction. 7 For more recent examples of the successes of density functional theory, the reader is referred to Cramer and Truhlar. 3  

Similarly, ab initio molecular mechanics, a hybrid approach that uses quantum mechanical methods to extract atomic forces from electronic structure calculations and integrates Newton's equations of motion to update the positions of atomic nuclei, is limited in scale to hundreds of atoms. 8 In contrast, classical MD 9–12 and Monte Carlo (MC) simulations 13–15 make use of approximate atomic and molecular interactions instead of undertaking electronic structure calculations. Advances in computational power and parallel computing have made million-atom simulations with explicit atomic interactions possible. 16–19 Methods to accelerate simulations or increase the system size, such as coarse graining, have also proven very useful in the simulation of polymers, 20 biological systems, 21,22 and plastic deformation of amorphous solids. 23 Despite the popularity and usefulness of coarse graining, in what follows the focus is on explicit atom simulations.

Regardless of whether all atoms are included, or coarse-graining methods are used, molecular simulations require the calculation of the energy of the system as a function of particle position. This enables the determination of the total energy of the system required for MC or MD simulations and the force acting on each particle for MD simulations. This is achieved by approximating quantum mechanical reality with analytic functions that represent the effective interactions between particles as a function of their spatial coordinates. These potential energy functions vary in their degree and manner of approximation and so also vary in their complexity and accuracy. Additionally, unique potential energy functions have been developed to model various materials and phenomena and, thus, have different strengths and weaknesses. Therefore, the choice of the potential energy function to be used in any simulation is governed by careful consideration of the processes and materials to be modeled, and the quantities that one hopes to analyze from the simulations. These considerations need to be balanced against available computational resources. An effective analytic potential will have the simplest functional form possible, while capturing the essence of the underlying quantum mechanical interactions. To best match reality, potential energy functions typically have several adjustable parameters.

Interatomic potentials are typically parameterized to reproduce quantities obtained from experiments or quantum mechanical calculations, with the choice of fitting parameters based on the ultimate end use of the potential. For example, potentials developed to model solids, e.g., metals and covalent materials, are typically fit to properties such as cohesive energies, elastic constants, lattice constants, and surface energies. For liquids, properties such as heat of vaporization and density are typically used to fit parameters for non-bonded interactions, and bonded parameters such as bond lengths, equilibrium angles, and force constants are matched to quantum mechanical calculations. Molecular properties, such as atomization energies and bond energies, have also proven very useful in the development of interaction models. To model chemical reactions, careful attention must be paid to fitting forces required to break bonds. A useful potential function will also have some degree of transferability, having the ability to describe structures not included in the fitting database, at least in a qualitative sense. This is all pursued while keeping in mind that an acceptable interaction model must be computationally efficient enough to be practically used to explore the systems of interest. This process of developing a potential energy function, sometimes referred to as an art as well as a science, requires a combination of chemical insight, trial and error, and tenacity on the part of its developer. 24  

In what follows, some of the commonly used interatomic potentials for computational materials research are reviewed, which are summarized in Table I . The discussion begins with nonreactive potentials typically used for simulations involving biological molecules, such as proteins. This is followed by a discussion of more sophisticated, reactive potentials for metals and covalent materials. Because the potentials discussed herein have been used extensively, this review paper does not attempt to discuss all of the published simulations that utilize these potentials. Rather, a few relevant example simulations are discussed and references to other relevant work are provided.

Comparison of strengths and weaknesses of several potential energy functions.

The Lennard-Jones (LJ) potential 25 is one of the simplest interatomic potentials, yet over the years it has found application in a multitude of simulations. Used as a stand-alone potential, or in conjunction with other potentials, it is responsible for describing general phenomena in a wide range of fields including simulations of liquids and polymers, and interactions between solids including friction and wear. The potential is particularly useful for simulation of adhesion and cohesion between non-reacting atoms, solids, and molecules. One general form of the LJ potential is as follows:

where r i j is the interatomic spacing between two atoms i and j . The parameters σ i j and ϵ i j represent the interatomic distance at an energy of zero and the minimum interaction energy, respectively. Thus, this potential is a two-body potential depending only on distance between pairs of atoms. The potential is made up of an attractive term that describes the van der Waals polarization energy and a repulsive term that comes into play only at short distances. As two atoms approach each other, the 1 r i j 6   term dominates and the LJ function is negative (attractive), with the depth of the attractive portion given by ε ij . As separation decreases further, the 1 r i j 12   term dominates and the potential switches from negative to positive at r ij = σ ij , with a steeply increasing repulsive wall at separations less than σ ij .

This LJ potential described above is sometimes referred to as a 12–6 potential denoting the exponents contained in the repulsive and attractive terms, respectively. The origin of the exponent 6 in the attractive term can be derived directly from first principles by considering the effects of polarization on electron density; however, determination of the repulsive exponent cannot be determined uniquely. While 12 is a common choice for the exponent and is used to describe intermolecular interactions in a number of potentials, other variants of the LJ potential exist, notably 9–3 which is commonly used to describe interactions between particles and a smooth wall. 26,27

The LJ potential has the distinction of being one of the first interatomic potentials used in MD simulations. In 1964, Rahman 28 presented one of the first simulation studies of liquid Ar. In that work, and in subsequent work by Verlet, 29 MD simulations using the LJ potential were compared to neutron scattering data, in part to provide proof of concept for this emerging scientific method. Remarkably, despite the fact that the potential is not fully tethered to quantum mechanical principles, this simple two-body potential proved capable of reproducing experimental results, such as the radial distribution function and self-diffusion coefficient, as a function of temperature and density. These early simulations also illustrated that as a first approximation, in the case of non-reacting atoms and molecules, many-body forces may be neglected. The methods employed to calculate properties, such as time correlations, self-diffusion coefficient, radial distribution functions, and structure factors, as well as other thermodynamics properties including temperature and pressure are the foundation of modern MD simulation.

Despite its simplicity, the LJ potential has proven to be extremely useful in the simulation of materials. For example, it was used to model interacting surfaces and has explained a number of fundamental questions in the field of tribology. 30–35 Amontons' law of friction states that the static friction coefficient between two contacting bodies is proportional to applied load and is independent of contact area. However, simulations of friction between bare surfaces show that friction vanishes except in cases of lattice commensurability as the area of contact grows. Müser and Robbins 30 provided a simple microscopic explanation for this deviation from Amontons' law and showed that an area-independent friction coefficient is produced for any surface geometry if an adsorbed layer of atoms is introduced into the interface. In their simulations, the atoms in the solid were anchored to their equilibrium positions via stiff springs, while the interactions between atoms in adjacent surfaces and adsorbed atoms were described using an LJ potential. This simulation showed that in the absence of absorbate atoms, friction disappeared as contact area was increased because the potential energy surfaces (PES) do not align and the two surfaces stay separated. When adsorbate atoms are present, they tend to sit in opposing valleys in the potential energy surface. In this configuration, the adsorbate atoms pin the two surfaces together and prevent sliding to produce a friction coefficient that is independent of contact area.

Lennard-Jones potentials coupled with MD simulations have also been used to test the applicability of continuum mechanics to nanoscale contacts. 32,33 It is clear that real materials, with surface roughness, anisotropic elastic properties, and non-continually varying strain fields, deviate from the assumptions of continuum mechanics. These simulations demonstrated that while atomic discreteness within bulk of the solids did not have a significant effect, atomic-scale surface roughness produced significant deviations from continuum theory with contact areas and stresses being potentially changed by a factor of 2. Specifically, Luan and Robbins 32 showed that spherical, amorphous, and stepped tips produce very different pressure distributions when in contact with a flat substrate (Fig. 1 ). Atomic-scale roughness and changes in geometry significantly impact the pressure distribution regardless of the presence of adhesion between the tip and the substrate. The most dramatic differences between continuum predictions and atomistic simulations were observed in friction and lateral contact stiffness, which changed by an order of magnitude. The above results have important ramifications for the interpretation of scanning probe microscopy experiments and potentially for macroscopic contacts.

(Top row) Spheres represent atoms in spherical, amorphous, and stepped tips. Pressure distributions for the tips without adhesion (second row) and with adhesion (third row). Reproduced with permission from B. Luan and M. O. Robbins, Nature 435 (7044), 929 (2005). Copyright 2005 Springer Nature.

In this section, non-reactive force fields commonly used for (bio)molecular simulations, whose flexibility enables them to be applied to a variety of systems, are discussed. These force fields, which use variations of the so-called Class I potential energy function, are commonly used to simulate a wide variety of materials, including proteins, surfactants, lipids, polymers, and other condensed-phase or aqueous systems. For the sake of brevity, only five widely used force fields are discussed: Amber, 36 Chemistry at HARvard Macromolecular Mechanics (CHARMM), 37,38 GROningen MOlecular Simulation (GROMOS), 39 all-atom optimized potential for liquid simulations (OPLS-AA), 40,41 and Transferable Potentials for Phase Equilibria (TraPPE). 42 The focus is on recent development and applications; more details of the history, development, and parameterization strategies of these force fields can be found elsewhere. 43,44

The functional forms of these potentials are similar and not based on any quantum mechanical rationale, with terms representing bond stretching, angle bending, rotation about dihedrals, improper dihedrals, and non-bonded interactions

Bond stretching, angle bending, and improper dihedrals are all modeled by harmonic functions. The first term in Eq. (2) is a summation over all bonded pairs of atoms, with b , b 0 , and K b representing the length, equilibrium length, and stiffness of each bond, respectively. The second summation in Eq. (2) is over all triplets of bonded atoms, i.e., atoms i , j , and k , where i is bonded to j and j is bonded to k . The angle formed by the two bond vectors is θ , with θ 0 and K θ representing the equilibrium angle and the angle stiffness, respectively. The dihedral term sums over quadruplets of atoms i , j , k , and l , with bonds between i and j , j and k , and k and l , and models rotation about the j-k bond with ϕ ⁠ , δ n ⁠ , and K ϕ , n representing the dihedral angle, phase of the dihedral angle, and a force constant, respectively. Because rotation about bonds is periodic, a cosine Fourier series is used. In the term that models the improper dihedrals, ϕ ⁠ , ϕ o ⁠ , and K ϕ represent the improper dihedral angle, equilibrium dihedral angle, and force constant, respectively. This function is used to maintain geometry at sites where an atom is bound to three other atoms (i.e., to maintain planarity around a central sp 2 hybridized atom). Not all force fields use a harmonic term to model improper dihedrals, some use the same periodic function as the proper dihedrals. The bond, angle, dihedral, and improper terms comprise the bonded interactions of the potential energy function. Coulomb's law models electrostatic interactions, with q i and q j representing the fixed charges on atoms i and j , and r ij their interatomic distance. Attractive dispersion and repulsive Pauli exclusion interactions are modeled by the LJ potential [(Eq. (1) ], which is sometimes referred to as the van der Waals term.

The parameters ε ij and σ ij depend on the atom types of both i and j . In practice, ε and σ are typically defined for each atom type, and mixing rules are used to determine the pairwise interactions. All five of the force fields discussed here use geometric combining rules for ε ij , that is, ε ij = ( ε i ε j ) 1/2 . The OPLS-AA and GROMOS potentials also use the geometric mean for σ ij , while Amber, CHARMM and TraPPE use the arithmetic mean, σ ij = ( σ i + σ j )/2. Other small differences between the functional forms of these potentials are discussed in the sections below. The Coulomb and LJ terms comprise the non-bonded portion of these potential energy functions. The summation includes all intermolecular interactions, as well as intramolecular interactions in which atoms are separated by at least 3 bonds.

One advantage of these types of force fields (with the exception of TraPPE) is that they typically include parameters for a variety of atoms, including carbon, nitrogen, phosphorous, oxygen, sulfur, and halogens, in a wide variety of chemical environments. Therefore, they are able to simulate complex biomolecules, e.g., proteins and nucleic acids, and can be used in applications such as drug discovery. 45 It is also relatively easy to add additional atom types and the relative simplicity of the function makes it feasible to use these force fields for simulations with hundreds of thousands or millions of atoms over time scales of microseconds. For example, Perilla and Schulten simulated a solvated human immunodeficiency virus (HIV-1) capsid, consisting of 64 × 10 6 atoms, for 1.2 μ s using the CHARMM36 force field. 16 One disadvantage of this class of force fields are that the system topology must be specified at the beginning of each simulation and bonds cannot be broken or formed. Therefore, chemical reactions cannot be modeled. While constant-pH methodologies have been developed to model proton titration, 46 other reactions cannot be studied. Another limitation is that atomic charges are fixed, so that redistribution of charge in response to a changing electric field cannot be captured. Polarizable force fields, such as the CHARMM Drude polarizable force field, 47 show promise, but are more limited in terms of available atom types, and are more computationally demanding. Finally, this class of force fields is generally targeted to liquid-phase systems. While they are suitable for a variety of liquid and aqueous systems, Amber, CHARMM, OPLS, and GROMOS are generally unable to model solids and are not parameterized to reproduce gas-phase data. In contrast, TraPPE is fit to vapor-liquid equilibrium (VLE) data. Therefore, it is better able to simulate gases and predict critical points, but it is more limited than the other force fields in terms of available atom types.

Originally developed for use in the Assisted Model Building with Energy Refinement program, 36 the Amber force field has the functional form given in Eq. (2) . Many protein parameter sets exist, detailed discussion of which is outside the scope of this review. In the widely used ff99SB parameter set, 48 high-level ab initio quantum mechanical calculations were used to improve the existing protein backbone dihedral parameters. Side-chain parameters were updated and validated using microsecond-timescale simulations in the ff99SB-ILDN (where ILDN corresponds to the one-letter codes of the four side-chains whose parameters were modified) force field. 49 Recently, backbone dihedral parameters were adjusted and side-chain dihedral parameters were completely refit in ff14SB. 50 In addition to proteins, there are parameter sets for nucleic acids, 51–55 carbohydrates, 56 and lipids, 57 as well as the general Amber force field (GAFF) for organic molecules. 58 GAFF is compatible with the other Amber force fields and has parameters suitable for modeling most organic molecules composed of H, C, N, O, S, P, and halogens. In addition to drug design, it can be used in simulations of large molecules, such as polymers. Baker et al. 59 simulated supramolecular polymers (polymers whose monomers are not covalently linked) consisting of 1,3,5-benzenetricarboxamide (BTA) units and found that small changes in BTA molecular structure (presence or absence of a stereogenic methyl group) result in self-assembled fibers that have almost identical structure but whose equilibrium dynamics differ by an order of magnitude. MD simulations using GAFF revealed that the fibers consisting of chiral units had increased internal order and hydrogen bonding compared to the achiral fibers (Fig. 2 ). 59  

Summary of simulation results outlined in Ref. 59 by Baker et al. that examine achiral (1) and chiral (2) polymer fibers in water. Starting system and equilibrated system are shown in (a) and (b), respectively. Panel (c) demonstrates that the water solution is able to penetrate the core of the polymer. Simulated (lines) and experimental (points) small-angle X-ray scattering results are compared in (d). Radial distribution functions (g(r)) for polymer cores (e) show the chiral system is more ordered.

The CHARMM (Chemistry at HARvard Macromolecular Mechanics) potential energy function uses the functional form of Eq. (2) with a Urey-Bradley term applied to the distance between the terminal atoms in a limited number of angles to improve agreement with vibrational spectra. 60 This term takes the form

where r 1,3 is the distance between the 1,3 atoms, r 1,3 ; 0 is the equilibrium distance, and K U B is the force constant. For simulations of proteins, a grid-based energy correction map (CMAP) is applied to the protein backbone ϕ and ψ dihedral angles to improve secondary structure sampling. 61,62

The current version of the force field, designated CHARMM36, has parameters for proteins, 63 nucleic acids, 64–66 carbohydrates, 67–69 lipids, 70 and general organic molecules (the CHARMM general force field, or CGenFF). 71 Additional developments in the CHARMM force field can be found in a recent review. 47 As with Amber, the various CHARMM parameter sets are designed to be compatible with each other, enabling simulation of heterogeneous systems, such as proteins in a phospholipid membrane. For example, Dror et al. simulated the M2 muscarinic acetylcholine receptor, a prototypical G-protein-coupled receptor (GPCR), embedded in a hydrated palmitoyloleoylphosphatidylcholine bilayer. 72 Ligands were placed in the aqueous part of the simulation box, away from the receptor, and allowed to freely diffuse before spontaneously binding. The simulated binding modes, which were validated experimentally, provide a structural basis for the rational design of drugs targeting M2 and other GPCRs. 72  

The all-atom optimized potential for liquid simulations (OPLS-AA) also uses the functional form of Eq. (2) with a four-term Fourier series to model dihedrals 40,41,73

where V 1 ,   V 2 ⁠ , V 3 ⁠ , and V 4 are fitting parameters. This potential can be used to model a wide variety of organic molecules, and additional parameter sets specifically for carbohydrates 74 and proteins 73 have been published. Recently, high-level quantum chemical methods were used to develop new peptide dihedral parameters (OPLS-AA/M), 75 and the OPLS3 force field for drug-like small molecules and proteins was released. 76 Other groups have also published parameters aimed at improving the description of phospholipids. 77,78

The original OPLS-AA parameters for hydrocarbons were unable to reproduce the properties of long-chain alkanes, such as dodecane, because they were developed using short alkanes. Siu et al. modified torsion parameters, partial charges, and the methylene hydrogen LJ well-depth to match experimental densities and heats of vaporization. 79 This parameter set, termed L-OPLS, had improved properties for long-chain hydrocarbons and was later extended to alcohols, esters, and monoolein bilayers. 80 This revised potential was recently used to examine the physical properties at ambient conditions for binary surrogate fuel mixtures of n -hexadecane with n -alkylbenzenes, with n-alkyl chains up to four carbon atoms. 81 Benzene and toluene mixtures with n -hexadecane had an interesting minimum in isoentropic bulk modulus as a function of n -alkylbenzene content not present in the mixtures containing the longer chain n -alkylbenzenes. Angular radial distribution functions demonstrated that the minimum in the bulk modulus data was linked to molecular-scale packing differences in mixtures with short-chain alklybenzenes compared to mixtures with the longer chain n -alkylbenzenes.

GROMOS differs from the force fields discussed above in that it uses a united-atom representation of aliphatic CH, CH 2 , CH 3 , and CH 4 groups. That is, a single interaction site is used to represent these groups, with no explicit hydrogen (EH) atoms. As with other nonreactive force fields, it has been applied to a wide range of biomolecules, including proteins, nucleic acids, 82 carbohydrates, 83 and lipids. The most recent parameter set is denoted as 54A8. 84  

The Transferable Potentials for Phase Equilibria (TraPPE) family of force fields, unlike the four previously discussed, are not geared toward biomolecular simulations at ambient conditions. Rather, TraPPE uses phase equilibrium data, typically the vapor-liquid coexistence curve, to fit intermolecular interaction parameters. This enables it to accurately reproduce properties over a wide range of temperatures and state points. TraPPE also differs from the above force fields in that bond lengths are fixed rather than modeled with a harmonic potential. The initial united-atom version, designated TraPPE-UA to distinguish it from the other TraPPE variants, began with parameters for normal alkanes, 42 then branched alkanes, 85 alkenes, and alkylbenzenes. 86 It was later extended to include a variety of molecule types, including alcohols; 87 ethers, glycols, ketones, and aldehydes; 88 amines and amides; 89 thiols, sulfides, disulfides, and thiophene; 90 and five- and six-membered cyclic alkanes and ethers. 91 The TraPPE-Explicit Hydrogen (TraPPE-EH) force field 92 includes interaction sites for hydrogen atoms. TraPPE-EH was originally developed for n-alkanes 92 and also includes parameters for benzene and five- and six-membered heterocyclic aromatic compounds, 93 and substituted benzenes and polycyclic aromatic compounds. 94 The TraPPE force fields are generally more accurate than Amber, CHARMM, OPLS-AA, and GROMOS in predicting vapor-phase properties and critical points due to the different fitting strategy. However, unlike those force fields, TraPPE cannot model complex biopolymers, such as proteins and polysaccharides.

The embedded-atom method (EAM) first detailed in the original work of Daw and Baskes 95 is a many-body interaction potential that is particularly well suited to the modeling of metals. The term “embedded” refers to the idea that each atom (technically, each ion because each atom donates valence electrons to the metal) is embedded in the environment established by all other atoms in the system. That environment is comprised of two distinct parts: (1) atoms of the host material interact pairwise with the embedded atom; and (2) valence electrons that are shared throughout the metal interact with the embedded atom in a manner that is fundamentally many-bodied in nature (meaning, not reducible to an equivalent set of pairwise interactions).

The total energy for a system containing N atoms using the EAM potential is

where F i is the energy required to embed atom i at its location, ρ h , i   is an approximate local electron density built by adding contributions from all atoms in the local region except the atom being embedded, and ϕ i j is the core-core repulsion between atoms i and j . The manner in which contributions from non-localized valence electrons of the metal are handled is what distinguishes the EAM from other many-body potentials. The forms of the many-body terms often address specific kinds of bonds and sequences of bonds, so in a sense one still views interactions as pertaining directly to a particular sequence of atoms. In contrast, the EAM potential makes use of the local arrangement of atoms in the embedded atom's neighborhood as a means to approximate the local charge density of valence electrons that are no longer associated with their parent atoms. That approximation to the local density is then converted to a contribution to the potential energy through an embedding function. This may seem to be a subtle distinction, but it is important to keep in mind that to the lowest-order, the valence electrons within a metal behave as a free-electron gas; the local arrangement of atoms is being used as means to approximate small local deviations from the density associated with the free-electron gas.

For MD simulations, the EAM is simple to implement and because it can use the same kinds of cutoffs used to manage long-range pairwise interactions, such as LJ and Morse potentials, it does not come at much additional computational cost. Also, by making use of cubic splines and tabulated data, extensive use can be made of experimental data and ab initio calculations (predominantly via DFT) to create optimized potentials that address a wide range of bulk and surface properties for many face-centered-cubic (fcc) metals including alloys. 96–100 EAM potentials were subsequently developed for body-centered-cubic (bcc) 18,98,101 and hexagonal-closed-packed (hcp) metals. 102,103

Very closely related to the EAM is the modified embedded-atom method (MEAM), which builds a directional dependence into how the electron density at the site of the embedded atom is calculated. While originally developed for fcc metals, 104 potentials for hcp 105 metals have also developed. The original MEAM approach 104 works with only first-nearest neighbors (NN), while later work by Lee 106 also considers directional effects with second-nearest neighbors (2NN MEAM). MEAM potentials with cutoffs between second and third neighbors for bcc metals and the fifth and six neighbor distances for fcc potentials have also been developed. 107 In principle, such modified approaches open up the range of materials to which the embedded atom formalism can be applied but, in practice, such directional considerations are more commonly associated with covalent bonding potentials. For example, Baskes developed a potential for Si within this formalism. 108 The resulting MEAM potentials are more fully empirical and less physically motivated. 104 This raises the question of how appropriate such an approach is when typically the large-scale systems of interest in many MD simulations are at some level removed from the smaller systems used to help fit potential parameters. Consequently, the MEAM is not necessarily meant to be taken as better, rather just different and possibly appropriate in certain contexts.

Additional variants of these embedded-atom models include the analytic EAM 109–111 and analytic MEAM for bcc metals. 112 Johnson's simple analytic fcc EAM was developed to study the relationships between the model parameters and system energies. 109 In that work, the average shear modulus was shown to play a dominant role in determining vacancy formation energies and the energies of low-index plane surfaces due to the mathematical form of the embedding function and its slope at the equilibrium electron density. It was also shown that the embedding functions used in several models are essentially equivalent. The MEAM for bcc metals utilizes a new pair potential. As a result, this potential is able to reproduce the shear moduli of all bcc transition metals. 112  

A very useful study of the general kinds of behavior that one can capture with the EAM, as opposed to a strictly pairwise potential, was conducted by Holian. 113 The context was the study of plastic flow under a variety of conditions. The manner of implementation was intentionally kept very simple with the focus being on establishing the qualitative differences that can result with the EAM compared to pair potentials. The study clearly demonstrates that the EAM enhances plastic flow and suppresses some of the behavior associated with brittle materials that are seen when using a pair potential. The key underlying reason for this is that the EAM captures the essential feature that the energy associated with creating a vacancy is considerably less than the cohesive energy per particle for the material as a whole. A full description of Holian's study is contained in Rapaport 114 with the added benefit that Rapaport includes a fully coded implementation. In what follows, we discuss Holian's chosen form of the EAM potential as a way of introducing the basics of the EAM formalism. For an N atom system, E t o t is given by

The factor χ specifies fractionally the contribution that the pair potential gives to the cohesive energy (per atom) of a close-packed solid in the nearest-neighbor approximation; this cohesive energy is expressed as

where d is the dimensionality of the system with the factor d d + 1 representing the number of nearest neighbors in a close-packed solid, and ϵ is the contribution associated with a single nearest-neighbor pair. So that all choices for χ yield this cohesive energy, parameters must be set to give this cohesive energy for both χ = 1 (only pairwise Lennard-Jones interactions) and χ = 0 (only many-bodied interactions via the embedding function). The use of the fractional parameter χ is specific to Holian's study: different choices for χ give different interaction potentials that can be directly compared against one another.

The pairwise only case ( χ = 1 ):   ϕ L J r i j is the LJ pair potential up to a radius r s p l beyond which a cubic spline is used to smoothly turn off the potential at radius r m a x ( ⁠ r s p l and r m a x are not independent and so only one can be freely set). The cohesive energy for the close-packed solid is obtained by setting all nearest neighbors at the minimum energy radius r = r 0 = 2 1 6 σ   where the LJ potential evaluates to − ϵ ⁠ .

The many-body only case ( ⁠ χ = 0 ⁠ ): U ρ i is the non-linear embedding function that describes the energy associated with embedding atom i to a location where the electron density has a value of ρ i ⁠ . The electron density ρ i is viewed as having attained its value by contributions from each atom in the local vicinity of the embedding site. Valence electrons contributed to the metal from each atom are not localized and so, in principle, contributions could be taken from every atom in the material, but practically speaking it is reasonable to only treat a local neighborhood as is done with pairwise interactions. To smoothly implement such a cutoff, each host atom contributes to ρ i through a weighting function w r i j

Here, e is the base of the natural logarithm. For a close packed solid with all nearest neighbors at r = r 0 ,   this gives a density ρ i = 1 e ⁠ . The non-linear embedding function is chosen to be

Substituting the density for the case of the close-pack solid, ρ i = 1 e ⁠ , does give the correct cohesive energy. With both χ = 0 and χ = 1 giving the correct cohesive energy, any value 0 ≤ χ ≤ 1 can be freely explored without adjusting any other parameters. Values in the range 0.2 ≤ χ ≤ 0.4 are appropriate for most metals with the interpretation that χ represents the ratio E v a c E c o h ⁠ , where E v a c is the energy needed to create a vacancy.

The above model nicely encapsulates the basic formalism of the EAM. An approach to create highly optimized potentials for specific metals opens up the functions ϕ r ⁠ , w r ⁠ , and U ρ as functions to be fitted through the adoption of as many parameters as is appropriate for the variety of properties for which the potential is to be optimized. For instance, in the work of Sheng, 100 potentials for 14 fcc elements were developed, and in each case, 36 parameters were used to fit the functions ϕ r ⁠ , w r ⁠ , and U ρ ⁠ . The wide range of properties for which the fits were optimized included lattice dynamics, mechanical properties, thermal behavior, energetics of competing crystal structures, defects, deformation paths, and liquid structures. As is the case generally with fitting parameters, caution must be exercised when using highly optimized potentials in contexts removed from the training environment in which the fits were developed. A large part of the attraction of the EAM is the ease and directness with which it connects to DFT calculations, and so the hope is that EAM potentials can confidently be employed in large-scale non-equilibrium molecular dynamics simulations for which ab initio approaches still remain far out of reach.

More recent work has focused on the development of a MEAM potential for complex materials such as metallic alloys 115 and even saturated hydrocarbons. 116 Subsequently, the MEAM parameters for saturated hydrocarbons were updated and a bond order (MEAM-BO) term was added to describe the energetics of unsaturated hydrocarbons, i.e., double and triple carbon bonds. 117 Unfortunately, while the MEAM formalism is a flexible formalism, a multi-component potential does not result solely from a suitable combining of single-element potentials. In the MEAM approach of Jelinek, 115 nine adjustable parameters are introduced for each two-element pairing from the single elements aluminum, silicon, magnesium, copper, and iron. The formalism does extend through the use of screening parameters to systems of any combination of the five chosen elements, but the focus of this work is on the rigorous testing of all the possible binary systems from this set of single elements. The recent work of Nouranian 116,118 demonstrates that the MEAM formalism can be used to model saturated hydrocarbons with an accuracy comparable to traditional hydrocarbon potentials such as reactive empirical bond-order (REBO) 119,120 and ReaxFF. 121 This opens the door to the possibility of using the MEAM framework to model complex systems involving metals and organics.

Recently, the MEAM potential and MD simulations were coupled with phase-field crystal model to investigate the solid-liquid properties of Fe. 122 In that work, a second nearest-neighbor parameter set for the MEAM potential was developed for Fe. This parameter set led to predictions of melting point, latent heat, and the expansion in melting that were very close to the experimental values. In addition, properties such as elastic constants, surface free energies, and vacancy formation energies were also well reproduced. These simulations were used to model the liquid-solid phase boundary (Fig. 3 ).

Snapshot of an MD simulation using the MEAM potential of solid-liquid coexisting structure for Fe in the ⟨001⟩{110} orientation. Spheres represent atoms that are colored based on their order parameter with blue and red representing solid and liquid, respectively. The transition from red to blue represents the solid-liquid interface. Reproduced with permission from Asadi et al. , Phys. Rev. B 91 (2), 024105 (2015). Copyright 2015 American Physical Society.

One of the first potential energy functions that attempted to model covalent systems was the Abell bond-order potential (BOP), which stipulates that the binding energy of a system may be calculated by adding the attractive and repulsive interactions between neighboring atoms in a system 123  

In this equation, the attractive E A and repulsive E R terms are modified as a function of the distance between atoms i and j and the cutoff function f c r i j has a value of zero if the bond distance r i j is greater than r m a x and one if it is less than r m i n ⁠ . This function is designed to have the energy smoothly go to zero beyond r m a x ⁠ . For values in between r m a x and r m i n ⁠ , the cutoff function is given by

The attractive and repulsive functions in Eq. (11) are analytic functions that account for all electronic contributions to the potential energy, without performing quantum mechanical calculations. Using appropriate attractive and repulsive functions allows the calculation of bulk-phase properties much more quickly than corresponding quantum mechanical calculations. The function b ij is the bond-order term, which measures the number of bonds between neighboring atoms. Originally, the bond order was calculated from higher-level quantum mechanical calculations. Tersoff developed a formalism that allows the bond order to be calculated from parameters, which includes the number of neighbors of each atom as well as the bond angles between atoms, and changes based on the local bonding environment. 124–127 Tersoff's work in parameterizing the bond-order term was then applied to systems composed of hydrogen and carbon, 119,128,129 hydrogen and oxygen, 130 silicon and hydrogen, 131 silicon and carbon, 132–134 and other covalently bound systems. 135,136

Analytic potentials with the highest degree of transferability are based on sound quantum-mechanical bonding principles. 24 The BOPs have been shown to be based on the second-moment approximation. 24,137 According to the moments theorem, the n -th moment of the local density of states on an atom i is determined by the sum of all paths between n neighboring atoms that start and end at atom i . A good estimate of the bond energy can be obtained by knowing only the second moment, which is related to the energy beginning on atom i and summing over the nearest neighbors. 137,138 Thus, the local electronic bond energy for each atom from molecular orbitals is proportional to the square root of the number of neighbors. This approximation was used to derive an expression for the energy of a solid in terms of pair-additive interactions between atoms. 138 The resulting expression can be rearranged into an expression that is similar to the bond-order potentials. 139 Thus, the bond-order expressions of the BOPs capture the essence of quantum-mechanical bonding. It is also possible to derive other analytic forms based on higher moments. For example, Pettifor and coworkers have developed a BOP that contains terms for σ and π bonding from the tight-binding model using the moments theorem. 140–143  

Working independently, Beardmore and Smith 144 and Dyson and Smith 132 made use of Brenner's hydrocarbon reactive empirical bond-order (REBO) potential 119,129 to develop a bond-order potential for C-Si-H. Subsequently, one of these potentials was modified slightly to improve interactions between organic molecules and Si surfaces. 145 Because these potentials utilize the REBO formalism, they inherit its strengths and weaknesses. The REBO potential and the BOP of Tersoff 125 are unable to reproduce the elastic constants of diamond and graphite; as a result, the C-Si-H potentials based on this formalism also have poor agreement with experimentally determined elastic properties. 146,147 In an effort to correct the inability to reproduce the elastic constants and other discrepancies, the functional form of the REBO potential was updated to create the second-generation REBO potential, or REBO2. 24,120 A larger database of molecules was included in the parameterization of the REBO2 potential, leading to improvements in the energies of hydrocarbon molecules, elastic constants, interstitial defect energies, and surface energies of diamond.

To allow the REBO2 potential to model chemical reactions during a simulation, the bond-order term, b ij , can be calculated based on the local coordination in which an atom is found in and is given by

The values for b i j σ π and b j i σ π   are different due to the fact that atoms i and j may be in different bonding environments in the system. The terms π i j R C and π i j D H   are sometimes summed and written as b i j π ,   which comes from determining if an atom is found in a conjugated system and if there is a dihedral angle present in carbon-carbon double bonds. This REBO2 potential was developed to model both molecular and solid-state hydrocarbon systems. Comparisons of bulk-terminated diamond surfaces show that the first-generation REBO potential incorrectly predicts the energy per atom of the π-bonded surface in a reconstructed diamond (111) surface compared to DFT results, while the REBO2 potential agrees with DFT results and shows good agreement in bond lengths and force constants for a wide range of hydrocarbon molecules for which it was not specifically parameterized, including conjugated and branched hydrocarbons. In addition, interstitial defect energies and surface energies for diamond are fairly well reproduced. Despite being able to reproduce the zero-Kelvin elastic properties of diamond and graphite, subsequent work demonstrated that the elastic constants of diamond as a function of temperature were not well reproduced. 148  

Recent work utilizing the REBO2 potential has focused on the fracture toughness of graphene 149 and carbon nanotubes. 150 The friction of self-mated diamondlike carbon (DLC) contacts as a function of hydrogen content was examined by Schall et al. using MD and the REBO2 potential. 151 Decreasing the hydrogen content was shown to increase the friction and the formation of interfacial covalent bonds. In addition to the formation of interfacial bonds when sliding against hydrogen-free films, sp 3 -hybrized carbon atoms in the upper DLC surface were converted to sp 2 -hybrized atoms as a result of the sliding (Fig. 4 ). This change in hybridization propagated 5 to 6 Å into the upper DLC surface during the 200 ps of sliding.

Coordination number for carbon atoms as a function of position normal to the loading axis ( y -axis) and MD simulation time ( x -axis) for the hydrogenated (a) and non-hydrogenated DLC films (b) in sliding contact with a mostly sp 3 -hybridized DLC counterface. The normal load is 60 nN on the upper DLC counterface. The interface between the two DLC surfaces is initially at a z coordinate of approximately 36 Å. The coordination number of the carbon atoms within the DLC is indicated by the color bar. A coordination number of 4, 3, and 2 corresponds to sp 3 -, sp 2 -, and sp-hybridization, respectively. Reproduced with permission from Schall et al., J. Phys. Chem. C 114 (12), 5321–5330 (2010). Copyright 2010 American Chemical Society.

Additional atom types, including C, O, and H, 152,153 C, F, and H, 154 and Si, 147 have also been parameterized within the REBO2 potential. A new BOP for systems containing Si, C, and H (2B-SiCH) was also recently parameterized. 155 This 2B-SiCH potential makes use of the C-C, C-H, and H-H parameters in REBO2 and the previously published Si parameters. 147 The third atom type, Si, was integrated into REBO2 via the use of tricubic splines, which provide maximum flexibility in the fitting with the result being only a 3.3% root-mean-square (RMS) error in the predicted energies of 105 Si-C-H containing molecules. The errors in the predicted energies do not increase with increasing molecule size as they do with the Si-C-H potential of Sbraccia et al. based on the first-generation REBO formalism that utilized two-dimensional spline functions. 145 It should also be noted that a number of empirical potential energy functions for silicon have been developed. Various formalisms have been used including the BOPs of Tersoff, 124,125,156 the Stillinger-Weber (SW) potential, 157 and the MEAM potential. 108,158 Balamane et al. 159 performed a detailed comparison of six Si potentials, including the most popular SW and Tersoff potentials. They concluded that each of these potentials has strengths and weaknesses and none are totally transferable. Subsequently, the environment-dependent interatomic potential (EDIP), 146 which was fit to ab initio data, and a BOP based on tight binding 160 were both developed with the aim of being more transferable and both do a much better job reproducing elastic properties and the Cauchy pressure of diamond-cubic Si.

While the REBO2 potential was sufficient for studying covalently bound solid-state systems such as diamond, graphite, and carbon nanotubes, the potential lacks a treatment for non-bonded interactions. This meant that simulations of liquid hydrocarbons and the intermolecular forces in graphite would not be sufficiently modeled using REBO2, or any REBO-type potential. Stuart et al. added long-range intermolecular interactions to REBO2 through the use of an adaptive LJ term to create the adaptive intermolecular reactive empirical bond-order (AIREBO) potential. 161 In addition to intermolecular interactions, torsional interactions, which explicitly model the rotation around single bonds, were added to the AIREBO potential. The total energy of the system in the AIREBO potential can be written as

For the long-range interactions, E L J − AIREBO is an adaptive LJ potential discussed below, E REBO is given by Eq. (11) , and the torsional potential E tors is given below. 161  

Because the AIREBO potential models bond breaking and formation, the inclusion of intermolecular interactions requires a switching function that modulates the strength of the intermolecular interaction, effectively turning it “off” and “on” according to the environment of the atom. This function changes based on three criteria: the distance between atoms, S d i s tan c e ⁠ , the strength of the bonding interaction between the atoms, S bond ⁠ , and whether or not both atoms are part of a network of bonds in the same system, S connectivity and is given by 161  

where E L J is a 12–6 LJ potential [(Eq. (1) ].

Because a particular atom may have a change in one of these three switching function values over the course of a simulation, checks are made each timestep of the simulation, allowing the intermolecular interactions to change adaptively. Each one of these switching functions may turn off the LJ interaction partially or entirely. There are parameters for carbon-carbon, carbon-hydrogen, and hydrogen-hydrogen interactions that allow the switching function to vary based on the type of atoms in a bond.

The torsional interactions within the AIREBO formalism are calculated using a cosine-power series, which produces a three-fold symmetric function for a bond between identical sp 3 carbon atoms, such as in ethane. This torsional potential is written as

where ω is the torsional angle and ε ′ is the height of the energy barrier to rotation. The need for the unique form of this torsional potential arises from the ability of the potential to model chemical reactions. The AIREBO potential has been used for simulations where intermolecular forces are important in maintaining the structure of hydrocarbon systems, such as simulations of graphene, 162 the chemisorption of hydrocarbon molecules on carbon nanotubes, 163 the compression of C 60 , Ne, or CH 4 filled carbon nanotubes, 164 and the friction of self-assembled monolayers. 165–168  

Recently, Vahdat et al. performed complementary tapping-mode atomic force microscopy (AFM) experiments and MD simulations that examined wear mechanisms between a DLC tip and an ultrananocrystalline diamond (UNCD) substrate. 169 Finite element (FE) simulations were used to determine the optimum tip shape and size for the MD simulations that would accurately capture deformations, stresses, and the contact area of the experiment. Repeated indents of the DLC and the UNCD were simulated and force curves were generated (Fig. 5 ) and atomic-scale wear mechanisms (Fig. 6 ) identified. Repeated indentation results in a hysteresis in the force curves, which is the result of the formation of covalent linkages between the tip and the surface (Fig. 5 ). Moreover, the sharp discontinuities or step-like structures apparent in the pullback curves correspond to the breaking of these covalent linkages upon pullback. By examining interfacial bond formation in the MD simulations, insights into the wear of DLC observed in the tapping mode AFM were obtained. Example bond formation events that result from the first, second, and third indents are depicted in Fig. 6 . Broadly speaking, interfacial chemical bonds are formed between unsaturated carbon atoms that are in close proximity as a result of the indentation. It should be noted, however, that proximity of unsaturated atoms is necessary for bond formation but not sufficient. Recent simulations demonstrate that there is a random component to bond formation. 170 Progressive indent and pullback cycles appear to cause increased roughness of the surface of the DLC.

( Left ) Close up of a DLC tip ( upper surface ) pulling back from a UNCD surface ( lower surface ) after contact. The maximum load was 50 nN before pullback. Carbon atoms are large grey and blue spheres and hydrogen atoms are small white spheres. Blue spheres represent atoms in the grain boundary of the UNCD. (Reproduced with permission from Vahdat et al. ACS Nano 8 (7), 7027––7040 (2014). Copyright 2014 American Chemical Society.) ( Right ) Force on the DLC tip as a function of separation distance for three indent and pullback cycles of the system shown at left . Arrows demonstrate force discontinuities that correspond to the severing of covalent linkages between surfaces.

Close up of interface between the DLC tip and UNCD surface before and after the first [(a) and (b)], second [(c) and (d)], and third [(e) and (f)] indent and pullback cycles. Color coding is the same as Fig. 4 except for red and green carbon atoms, which represent atoms participating in the chemical reactions. Arrow and circles show atoms that form a covalent linkage upon compression. (Based on work described in Ref. 169 .)

The AIREBO potential has the same values of σ and ε for all carbon-carbon interactions, regardless of the hybridization of the carbon atom or the type of atoms connected to it, e.g., the carbon atom in a –CH 2 - versus –CH 3 group. This means that no distinction is made for the chemical environment of the carbon atoms, i.e., the interaction between two sp 3 – hybridized carbon atoms is the same as the interaction between two sp-hybridized carbon atoms. To add more functionality and accuracy to the AIREBO potential, values of σ and ε for carbon atoms were modified to account for the hybridization of the carbon atom, i.e., the number of atoms bound to the carbon of interest. In other words, the values of σ and ε are now functions of the number and type of atoms bound to a particular carbon atom. The equations describing the new values for these parameters are

In these equations, the values of c k l σ and c k l ε represent particular coefficients in a rectangular region of space bound by four endpoints, N m i n C ,   N m i n H ⁠ , N m i n C ,   N m a x H ⁠ , N m a x C ,   N m i n H ⁠ , and N m a x C ,   N m a x H ⁠ . A bicubic spline is used to generate a smooth pathway between each endpoint. Additionally, the use of a bicubic spline allows the intermolecular interactions to be calculated during chemical reactions, where the coordination number is not a single, integer value. More information on the calculation of these adaptive intermolecular parameters and the spline functions for σ and ε can be found in the work by Liu and Stuart. 171 The use of these adaptive parameters based on coordination number allows this modified LJ portion of the AIREBO potential, or the mod-LJ AIREBO, to more accurately determine densities and enthalpies of vaporization of a wide range of organic molecules containing carbon atoms with various hybridization states. The mod-LJ AIREBO potential was also used to model mixtures of hydrocarbons in a study of molecules that may represent suitable renewable fuels, 172 with calculated densities and enthalpy of vaporization values comparing favorably to experiment.

Stuart and coworkers also recently took the step to make calculations of the bond order consistently bond-centric, as opposed to being a mixture of bond- and atom-centric approaches. 173 In the AIREBO potential, the b i j σ π and b j i σ π terms [Eq. (13) ] are calculated based on the types of neighboring atoms found around atoms i and j , where

where g i θ j i k   is a three-body term that is a function of the angle between bonds i-j and i-k. , the e λ j i k is a correction term added to improve the potential energy surface for the abstraction of hydrogen from diamond surfaces, and ω i k is a bond-weighting term. The P i j term is a bicubic spline that serves as a correction term to the bond order. This definition only looks at the types of atoms bound to atom i . To modify this equation and make it more bond-centric, the updated definition is

The terms N ¯ i j C and N ¯ i j H are now calculated as an average of the coordination numbers from both sides of a particular ij bond: N ¯ i j C = 1 2 N ¯ i j t + N ¯ j i t ⁠ , whereas the old calculation was single-sided based on atom i . The use of this new bond-centered approach provides a lower RMS error in energies when compared to DFT calculations for a variety of hydrocarbons, including structural isomers of C 28 and C 20 . 173  

While the LJ potential that models the intermolecular interactions works well at ambient pressures, it has been shown to be insufficient for reactions taking place at high pressures. 174 To better model systems at high pressures, the AIREBO-M potential models the long-range interactions using a Morse potential instead of an LJ potential, which is known to be too repulsive at close interatomic distances. 175 The new intermolecular Morse potential is written as

In this equation, ε i j and r i j e q describe the minimum energy and location of the potential energy well, while the parameter α adjusts the shape of the potential energy well at the minimum energy. This potential retains the same cutoffs for the intermolecular interactions as the original AIREBO potential, with the α term providing the ability to modify the shape of the potential energy well, and therefore the steepness of the repulsive wall, based on the types of atoms that are interacting. The AIREBO-M potential has much better agreement with experiment for the pressure dependence of the interlayer separation in a graphene bilayer than the original AIREBO potential. 175 Additionally, the lattice parameters in orthorhombic polyethylene produced using the AIREBO-M potential compared favorably with the experimental values.

In an effort to save computational time and because the REBO2 was developed originally to model nonpolar, covalent materials, a switching function was used in the REBO2 potential to cut off the potential between 1.7 and 2.0 Å for carbon-carbon interactions. 120 In their study of grain boundary fracture in diamond, Shenderova et al. 176 concluded that this switching function was problematic as carbon-carbon bonds are stretched beyond the 1.7 Å transition because it significantly influenced the forces in the vicinity of the inflection point, around 1.85 Å, in diamond for the ⟨111⟩ direction. To avoid this problem, they extended the C-C cutoff distance far beyond the inflection point. Later, Mattoni et al. examined brittle fracture in elemental and group-IV materials. 177 Not surprisingly, they also concluded that short-range potentials, such as Tersoff, 124–127,156,178 Stillinger-Weber, 157,179 and EDIP, 146,180 are not able to reproduce fracture-related properties due to their short-range cutoffs, which result in the overestimation of the forces to break bonds.

Subsequently, Pastewka et al. 181 examined fracture-related properties in single crystal diamond and carbon nanotubes using the REBO2 potential. As expected, the short-range nature of the hydrocarbon REBO2 caused the force to break (and conversely form) covalent bonds to be significantly larger than that predicted by DFT calculations, which resulted in non-physical crack propagation behavior (Fig. 7 ). This is a natural consequence of matching the bond energy to the experimental value while imposing a short-range cutoff. This shortcoming is shared by all bond-order potentials based on the REBO-style formalism with short-range cutoffs. 132,145,147,155,182,183 To address this problem, Pastewka et al. 181,184 altered the form of the cutoff function, f c r i j ⁠ , in Eq. (12) so that first nearest-neighbor interactions could be dynamically adjusted depending upon the local atomic environment of the bond. In the screened REBO potential, or the REBO-S, f c r i j takes the form

where the function f c 12 r i j is identical to Eq. (12) , f c 34 r i j is an additional cutoff function, and S i j is the total screening function. This screening function was originally introduced by Baskes et al. 185 and later by Kumagai et al. 186 The screening function uses an additional empirical function C i j k to screen bonds 181,184,187 and is given by

where C m a x and C m i n are parameters. 181 The screening function works by identifying atoms in the vicinity of a bond. If the atom is close to the bond, the bond is screened; if it is far away, the bond it is unscreened. This is in contrast to traditional methods of counting atoms within a distance of a bond. The resulting REBO-S potential is able to reproduce the correct crack propagation behavior in diamond (Fig. 7 ) and brittle fracture of carbon nanotubes.

Propagation of a 110 crack through a diamond (111) plane. ( Left ) Propagation simulated using the REBO2 potential, which shows tip blunting, while the REBO-S potential ( right ) produces more realistic propagation behavior. Figures provided by Lars Pastewka based on work outlined in Ref. 181 . Reproduced with permission from Pastewka et al. , Phys. Rev. B 78 (16), 161402R (2008). Copyright 2008 American Physical Society.

Because second-nearest neighbors are not well separated from first neighbors in amorphous structures, the short-range cutoff of the REBO-type potentials can also lead to anomalously high amounts of sp 2 -hybridized atoms in amorphous carbon (a-C). The screening function employed in the REBO-S hydrocarbon potential results in the percent sp 3 carbon as a function of density more closely resembling the results from DFT calculations. 181 Interestingly, it was possible to obtain a-C with a large percentage of sp 3 -hybridized carbon using the original parameterization of the REBO potential. 119,128 The use of a different functional form in REBO2 along with the other changes resulted in amorphous carbon with very high amounts of sp 2 carbon, much larger than DFT calculations predicted. 120  

The REBO-S potential and MD simulations have been used to study the anisotropy of wear on unsaturated diamond 188 and wear and plasticity in tetrahedral amorphous carbon (ta-C). 189 The wear of diamond was shown to proceed via an sp 3 to sp 2 transition that resulted in the formation of an amorphous layer. The growth of this layer was shown to strongly depend on the surface orientation and sliding direction, thus leading to the anisotropy observed in the polishing of diamond. MD simulations of wear in self-mated tetrahedral amorphous carbon (ta-C) films occurred via the formation of an amorphous carbon (a-C) layer with increased sp 2 content. 189 A detailed analysis of the underlying rehybridization mechanism revealed that the sp 3 to sp 2 transition is triggered by plasticity and occurs in a region that has not experienced plastic yield. This was in marked contrast to the wear of diamond which occurred in the interfacial region. 188 The main strength of the screened REBO potentials, compared to their short range antecedents, is in modeling processes where bonds are broken and formed, particularly in amorphous materials.

While the empirical bond-order potentials allow for reactions to occur between atoms of similar electronegativity, if there is a large difference in electronegativity, then a method where the charge is able to change and be equilibrated is needed if polarization and chemical reactions are to be modeled. One of the first methods that was developed to deal with charges dynamically changing is the electronegativity equalization method (EEM), 190 which uses parameters derived from quantum mechanical calculations of individual atoms. This method has been used by several groups in molecular simulations and is also referred to as charge equilibration (QEq). 191 In the EEM and QEq methods, the charge of an atom is determined iteratively, is self-consistent, and is continuously calculated throughout a simulation. The charges in the system distribute themselves so that the electronegativity at each atomic site is equalized for a given nuclear configuration. This requires the charges at each step in the simulation to reach equilibrium and traditionally is calculated by performing an inversion of an N × N matrix. As the size of the system grows by N atoms, this procedure becomes highly time-intensive, scaling as O ( N 2 ). 192 To avoid the matrix inversion, an approximate extended-Lagrangian approach to treat charges as dynamical variables that evolve in time was developed. 193 As a result, the simulation time required to calculate charge decreased significantly.

In a system where atoms or ions are interacting with one another, the total potential energy of the system comes from the sum of the self-energy terms of each atom along with the interaction energy between atoms. 192 The self-energy term is defined by the QEq method as

In the above equation, χ i is the electronegativity value of atom i , q i is the charge of atom i , and J i is the Coulombic self-repulsion term. The Yasukawa potential 194 was developed by adding charge transfer interactions to the Tersoff potential. The general form of the Yasukawa potential is

The potential energy of the atom i itself is defined by E i S , while the intermolecular interactions ( ⁠ E i j ) between atoms i and j are based on the both the distance between atoms i and j ( r ij ) and the charge on each atom ( q i and q j ). The total potential energy E i j comes from the interaction between atoms i and j that are described by four different terms: a repulsive interaction [ E i j R r i j ] ⁠ , a short-range attractive interaction ( ⁠ E i j A ) ⁠ , an ionic interaction ( ⁠ E i j I ) ⁠ , and a van der Waals interaction ( ⁠ E i j v d w )

A full description of each component of the overall potential energy equation is provided by Yasukawa and others. 192,194,195 Tests of the Yasukawa potential with applications to the Si/SiO 2 system showed that cohesive energies for various silica polymorphs produced the incorrect ordering based on energy per atom. 192  

To produce accurate energies for silica polymorphs, as well as physical properties such as lattice constants and bond lengths, the Charge-Optimized Many Body (COMB) potential was developed by Yu et al. 192 With the success of this first-generation potential, future iterations of the COMB potential were developed and have increased the transferability of the potential, with systems such as Cu, 196 Pt and Au, 197 Zr, 198 water adsorbed on copper, 199 and water during electrochemical reactions 200 being studied recently. A brief description of the most recent version of the COMB potential (COMB3) is provided below. Differences among the various versions of the COMB potential have been tabulated and discussed in a previous publication by Liang et al. 195  

The total potential energy function for the COMB potential may be written as

In the third-generation COMB potential, the electrostatic energy (E es ) term consists of self-energies as well as the charge-charge interactions between two atoms (E qq ), the interaction between the charge and the nucleus (E qZ ) and a term associated with the polarizability (E polar )

with E self given by

The self-energy is a Taylor series expansion based on the charge of an atom. In the self-energy equation, χ i and J have the same meaning as in Eq. (22) . The field term in Eq. (28) is related to the electronegativity and hardness of an atom in a particular environment. The atomic hardness of an atom may change as it moves from an atom in a small molecule to an embedded atom in a bulk system. This field-effect term allows both atomic and bulk properties to be calculated depending on an atom's environment.

In addition to the self-energy term, the electrostatic portion of the potential energy includes charge-charge interactions that are based on the electron density decay in an s -orbital, given by J i j q q

The interactions between the charge and the nucleus ( E qZ ) are dependent on the Coulombic attraction integral between electrons in the valence shell and the effective nuclear charge ( Z j )

Finally, the polarizability interaction is calculated by finding the electrostatic field generated by atomic charges

In the previous equation, μ i → is the dipole moment, E i q → is the electric field, and T ij is a damping factor that changes as the distance between atoms changes. Together, Eqs. (27)–(31) make up the electrostatic interactions in the third-generation COMB potential.

Short range interactions in the COMB formalism are similar to the short-range interactions found in the second-generation REBO potential

Here, both the repulsive and attractive energies are functions of the distance and charge on each atom in a bond. The F c (r ij ) term is a cut-off function, which is used to allow full interaction at less than a certain distance between atoms i and j , behaves like a cosine function between the minimum and maximum distance allowed between atoms i and j , and completely removes interactions if the distance is too large between atoms i and j . The cut-off function distances are parameters that are able to be changed on a system-by-system basis. The terms b ij and b ji are bond-order terms, which take into account the bonding environment around an atom.

The long-range van der Waal interactions, E i j v d w ⁠ , between atoms i and j are modeling using the Lennard-Jones 12–6 potential as in Eq. (1) . In the COMB potential, a spline function is used to smoothly turn off these long-range interactions at short interatomic distances and E i j v d w smoothly goes to zero at the cut-off distance for Coulomb interactions.

As currently implemented, the third-generation COMB potential may include additional correction terms, which are collected in the E corr r term. One type of correction utilizes Legendre polynomials to help improve bond-angle terms, in particular, systems. Another adds a bond-bending correction that takes into account the angle formed between two bonds. There are also corrections applied for repulsions at certain distances, over-coordinated atoms, and finally a charge barrier correction which helps in the electronegativity equalization. Depending on the system of study, some of these correction terms may not be required and are set to zero. However, they are retained in the functional form of the third-generation COMB potential.

With the modifications to the original Yasukawa potential, the first COMB potential was successful in predicting the structure and energies of various silica polymorphs with greater accuracy than the Yasukawa potential 192 and had good agreement with various EAM potentials 98,201 when reproducing lattice properties and polymorph energies of pure copper. 196 Second iterations of the COMB potentials have been used for Si/SiO 2 systems as well as amorphous silica 202 as well as Cu/Cu 2 O interfaces, and modeling oxidation of a Cu(100) surface. 203 The second-generation COMB potentials provided more terms in the potential energy function compared to the first COMB potential, and while it provided good agreement for systems for which it was parameterized, the results did not prove to be transferrable to systems of similar composition. For instance, while good agreement with DFT calculations existed for modeling atomic oxygen, the second iterations of the COMB potential could not model ozone sufficiently. 195 The third-generation COMB (COMB3) potential has been applied to organic molecules reacting on the Cu (111) surface, 204 titanium and titania, 205 aluminum oxide, 206 uranium and uranium oxide, 207,208 and zirconium. 209 The COMB3 potential has also been used to examine how indentation changes the distribution of charges within ZrO 2 (Fig. 8 ). 198  

Simulated charges of the (a) [110] and (b) [001] plane of ZrO 2 during a nanoindentation simulation, calculated using COMB3 potential. Oxygen atoms are more electronegative than Zr atoms, but during deformation the charges on the Zr atoms become less positive while the charges on the O atoms become less negative. Color bar indicates charges on the atoms. Reproduced with permission from Lu et al. , J. Nucl. Mater. 486 , 250–266 (2017). Copyright 2017 Elsevier B.V.

Because hydrocarbons are nonpolar molecules and carbon and hydrogen have similar electronegativities, the electronegativities of carbon and hydrogen were ignored in the parameterization of the AIREBO potential (and in the previously developed REBO2 potential). To broaden the systems that can be examined using the AIREBO potential, oxygen was added creating the qARIEBO potential. 210,211 The addition of this highly electronegative element required not only parameterizing covalent interactions for oxygen with carbon and hydrogen but also the inclusion of terms in the potential to model the partial charges of all three atoms.

The qAIREBO potential adds an electrostatic energy term, E QEq ⁠ , to an AIREBO-type potential in the total energy calculation

where E REBO has the form of the REBO2 potential. 120 The covalent bonding parameters for carbon and hydrogen remain unchanged from the REBO2 potential, while the parameters related to oxygen were adapted from the previous version of the REBO potential that contained covalent terms for oxygen but no charge or intermolecular interaction terms. 152 For optimum flexibility and accuracy, a tricubic spline 212 instead of bicubic splines 152 was adopted for the P ij corrections [Eq. (18) ] to the bond order. The additional flexibility of this choice allowed the atomization energies of a broad range C-, H-, and O-containing molecules to be well reproduced.

The terms E tors and E LJ use the AIREBO 161 formalism for treatment of the torsional and LJ contributions to the potential, respectively, but additional oxygen-related interactions needed to be fit. The term E QEq encapsulates the charge contributions to the potential energy. The simplest way to simulate different electronegativities of atoms is to have fixed partial charges on each atom at the start of the simulation. This would be problematic in an environment where chemical reactions are to be simulated. Specifically, if bonds are allowed to break and form in the course of a simulation, the associated partial charges in these bonds should also change.

As discussed above, the EEM 190 or the QEq method 191 allows the partial charges to be dynamically altered over the course of a simulation. Modeling fluctuating charges using the extended Lagrangian approach 193 with EEM seems ideally suited to examine chemical reactions using MD. However, the use of this approach requires the definition of a charge-neutral entity, such as a molecule. Using the molecule as the charge neutral entity in molecular systems can lead to atomic charges and molecular dipole moments that are too large. 213,214 Recently, the EEM method was shown to underestimate the dipole moment for large biological molecules. 215 One way to address this problem is to constrain sub-entities of large molecules to be charge neutral. Utilizing this constraint would be problematic if chemical reactions are to be simulated. In light of this, a new algorithm was developed for use with the qAIREBO potential that merged charge equilibration with bond-order potentials in a way that did not require the molecules to be charge neutral entities. This method is known as the bond-order potential/split charge equilibration (BOP/SQE). 211 This BOP/SQE method ties the amount of charge that can be transferred to the value of the bond order.

The electrostatic energy of the system, E QEq ⁠ , is composed of a self-energy term and a shifted Coulomb interaction evaluated using Slater orbitals, which adopts a spherical cut-off function. 216 The bond-centered approach to QEq discussed above is used to smoothly evolve the atomic charges in the context of the dynamic system, where bonds can form and break. In this implementation, the charge on atom i, Q i ⁠ , is defined as the sum of all the charges transferred across each of its bonds, q ¯ ij ⁠ , and is given by

where f ¯ ij is the fractional charge transferred across the bond, while q ¯ ij max limits the amount of charge transferred across the bond and is equal to q ¯ ij max =     min ( b ij , b ji ) ⁠ . By linking the fractional charge to the bond order, the charge dynamics evolve yielding a fractional-charge transfer as long as a bond exists. The fractional charge is updated solely by means of a charge acceleration, which contains terms to evolve the atomic electronegativities in time and a penalty term that limits the maximum fractional charge transferred to | f ¯ ij | < 1 ⁠ . In principle, any function that goes to zero as f ¯ ij goes to zero and diverges to infinity as | f ¯ ij | approaches one can be used for the penalty term. 211 It is the inclusion of this term that prevents charges from propagating long distances and growing too large during the QEq process.

The BOP/SQE method was used in conjunction with the qAIREBO potential for C, O, and H to simulate the confinement of sub-monolayer amounts of water (64 and 128 molecules) between two nanostructured (ultrananocrystalline diamond or UNCD) surfaces at 300 K. 217 At large surface separations, interfacial water molecules formed hydrogen-bonded droplets or aggregates that were localized to regions on both the upper and lower UNCD surfaces. Increasing the number of water molecules to 256, which corresponds to approximately one monolayer, leads to the formation of “bridges” of water molecules between the two UNCD surfaces [Figs. 9(c) and 9(d) ] during equilibration of the system. Clusters of water molecules on the UNCD surfaces and in the “bridging” structure are hydrogen bonded. In addition, the reactive nature of the qAIREBO potential allowed for water molecules to react with unsaturated sites on the UNCD surface [Fig. 9(d) ].

MD simulations using the qAIREBO potential of the equilibration of 256 water molecules at 300 K confined between two UNCD surfaces. The equilibration was carried out for a total of 25 ps. Snapshots in (a), (b), (c), and (d) correspond to 0 ps, 12.4 ps, 17.6, and 23.3 ps, respectively. In the UNCD substrates, blue and gray spheres represent carbon atoms in the grain boundary and in the diamond grain, respectively. For clarity, portions of the upper and lower UNCD surfaces are not shown. Red and yellow spheres represent oxygen and hydrogen atoms, respectively, in water molecules. The green dotted circle indicates locations where water molecules have reacted with the UNCD substrate. Simulation details are the same as those outlined in Ref. 217 .

Depending upon the amount of confined water present, compression can lead to the layering and “squeezing out” of successive layers of water. For sub-monolayer coverages of water, compression caused the formation of a single layer of water without passing through conformations with two or more layers. 217 In contrast, with enough water to form a monolayer present between the surfaces, compression leads to the formation of two distinct layers of water as shown in the density profile in Fig. 10(a) . The force on the upper UNCD surface as a function of the separation between the two surfaces is shown for the three water coverages in Fig. 10(b) . The region between the two dotted vertical lines corresponds approximately to the region where two layers of interfacial water are present. After approximately 20 nN of normal force is applied to the system with 256 molecules, two layers become one layer [Fig. 10(a) ], which is also marked by a change in slope of the force curve [Fig. 10(b) ].

(a) Density profiles of oxygen atoms in the direction perpendicular to the plane containing the UNCD surfaces at three different loads for the simulation system containing 256 water molecules. The x-axis represents the separation of the surfaces with zero corresponding to the midpoint between the two UNCD surfaces. (b) Force curves (load versus separation distance) for the compression of water molecules between UNCD surfaces. The number of water molecules is indicated in the legend. Simulation details are the same as those outlined in Ref. 217 .

The nature of MD simulations also allows for an examination of the ordering of the water molecules in contact with the UNCD surfaces as a function of separation. For sub-monolayer coverages of water on UNCD, clusters of water formed preferentially over the crystalline regions of the UNCD. 217 Snapshots of the MD simulation for the compression of a monolayer of water at various points in the simulation are shown in Fig. 11 . For loads between 0 and 10 nN, it is clear that water molecules still aggregate over the crystalline regions of the UNCD [Figs. 11(a) and 11(b) ]. At approximately 30 nN of load, there is only one distinct layer of interfacial water and these molecules are more organized on the underlying UNCD surface. It is also clear that the most of the water molecules are arranged to maximize hydrogen-bonding [Fig. 11(c) ]. It is worth noting that there still are regions on the UNCD that are not covered by water, which means that there are a few water molecules [lower left quadrant of Fig. 11(c) ] that are still “stacked up.” It takes the application of approximately 55–60 nN to cause the formation of a single monolayer of water between the UNCD surfaces, which is accompanied by another change in slope of the force curve [Fig. 10(d) ]. The water molecules in this monolayer form a well ordered, hydrogen-bonded structure [Fig. 11(d) ] that covers nearly the entire UNCD surface.

MD simulations using the qAIREBO potential of the compression of 256 water molecules at 300 K confined between two UNCD surfaces. For clarity, the upper UNCD surface is not shown and the simulation is viewed along the compression direction. The plane of the UNCD surfaces that is perpendicular to the compression direction is 40 Å × 40 Å. Snapshots in (a), (b), (c), and (d) correspond to 0 nN, 10.0 nN, 30 nN, and 90 nN, respectively. Color coding is the same as in Fig. 8 . Simulation details are the same as those outlined in Ref. 217 .

In 2001, a new reactive force field (ReaxFF) for hydrocarbons was introduced, which aimed to give a general bond-order-dependent potential in conjunction with polarizable charge and bonded and non-bonded interactions. 121 The ReaxFF potential allows for exploration of reactive chemistry at a size and scope previously unobtainable through computationally expensive quantum mechanical (QM) calculations. This potential is complicated and an extensive review of it has been published very recently. 218 The current implementation of the ReaxFF is described in detail by Russo and van Duin. 219 To date, parameterizations of the ReaxFF potential have been extended across much of the periodic table. Full functional forms can be found in the Supplementary material section of Chenoweth et al. 220 In light of the recently published review, a brief description of the ReaxFF is outlined here and the reader is directed to the published work for additional details.

The ReaxFF potential adopted the central-force formalism where all atom pairs have non-bonded interactions, which leads to smooth bond dissociation without the need for switching functions. The central-force formalism also includes added local variations that include bond stiffness, angular terms, torsion, etc., needed to describe complex molecules. The system energy for the ReaxFF potential is divided up into various contributions to the total potential energy, E t o t ⁠ , as shown in the following equation:

The first term, E bond ⁠ , assumes that the bond-order is determined directly from the interatomic distance and valence (one for hydrogen, four for carbon). The over-coordination term, E over ,   applies penalties when the coordination number of a given atom is outside the range of its ideal valence. The valence angle and torsion terms, E angle and   E t o r ⁠ , are added to calculate the energy due to variations of the valence and torsion angles from their equilibrium values. To account for non-bonded van der Waals interactions, a distance-corrected Morse-type potential term, E v d w ⁠ , is used. This term includes a shielded interaction to avoid excessively high repulsions between bonded atoms and atoms within a valence angle triplet. Coulomb interactions are taken into account for all atom pairs through the E Coulomb term and a shielded Coulomb potential. The ReaxFF potential utilizes the electronegativity equalization method 221 to distribute charges according to differences in atomic electronegativities. The E specific term represents system-specific terms that are not usually included except in special cases such as lone-pair conjugation, hydrogen bonding, and so on.

The original parameter set for the carbon and hydrogen ReaxFF potential was optimized through a successive one-parameter search using appropriate QM data. Addition of other atom types, notably Si and O 222 as well as N 223,224 led to minor changes in the functional form that enabled the potential to model lone pairs in O and a three body conjugation term for –NO 2 groups in nitramine compounds. The current form of ReaxFF, as described by Russo and van Duin, has shown considerable transferability. 219 There are currently three main branches of the ReaxFF potential. These include the original potential designed for combustion, a branch for aqueous simulations, and various other independent parameterizations. The combustion branch, originally developed to model gas-phase reactions, has been extended to include a number of oxide materials and catalysts; however, it is limited in that it was not originally designed to model liquid-phase water. The so-called aqueous branch was developed subsequently to provide a solution to this issue. The general functional form was maintained but parameters were refit to provide a reasonable description of liquid-phase water. There are also several independent branches developed to study specific systems many of which have now been merged back into the combustion branch.

One of the primary advantages of the ReaxFF potential is its ability to model reactive chemistry with an accuracy on par with QM calculations with a much lower computational cost. This has enabled the study of systems of scope and scale well beyond those accessible to QM. Despite this advantage, it is worth noting that ReaxFF is the most computationally expensive potential discussed in this review. ReaxFF has proven especially useful and advantageous over QM in modeling catalytic processes. For example, ReaxFF has been utilized in modeling Ni-catalyzed carbon nanotube growth. The growth process requires accurate modeling of bond breaking and bond formation but also diffusion of growth precursors on the Ni-catalyst surface which requires times intractable for QM dynamics simulations. 225,226 ReaxFF has also been used to study the catalytic properties of mixed metal oxides. The crystallographic structure of such materials is complex and varied with metal ions sitting in partial, mixed, and irregular positions on the crystal lattice, see Fig. 12 . The unit cell required to model such a structure is prohibitively large for a QM calculation. 220,227 Several parameter sets have been published for C/O/N/H interactions within the ReaxFF formalism. Recently, three of these parameter sets were tested by comparing the simulation of the decomposition of polyvinyl nitrate under shock and thermal loading to experimental results. 228 The authors found that ReaxFF-2016 229 predicts the onset of chemical reactions and the time scale associated with the formation of NO 2 in excellent agreement with shock experiments.

Catalytic reactions on mixed metal oxides (MMO) such as the MgAl 2 O 4 spinel illustrated here are difficult to simulate using ab initio techniques due to the size and complexity of the unit cell. Purple, yellow, and red spheres represent Mg, Al, and O, respectively. ReaxFF has been successful in modeling surface reactions in such systems. Development of new parameterizations for MMO catalysis is an active area for ReaxFF development. (See, for example, Ref. 218 .)

Recently, explicit electrons were added within the framework of the ReaxFF potential using a pseudoclassical description. This eReaxFF MD method was used to calculate electron affinities for several types of radicals and the results compared favorably to Ehrenfest dynamics. This method may prove to be a useful tool for studying large-scale systems with explicit electrons as an alternative to quantum mechanical methods. 230  

The classical potential energy functions described above all attempt to approximate the quantum-mechanical potential energy surface (PES) of a system. As seen from the above discussion, the development of accurate force fields is a laborious task and the accuracy of a potential depends on both the chosen functional form as well as the values of adjustable parameters. Recently, an increasing amount of research has focused on using machine learning (ML) methods to construct interatomic potentials. 231,232 Using flexible functional forms, which generally do not have a direct physical meaning, ML techniques such as artificial neural networks fit a PES to electronic structure calculations by adjusting a large number of parameters. 232 Such potentials can predict molecular energies with accuracy comparable to quantum mechanical methods, but at speeds up to five orders of magnitude faster. 233 Additionally, ML potentials are applicable to all types of systems and do not require the specification of bonds or atom types as in empirical potentials. Recent examples of applications of ML potentials include modeling surfaces such as SnO 2 , 234 organic molecules containing H, C, N, and O, 233 and ionic systems such as NaCl. 235 Machine learning potentials are often able to match DFT results where empirical potentials fail. For example, Si surfaces undergo complex reconstructions due to quantum mechanical effects. A recent ML potential was able to correctly predict the ordering of the stability of Si(111) reconstructions (Fig. 13 ), while the MEAM, ReaxFF, SW, and Tersoff potentials all incorrectly predict the unreconstructed surface to have the lowest energy. 236  

Catalytic SOAP-GAP predictions for silicon surfaces. (a) The tilt angle of dimers on the relaxed reconstructed Si(100) surface is the result of a Jahn-Teller distortion and is predicted to be about 19° by DFT and SOAP-GAP. Empirical force fields show no tilt. (b) The Si(111)-7 × 7 reconstruction is a complex structure which emerges from the interplay of different quantum mechanical effects. The SOAP-GAP–relaxed structure colored by predicted local energy error. (c) The energy ordering predicted on silicon surfaces by a SOAP-based ML model compared to DFT and commonly used force fields (inset). Force fields incorrectly predict the unreconstructed surface (dashed lines) to have a lower-energy state. Reproduced with permission from Bartok et al. , Sci. Adv. 3 , e1701816 (2017). Copyright 2017 Elsevier.

While ML potentials display promise, they suffer from several drawbacks. Construction of the potentials requires a sufficient training set, which usually involves performing a large number of high-level electronic structure calculations, and transferability to systems outside the training set is often limited. In addition, due to their complexity, ML potentials are typically one to two orders of magnitude slower than empirical force fields such as those discussed here. 232  

In-depth discussion of ML potentials is outside the scope of this review; further details may be found elsewhere. 231,232

The force fields described above are diverse. This diversity seems at odds with the idea that fundamentally quantum mechanics can describe any atomic system, so why such variety in these models? While the scale of systems that are treated directly by quantum mechanical methods will undoubtedly continue to grow, there will always be an accompanying interest in atomic systems much larger in size that will require the use of classical potentials to model atomic interactions. Even with a classical approach, the time scales and size scales that can be reasonably addressed are limited, in part by the complexity of the chosen interaction model. Choosing an interaction model can be largely driven by computational expense, but it is also driven by the nature of the questions being asked. A simple model, even something as simple as the LJ potential, can still sometimes be the most telling when tackling fundamental questions of a general nature. At the other extreme, very complex potentials are being used to tackle very specific questions, for instance, the ReaxFF potential provides a framework for tying DFT calculations to a parameter set that is fine-tuned for any reactive system we can build using atoms from across the periodic table (the running of DFT calculations and the optimization of parameters for particular systems can be a routine part of using ReaxFF).

Aside from the EAM potential (used primarily to describe metals which require a fundamentally different approach), most interaction potentials share a lot of common general structure that is connected to the basic conception of non-bonded and bonded interactions and in the case of reactive potentials, the description of bond breaking and forming. The specific structures found across these potentials, however, simply do not succeed all that well with regard to transferability: parameters are fitted to describe certain features and it is generally the case that these parameters do not do a great job when other properties or geometries are investigated. So, for instance, what makes the TraPPE potential distinctive is not really the chosen form for the potential, but rather its focus on phase equilibrium data to fit parameters. The problem of transferability casts another view on the issue of model complexity: more complex models are often tied more closely to the specific kinds of environments the model is trying to capture, and that complexity may be misleading if applied to systems that do not directly reflect the structure of the model systems from which parameters were fit. Simpler models, on the other hand, might be applied to a wider variety of systems though in such cases the focus will be more on general trends rather than specific details.

In short, computational expense and transferability can be viewed as the two main issues that drive the continuous development of a wide range of interaction models that serve a variety of purposes and exhibit a wide range in complexity. Most models do share a great deal of common structure, and so with a bit of experience, it is possible in many cases to choose a suitable potential to be used as is or possibly adapted to suit specific needs.

This work was supported by the Office of Naval Research (ONR) under Contract Nos. N0001416WX01648 and N0001417WX00892 and the Air Force Office of Scientific Research (AFOSR) under Contract No. F4FGA06055G002. B.H.M. also acknowledges the Research Office of the United States Naval Academy for partial support.

Citing articles via

Submit your article.

Sign up for alerts

• Online ISSN 1931-9401
• For Researchers
• For Librarians
• For Advertisers
• Our Publishing Partners  
• Physics Today
• Conference Proceedings
• Special Topics

pubs.aip.org

• Privacy Policy
• Terms of Use

Connect with AIP Publishing

This feature is available to subscribers only.

Sign In or Create an Account

Biomolecular force fields: where have we been, where are we now, where do we need to go and how do we get there?

• Perspective
• Published: 30 November 2018
• Volume 33 , pages 133–203, ( 2019 )

Cite this article

• Pnina Dauber-Osguthorpe 2 &
• A. T. Hagler   ORCID: orcid.org/0000-0002-7832-2136 1   nAff3  

2987 Accesses

74 Citations

12 Altmetric

Explore all metrics

In this perspective, we review the theory and methodology of the derivation of force fields (FFs), and their validity, for molecular simulations, from their inception in the second half of the twentieth century to the improved representations at the end of the century. We examine the representations of the physics embodied in various force fields, their accuracy and deficiencies. The early days in the 1950s and 60s saw FFs first introduced to analyze vibrational spectra. The advent of computers was soon followed by the first molecular mechanics machine calculations. From the very first papers it was recognized that the accuracy with which the FFs represented the physics was critical if meaningful calculated structural and thermodynamic properties were to be achieved. We discuss the rigorous methodology formulated by Lifson, and later Allinger to derive molecular FFs, not only obtain optimal parameters but also uncover deficiencies in the representation of the physics and improve the functional form to account for this physics. In this context, the known coupling between valence coordinates and the importance of coupling terms to describe the physics of this coupling is evaluated. Early simplified, truncated FFs introduced to allow simulations of macromolecular systems are reviewed and their subsequent improvement assessed. We examine in some depth: the basis of the reformulation of the H-bond to its current description; the early introduction of QM in FF development methodology to calculate partial charges and rotational barriers; the powerful and abundant information provided by crystal structure and energetic observables to derive and test all aspects of a FF including both nonbond and intramolecular functional forms; the combined use of QM, along with crystallography and lattice energy calculations to derive rotational barriers about ɸ and ψ; the development and results of methodologies to derive “QM FFs” by sampling the QM energy surface, either by calculating energies at hundreds of configurations, or by describing the energy surface by energies, first and second derivatives sampled over the surface; and the use of the latter to probe the validity of the representations of the physics, reveal flaws and assess improved functional forms. Research demonstrating significant effects of the flaws in the use of the improper torsion angle to represent out of plane deformations, and the standard Lorentz–Berthelot combining rules for nonbonded interactions, and the more accurate descriptions presented are also reviewed. Finally, we discuss the thorough studies involved in deriving the 2nd generation all-atom versions of the CHARMm, AMBER and OPLS FFs, and how the extensive set of observables used in these studies allowed, in the spirit of Lifson, a characterization of both the abilities, but more importantly the deficiencies in the diagonal 12-6-1 FFs used. The significant contribution made by the extensive set of observables compiled in these papers as a basis to test improved forms is noted. In the following paper, we discuss the progress in improving the FFs and representations of the physics that have been investigated in the years following the research described above.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Subscribe and save.

• Get 10 units per month
• Download Article/Chapter or eBook
• 1 Unit = 1 Article or 1 Chapter
• Cancel anytime

Price includes VAT (Russian Federation)

Instant access to the full article PDF.

Rent this article via DeepDyve

Institutional subscriptions

Similar content being viewed by others

Force field development phase ii: relaxation of physics-based criteria… or inclusion of more rigorous physics into the representation of molecular energetics.

Molecular Mechanics: Principles, History, and Current Status

Abbreviations

Arithmetic–geometric

Assisted model building with energy refinement

Atomic multipole optimized energetics for biomolecular applications

Ben Naim–Stillinger

• Consistent force field

Chemistry at HARvard Macromolecular Mechanics

Complete neglect of differential overlap

Condensed-phase optimized molecular potentials for atomistic simulation studies

Consistent valence force field

Density functional theory

Empirical conformational energy program for peptides

Extended Huckel theory

Force field

Fluctuating charges

GRO ningen MOl ecular S imulation

Hydroxyproline

Linear combination of atomic orbitals

Lennard-Jones

Least squares

Monte Carlo

Momany, Carruthers, McGuire, and Scheraga Force Field

Matsuoka–Clementi–Yoshimine

• Molecular dynamics

MDL drug data report

Molecular design limited

• Molecular mechanics

Merck molecular force field

N -methylacetamide

Optimized potential for liquid simulations

OPLS-AA/L OPLS all atom FF (L for LMP2)

Perturbative configuration interaction using localized orbitals

Protein data base

Potential energy function consortium (Biosym)

Quantum chemistry program exchange

• Quantum derivative fitting

Charge dependent polarizability

Quantum mechanics

Restrained electrostatic potential

Root mean square

Root mean square deviation

Self-consistent field-linear combination of atomic-molecular orbital (wave function)

Spectroscopically determined force fields (for macromolecules)

Simple point charge (water model)

Four point water model replacing Ben-Naim Stillinger (BNS) model

Slater-type atomic orbitals

Transferable intermolecular potential (functions for water, alcohols and ethers)

Tri- tert -butylmethane

Urey–Bradley

• van der Waals

Valence force field

Waldman–Hagler

Hofmann AW (1865) On the combining power of atoms. Proc R Inst 4:401–430

Google Scholar  

Pasteur L (1848) Recherches sur les relations qui peuvent exister entre la forme cristalline, la composition chimique, et le sens de la polarisation rotatoire. Anal Chim Phys 24:442–459

Le Bel JA (1874). Bull Soc Chim Fr 22:335–347

van’t Hoff JH (1874) Sur les formules de Structure dans l’Espace. Arch Neerl Sci Exactes Nat 9:1–10

Dreiding VAS (1959) Einfache molekularmodelle. Helv Chim Acta 48:1339–1344

Article   Google Scholar  

Koltun WL (1965) Precision space-filling atomic models. Biopolymers 3:665–679

Article   CAS   PubMed   Google Scholar  

Larson GO (1964) Atomic and molecular models made from vinyl covered wire. J Chem Educ 41:219

Article   CAS   Google Scholar  

Bernal JD, Fowler RH (1933) A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions. J Chem Phys 1:515–548

Wilson EBJ (1939) A method of obtaining the expanded secular equation for the vibration frequencies of a molecule. J Chem Phys 7:1047–1052

Barton DHR (1946) Interactions between non-bonded atoms, and the structure of cis-decalin. J Chem Soc 174:340–342

Mason EA, Kreevoy MM (1955) A simple model for barriers to internal rotation. J Am Chem Soc 77:5808–5814

Pitzer KS, Donath WE (1959) Conformations and strain energy of cyclopentane and its derivatives. J Am Chem Soc 81:3213–3218

Hendrickson James B (1961) Molecular geometry I. Machine computation of the common rings. J Am Chem Soc 83:4537–4547

Wilson EBJ, Decius JC, Cross PC (1955) Molecular vibrations: the theory of infrared and Raman vibrational spectra. McGraw-Bill, New York

Urey HC, Bradley CA (1931) The vibrations of pentatonic tetrahedral molecules. Phys Rev 38:1969

Simanouti T (1949) The normal vibrations of polyatomic molecules as treated by Urey–Bradley field. J Chem Phys 17:245–248

Dennison DM (1931) The infrared spectra of polyatomic molecules part I. Rev Mod Phys 3:280–345

Cross PC, Van Vleck JH (1933) Molecular vibrations of three particle systems with special applications to the ethyl halides and ethyl alcohol. J Chem Phys 1:350–356

Snyder RG, Schachtschneider JH (1963) Vibrational analysis of the n-paraffins-I. Assignments of infrared bands in the spectra C 3 H 8 through n-C 19 H 40 Spectrochim Acta 19:85–116

Schachtschneider JH, Snyder RG (1963) Vibrational analysis of the n-paraffins-II. Normal co-ordinate calculations. Spectrochim Acta 19:117–168

Williams JE, Stang PJ, Schleyer PVR (1968) Physical organic chemistry: quantitative conformational analysis; calculation methods. Annu Rev Phys Chem 19:531–558

Scott RA, Scheraga HA (1966) Conformational analysis of macromolecules. I. Ethane, propane, n-butane, and n-pentane. Biopolymers 4:237–238

Scott RA, Scheraga HA (1966) Conformational analysis of macromolecules. II. The rotational isomeric states of the normal hydrocarbons. J Chem Phys 44:3054–3069

Jones JE (1924) On the determination of molecular fields-II. From the equation of state of a gas. Proc R Soc London Ser A 106:463–477

Buckingham RAA (1938) The classical equation of state of gaseous helium, neon and argon. Proc R Soc London Ser A 168:264–283

Lifson S, Warshel A (1968) Consistent force field for calculations of conformations, vibrational spectra, and enthalpies of cycloalkane and n-alkane molecules. J Chem Phys 49:5116

Waldman M, Hagler ATT (1993) New combining rules for rare-gas van der Waals parameters. J Comput Chem 14:1077–1084

Hill TL (1946) On steric effects a reactive potential for hydrocarbons with intermolecular interactions on steric effects. J Chem Phys 14:465

Westheimer FH, Mayer JE (1946) The theory of the racemization of optically active derivatives of diphenyl. J Chem Phys 14:733–738

Dostrovsky BI, Hughes ED, Ingold CK (1946) Mechanism of substitution at a saturated carbon atom.I. The role of steric hindrance. J Chem Soc 173–194

Westheimer FH (1956) In: Newman MS (ed) Steric effects in organic chemistry, Chap. 12. Wiley, Hoboken

Allinger NL (1959) Conformational analysis. III Application to some medium ring compounds. J Am Chem Soc 81:5727–5733

Wiberg KB (1965) A scheme for strain energy minimization. Application to the cycloalkanes. J Am Chem Soc 87:1070–1078

Edsall JT et al (1966) A proposal of standard conventions and nomenclature for description of polypeptide conformation. J Biol Chem 241:1004–1008

CAS   PubMed   Google Scholar  

Ramachandran GN, Ramakrishnan C, Sasisekharan V (1963) Stereochemistry of polypeptide chain configurations. J Mol Biol 7:95–99

Ramakrishnan C, Ramachandran GN (1965) Stereochemical criteria for polypeptide and protein chain conformation II. Allowed conformations for a pair of peptide units. Biophys J 5:909–933

Article   CAS   PubMed   PubMed Central   Google Scholar  

Ramachandran GN, Venkatachalam CM, Krimm S (1966) Stereochemical criteria for polypeptide and protein chain conformations III. Helical and hydrogen-bonded polypepide chains. Biophys J 6:849–872

De Santis P, Giglio E, Liquori AM, Ripamonti A (1965) van der Waals interaction and the stability of helical polypeptide chains. Nature 206:456–458

Article   PubMed   Google Scholar  

Brant DA, Flory PJ (1965) The configuration of random polypeptide chains. II. Theory. J Am Chem Soc 87:2791–2800

Scott RA, Scheraga HA (1965) Method for calculating internal rotation barriers. J Chem Phys 42:2209

Scott RA, Scheraga HA (1966) Conformational analysis of macromolecules. III. Helical structures of polyglycine and poly- l -alanine. J Chem Phys 45:2091–2101

Momany FA, McGuire RF, Yan JF, Scheraga HA (1970) Energy parameters in polypeptides. 3. Semiempirical molecular orbital calculations for hydrogen-bonded model peptides. J Phys Chem 74:2424–2438

Lippincott ER, Schroeder R (1955) One-dimensional model of the hydrogen bond. J Chem Phys 23:1099

Schroeder R, Lippincott ER (1957) Potential function model of hydrogen bonds. II. J Phys Chem 61:921–928

Moulton WG, Kromhout RA (1956) Nuclear magnetic resonance: structure of the amino group. II. J Chem Phys 25:34–37

Chidambaram R, Balasubramanian R, Ramachandran GN (1970) Potential functions for hydrogen bond interactions I. A modified Lippincott–Schroeder potential function for NH–O interaction between peptide groups. Biochim Biophys Acta 221:182–195

Perutz MF et al (1960) Structure of haemoglobin—3-dimensional fourier synthesis At 5.5-Å resolution, obtained by X-ray analysis. Nature 185:416–422

Kendrew JC (1961) Three-dimensional structure of a protein molecule—way in which chain of amino acid units in a protein molecule is coiled and folded in space has been worked out for first time—protein is myoglobin, molecule of which contains 2,600 atoms. Sci Am 205:96

Lifson S (1972) In: Jaenicke R, Helmreich E (eds) Protein–proptein interactions. Springer, New York, pp 3–16

Chapter   Google Scholar  

Lifson S (1973) Recent developments in consistent force field calulations. Dyn Asp Conform Chang Biol Macromol 421–430

Bixon M, Lifson S (1967) Potential functions and conformations. Cycloalkanes Tetrahedr 23:769–784

Warshel A, Lifson S (1969) An empirical function for second neighbor interactions and its effect on vibrational modes and other properties of cyclo- and n-alkanes. Chem Phys Lett 4:255–256

Warshel A, Levitt M, Lifson S (1970) Consistent force field for calculation of vibrational spectra and conformations of some amides and lactam rings. J Mol Spectrosc 33:84–99

Warshel A, Lifson S (1970) Consistent force field calculations. II. Crystal structures, sublimation energies, molecular and lattice vibrations, molecular conformations, and enthalpies of alkanes. J Chem Phys 53:582–594

Ermer O, Lifson S (1973) Consistent force field calculations. III. Vibrations, conformations, and heats of hydrogenation of nonconjugated olefins. J Am Chem Soc 95:4121–4132

Ypma TJ (1995) Historical development of the Newton–Raphson method. SIAM Rev 37:531–551

Hirschfelder JO, Curtiss CF, Bird RB (1954) Molecular theory of gases and liquids. Wiley, Hoboken

Hagler AT, Stern PS, Lifson S, Ariel S (1979) Urey–Bradley force field, valence force field, and ab initio study of intramolecular forces in tri-tert-butylmethane and isobutane. J Am Chem Soc 101:813–819

Morse PM (1929) Diatomic molecules according to the wave mechanics. II. Vibrational levels. Phys Rev 34:57–64

Levitt M, Lifson S (1969) Refinement of protein conformations using a macromolecular energy minimization procedure. J Mol Biol 46:269–279

Momany FA, Vanderkooi G, Scheraga HA (1968) Determination of intermolecular potentials from crystal data. I. General theory and application to crystallne benzene at several temperatures. Proc Natl Acad Sci USA 61:429–436

McGuire RF et al (1971) Determination of intermolecular potentials from crystal data. II. Crystal packing with application to poly(amino acids). Macromolecules 4:112–124

Momany FA, Carruthers LM, McGuire RF, Scheraga HA (1974) Intermolecular potentials from crystal data. III. Determination of empirical potentials and applications to the packing configurations and lattice energies in crystals of hydrocarbons, carboxylic acids, amines, and amides. J Phys Chem 78:1595–1620

Momany FA, Carruthers LM, Scheraga HA (1974) Intermolecular potentials from crystal data. IV. Application of empirical potentials to the packing configurations and lattice energies in crystals of amino acids. J Phys Chem 78:1621–1630

Pople J, Beveridge D (1970) Approximate molecular orbital theory. McGraw Hill, New York

Momany FA, McGuire RF, Burgess AW, Scheraga HA (1975) Energy parameters in polypeptides.7. Geometric parameters, partial atomic charges, nonbonded interactions, hydrogen-bond interactions, and intrinsic torsional potentials for naturally occurring amino-acids. J Phys Chem 79:2361–2381

Nemethy G, Pottle MS, Scheraga HA (1983) Energy parameters in peptides. 9. Updating of geometrical parameters, nonbonded interactions, and hydrogen bond interactions for the naturally occurring amino acids. J Phys Chem 87:1883–1887

Nemethy G et al (1992) Energy parameters in polypeptldes. 10. Improved geometrical parameters and nonbonded interactions for use in the ECEPP/3 algorithm, with appllcatlon to proline-containing peptides. J Phys Chem 96:6472–6484

Allinger NL, Miller MA, Vancatledge FA, Hirsch JA (1967) Conformational analysis. LVII. The calculation of the conformational structures of hydrocarbons by the Westheimer–Hendrickson–Wiberg method. J Am Chem Soc 89:4345–4357

Allinger NL, Tribble MT, Miller MA, Wertz DH (1971) Conformational analysis. LXIX. An improved calculations of the structures and energies of hydrocarbons. J Am Chem Soc 93:1637–1648

Allinger NL (1976) Calculation of molecular structure and energy by force-field methods. Adv Phys Org Chem 13:1–82

CAS   Google Scholar  

Allinger NL (1977) Conformational analysis. 130. MM2. A hydrocarbon force field utilizing V1 and V2 torsional terms. J Am Chem Soc 99:8127–8134

Allinger NL, Hickey MJ (1975) Conformational analysis CVIII. The calculation of the structures and energies of alkanethiols and thiaalkanes by the molecular mechanics method. J Am Chem Soc 97:5167–5177

Allinger NL, Hickey MJ, Kao J (1976) Conformational analysis CXI. The calculation of the structures and energies of disulfides by the molecular mechanics method. J Am Chem Soc 98:2741–2745

Allinger NL, Kao J (1976) Conformational analysis XCIV. Molecular mechanics studies of sulfoxides. Tetrahedron 32:529–536

Allinger NL, Chung DY (1976) Conformational analysis. 118. Application of the molecular-mechanics method to alcohols and ethers. J Am Chem Soc 98:6798–6803

Allinger NL, Chang SHM (1977) Conformational analysis—CXXIII. Tetrahedron 33:1561–1567

Williams DE (1966) Nonbonded potential parameters derived from crystalline aromatic hydrocarbons. J Chem Phys 45:3770–3778

Williams DE (1967) Nonbonded potential parameters derived from crystalline hydrocarbons. J Chem Phys 47:4680–4684

Kitaigorodskii AI (1965) The principle of close packing and the condition of thermodynamic stability of organic crystals. Acta Crystallogr 18:585–590

Kitaigorodskii AI (1968) Studies on organic chemical crystallography in USSR. Sov Phys Crystallogr USSR 12:692

Kitaigorodskii AI, Mirskaya KV, Tovbias AB (1968) Lattice energy of crystalline benzene in atom-atom approximation. Sov Phys Crystallogr USSR 13:176

Kitaigorodskii AI, Mirskaya KV (1969) Atom-atom potential method in physics of molecular crystal. Acta Cryst A 25:S91

Kitaigorodskii AI, Mirskaya KV (1965) Quadrupole interaction in a molecular crystal. Sov Phys Crystallogr 10:121

Hagler AT, Huler E, Lifson S (1974) Energy functions for peptides and proteins. I. Derivation of a consistent force field including the hydrogen bond from amide crystals. JAm Chem Soc 96:5319–5327

Hagler AT, Lifson S (1974) Energy functions for peptides and proteins. II. The amide hydrogen bond and calculation of amide crystal properties. J Am Chem Soc 96:5327–5335

Hagler A, Lifson SA (1974) Procedure for obtaining energy parameters from crystal packing. Acta Cryst B30:1336–1341

Hagler ATT, Dauber P, Lifson S (1979) Consistent force field studies of intermolecular forces in hydrogen-bonded crystals. 3. The C=O⋯H–O hydrogen bond and the analysis of the energetics and packing of carboxylic acids. J Am Chem Soc 101:5131–5141

Murthy ASN, Rao CNR (1970) Recent theoretical treatments of the hydrogen bond. J Mol Struc 6:253–282

Pearlman Da et al (1995) AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules. Comput Phys Commun 91:1–41

Stockmayer WH (1941) Second virial coefficients of polar gases second virial coefficients of polar gases. J Chem Phys 9:398–402

Liquori AM (1960) The stereochemical code and the logic of a protein molecule. Q Rev Biophys 2:65–92

Berkovitch-yellin Z, Leiserowitz L (1980) The role of Coulomb forces in the crystal packing of amides. A study based on experimental electron densities. J Am Chem Soc 102:7677–7690

Spackman MA (1987) A Simple quantitative model of hydrogen bonding. Application to more complex systems. J Phys Chem 91:3179–3186

Spackman MA, Weber HP, Craven BM (1988) Energies of molecular interactions from Bragg diffraction data. J Am Chem Soc 110:775–782

Dinur U, Hagler AT (1992) The role of nonbond and charge flux in hydrogen bond interactions. The effect on structural changes and spectral shifts in water dimer. J Chem Phys 97:9161–9172

Shaik MS, Liem SY, Popelier PLA (2010) Properties of liquid water from a systematic refinement of a high-rank multipolar electrostatic potential. J Chem Phys 132:174504

Ren P, Wu C, Ponder JW (2011) Polarizable atomic multipole-based molecular mechanics for organic molecules. J Chem Theory Comput 7:3143–3161

Bakó I et al (2010) Hydrogen bonded network properties in liquid formamide. J Chem Phys 132:14506

Lifson S, Hagler AT, Dauber P (1979) Consistent force field studies of intermolecular forces in hydrogen-bonded crystals. 1. Carboxylic acids, amides, and the C=O⋯H-hydrogen bonds. J Am Chem Soc 101:5111

Hagler AT, Lifson S, Dauber P (1979) Consistent force field studies of intermolecular forces in hydrogen-bonded crystals. 2. A benchmark for the objective comparison of alternative force fields. J Am Chem Soc 101:5122–5130

Morozov AV, Kortemme T (2006) Potential functions for hydrogen bonds in protein structure prediction and design. Adv Protein Chem 72:1–38

McGuire RF, Momany FA, Scheraga HA (1972) Energy parameters in polypeptides. V. An empirical hydrogen bond potential function based on molecular orbital calculations. J Phys Chem 76:375–393

Sippl MJ, Ncmethy G, Scheraga HA (1984) Intermolecular potentials from crystal data. 6. Determination of empirical potentials for O–H⋯O=C hydrogen bonds from packing conflguratlons. J Phys Chem 7:6231–6233

Arnautova YA, Jagielska A, Scheraga HA (2006) A new force field (ECEPP-05) for peptides, proteins, and organic molecules. J Phys Chem B 110:5025–5044

Weiner SJ et al (1984) A new force-field for molecular mechanical simulation of nucleic-acids and proteins. J Am Chem Soc 106:765–784

Cornell WD, Kollman, PA (1995) A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J Am Chem Soc 117:5179–5197

Mccammon JA, Wolynes PG, Karplus M (1979) Picosecond dynamics of tyrosine side chains in proteins. Biochemistry 18:927–942

Brooks BR et al (1983) CHARMM: a program for macromolecular energy, minimization, and dynamics calculations. J Comput Chem 4:187–217

Neria E, Fischer S, Karplus M (1996) Simulation of activation free energies in molecular systems. J Chem Phys 105:1902–1921

Allinger NL, Yuh YH, Jenn-Huei L (1989) Molecular mechanics. The MM3 force field for hydrocarbons. J Am Chem Soc 11:8551–8566

Allinger NL, Chen K, Lii J (1996) An improved force field (MM4) for saturated hydrocarbons. J Comp Chem 17:642–668

Hagler AT (2018) The second phase of FF research introduces a fork in the road: relaxation of physics-based criteria through the introduction of correction factors… or Inclusion of more rigorous physics into the representation of molecular energetics. J Comput Aided Mol Des. https://doi.org/10.1007/s10822-018-0134-x

Pullman B, Pullman A (1959) The electronic structure of the purine-pyrimidine pairs of DNA. Biochim Biophys Acta 36:343–350

Pullman B (1969) Opening remarks: quantum-mechanical calculations of biological structures and mechanisms. Ann N Y Acad Sci 158:1–19

Hückel VE (1931) Quantentheoretische beitrgge zum benzolproblem. I. Die elekfronenkonfigurafion des benzols und verwandfer verbindungen. Zeitsch Phys 70:204–286

DelRe G, Pullman B, Yonezawa T (1963) Electronic structure of the alpha-amino acids of proteins. I. Charge distributions and proton chemical shifts. Biochim Biophys Acta 75:153–182

DelRe G (1958) Charge distributions in saturated organic molecules. A simple MO-LCAO method. J Chem Soc 4031–4040

Hoffmann R (1963) An extended hückel theory. I. Hydrocarbons. J Chem Phys 39:1397

Pople JA, Santry DP, Segal GA (1965) Approximate self-consistent molecular orbital theory. I. Invariant procedures. J Chem Phys 43:S129–S129

Hoffmann R, Imamurs A (1969) Quantum mechanical approach to the conformational analysis of macromolecules in ground and excited states. Biopolymers 7:207–213

Rossi A, David CW, Schor R (1969) Extended Hückel calculations on polypeptide chains. The alpha helix. Theor Chim Acta 14:429–431

Kier LB, George JM (1969) Extended Hiickel MO calculations of the conformation of several amino acids. Theor Chim Acta 260:258–260

Govil G (1970) Molecular orbital calculations on polypeptides and proteins. Part 1. Extended Huckel theory calculations on trans-dipeptides. J Chem Soc A 2464–2469

Govil G (1971) Molecular orbital calculations on polypeptides and proteins. Part 2 stability and conformational structure of cis-dipeptides. J Chem Soc A 386–388

Govil G (1971) Molecular orbital calculations on polypeptides and proteins. 3 conformation of side chains. J Indian Chem Soc 48:731

Govil G, Saran A (1971) Molecular orbital calculations on polypeptides and proteins. 4. Studies by EHT and CNDO on the influence of hydrogen bond and chain length on certain structures. J Chem Soc A 3624–3627

Govil G, Saran A (1972) Molecular orbital calculations on polypeptides and proteins 5. Stable conformations of oxygen-containing peptides. J Chem Soc Faraday Trans 68:1176–1180

Yan JF, Momany FA, Hoffmann R, Scheraga HA, McGuire RF (1970) Energy parameters in polypeptides 2. Semiempirical molecular orbital calculations for model peptides. J Phys Chem 74:420

Yan JF, Momany FA, Scheraga HA (1970) Conformational analysis of macromolecules. VI. Helical structures of 0-, rn-, and p-chlorobenzyl esters of poly- l -aspartic acid. J Am Chem Soc 92:1109–1115

Momany FA, McGuire RF, Yan JF, Scheraga H (1971) Energy parameters in polypeptides. IV. Semiempirical molecular orbital calculations of conformational dependence of energy and partial charge in di- and tripeptides. J Phys Chem 75:2286–2297

Diner S, Malrieu JP, Claverie P (1969) Localized bond orbitals and the correlation problem I. The perturbation calculation of the ground state energy. Theor Chim Acta 13:1–17

Maigret B, Pullman B, Dreyfus M (1970) Molecular orbital calculations on the conformation of polypeptides and proteins. 1. Preliminary investigations and simple dipeptides. J Theor Biol 26:321–333

Maigret M, Pullman B, Caillet J (1970) The conformational energy of an alanyl residue preceding proline: a QM approach. Biochem Biophys Res Comm 40:808–813

Pullman B, Pullman A (1974) Molecular orbital calculations on the conformation of amino acid residues of proteins. Adv Prot Chem 28:347–526

Boyd DB (2013) Pioneers of quantum chemistry: ACS symposium series, vol. 1122. American Chemical Society, Washington, DC, pp 221–273

Book   Google Scholar  

Clementi E (1962) SCF-MO wave functions for the hydrogen fluoride molecule. J Chem Phys 36:33

Veillard A, Clementi E (1967) Complete multi. configuration self-consistent field theory. Theor Chim Acta 7:133–143

Clementi E (1972) Computation of large molecules with the Hartree–Fock model. Proc Natl Acad Sci USA 69:2942–2944

Clementi E, Davis DR (1967) Electronic structure of large molecular systems. J Comp Phys 2:223–244

Clementi E, Andre JM, Andre MC, Klint D, Hahn D (1969) Study of the electronic structure of molecules. X ground state for adenine, cytosine, guanine and thyamine. Acta Phys Acad Sci Hungaricae 27:439–521

Clementi E, Mehl J, von Niessen W (1971) Study of the electronic structure of molecules. XII. Hydrogen bridges in the guanine–cytosine pair and in the dimeric form of formic acid. J Chem Phys 54:508

Popkie H, Kistenmacher H, Clementi E (1973) Study of the structure of molecular complexes. IV. The Hartree–Fock potential for the water dimer and its application to the liquid state. J Chem Phys 59:1325

Kistenmacher H (1974) Study of the structure of molecular complexes. VI. Dimers and small clusters of water molecules in the Hartree–Fock approximation. J Chem Phys 61:546

Matsuoka O, Clementi E, Yoshimine M (1976) CI study of the water dimer potential surface. J Chem Phys 64:1351

Barker JA, Watts RO (1969) Structure of water: a Monte Carlo calculation. Chem Phys Lett 3:144–145

Narten AH (1972) Liquid water: atom pair correlation functions from neutron and X-ray diffraction. J Chem Phys 56:5681

Lie GC, Clementi E (1975) Study of the structure of molecular complexes. XII. Structure of liquid water obtained by Monte Carlo simulation with the Hartree–Fock potential corrected by inclusion of dispersion forces. J Chem Phys 62:2195

Bartlett RJ, Shavitt I, Purvis GD (1979) The quartic force field of H 2 O determined by many-body methods that include quadruple excitation effects. J Chem Phys 71:281

Lie GC, Clementi E (1986) Molecular-dynamics simulation of liquid water with an ab initio fiexible water-water interaction potential. Phys Rev A 33:2679–2693

Corongiu G, Clementi E (1993) Molecular dynamics simulations with a flexible and polarizable potential: density of states for liquid water at different temperatures. J Chem Phys 98:4984–4990

Corongiu G, Clementi E (1992) Liquid water with an ab initio potential: X-ray and neutron scattering from 238 to 368 K. J Chem Phys 97:2030

Corongiu G, Clementi E (1993) Solvated water molecules and hydrogen-bridged networks in liquid water. J Chem Phys 98:2241–2249

Sciortino F, Corongiu G (1993) Structure and dynamics in hexagonal ice: a molecular dynamics simulation with an ab initio polarizable and flexible potential. J Chem Phys 98:5694

Clementi E, Cavallone F, Scordamaglia R (1977) Analytical potentials from ‘ab Initio’ computations for the Interaction between biomolecules. 1. Water with amino acids. J Am Chem Soc 99:5531–5545

Scordamaglia R, Cavallone F, Clementi E (1977) Analytical potentials from ‘ab initio’ computations for the interaction between biomolecules. 2. Water with the four bases of DNA. J Am Chem Soc 99:5545–5550

Clementi E, Corongiu G, Lelj F (1979) Analytical potentials from ab initio computations for the interaction between biomolecules. V. The phosphate group in nucleic acids. J Chem Phys 70:3726

Clementi E, Corongiu G (1979) Interaction of water with DNA single and double helix in the B conformation. Int J Quantum 16:897–915

Aida M, Corongiu G, Clementi E (1992) Ab initio force field for simulations of proteins and nucleic acids. Int J Quantum Chem 42:1353–1381

Corongiu G (1992) Molecular dynamics simulation for liquid water using a polarizable and flexible potential. Int J Quantum Chem 42:1209–1235

Clementi E, Corongiu G (1985) Computer simulations of complex chemical systems. Adv Biophys 20:75–107

Clementi E et al (1991) Selected topics in ab initio computational chemistry in both very small and very large chemical systems. Chem Rev 91:699–879

Pople JA (1970) Molecular orbital methods in organic chemistry. Acct Chem Res 3:217–223

Hehre WJ, Radom L, Schleyer PVR, Pople JA (1986) Ab initio molecular orbital theory. Wiley, Hoboken

Pople J (1999) Nobel lecture: quantum chemical models. Rev Mod Phys 71:1267–1274

Hehre WJ, Stewart RF, Pople JA (1969) Self-consistent molecular-orbital methods. I. Use of Gaussian expansions of slater-type atomic orbitals. J Chem Phys 51:2657

Hehre WJ, Lathan WA, Ditchfield R, Newton MD, Pople JA (1970) Gaussian 70. Quantum Chem Exch Progr 237

Del Bene J, Pople JA (1969) Intermolecular energies of small water polymers. Chem Phys Lett 4:426–428

Del Bene J, Pople JA (1970) Theory of molecular interactions. I. Molecular orbital studies of water polymers using a minimal slater-type basis. J Chem Phys 52:4858

Del Bene JE (1973) Theory of molecular interactions. III. A comparison of studies of H 2 O polymers using different molecular-orbital basis sets. J Chem Phys 58:3605

Hagler AT, Leiserowitz L, Tuval M (1976) Experimental and theoretical studies of barrier to rotation about N–C-alpha and C-alpha-C′ bonds (phi and psi) in amides and peptides. J Am Chem Soc 98:4600–4612

Bernstein J, Hagler AT (1978) Conformational polymorphism the influence of crystal structure on molecular conformation. J Am Chem Soc 100:673–681

Hagler AT, Bernstein J (1978) Conformational polymorphism 2. Crystal energetics by computational substitution—further evidence for sensitivity of method. J Am Chem Soc 100:673–681

Bernstein J, Hagler AT (1979) Polymorphism and isomorphism as tools to study the relationship between crystal forces and molecular-conformation. Mol Cryst Liq Cryst 50:223–233

Lommerse JPM et al (2000) A test of crystal structure prediction of small organic molecules research papers. Acta Cryst B 56:697–714

Van Eijck BP (2002) Crystal structure predictions using five space groups with two independent molecules. The case of small organic acids. J Comput Chem 23:456–462

Dunitz JD, Gavezzotti A (2009) How molecules stick together in organic crystals: weak intermolecular interactions. Chem Soc Rev 38:2622–2633

Bardwell D et al (2011) Towards crystal structure prediction of complex organic compounds—a report on the fifth blind test. Acta Crystallogr B 67:535–551

Hornak V et al (2006) Comparison of multiple amber force fields and development of improved protein backbone parameters. Proteins-Struct Funct Bioinforms 65:712–725

Feig M, MacKerell AD, Brooks CL (2003) Force field influence on the observation of π-helical protein structures in molecular dynamics simulations. J Phys Chem B 107:2831–2836

Bochevarov AD et al (2013) Jaguar: a high-performance quantum chemistry software program with strengths in life and materials sciences. Int J Quantum Chem 113:2110–2142

Kitano M, Fukuyama T, Kuchitsu K (1973) Molecular structure of N -methylacetamide as studied by gas electron diffraction. Bull Chem Soc Jpn 46:384–387

Kitano M, Kuchitsu K (1974) Molecular structure of N -methylformamide as studied by gas electron diffraction. J Bull Chem Soc Jpn 47:631–634

Kitano M, Kuchitsu K (1973) Molecular structure of acetamide as studied by gas electron diffraction. J Bull Chem Soc Jpn 46:3048–3051

Brock CP, Ibers JA (1975) The role of crystal packing forces in the structure of pentaphenylantimony coordinate molecules. In such molecules the trigonal. Acta Cryst A 31:38–42

Hagler AT, Leiserowitz L (1978) Amide hydrogen-bond and anomalous packing of adipamide. J Am Chem Soc 100:5879–5887

Guo H, Karplus M (1992) Ab initio studies of hydrogen bonding of N -methyiacetamide: structure, cooperativlty, and internal rotational barriers. J Phys Chem 96:7273–7287

Mennucci B, Martínez JM (2005) How to model solvation of peptides? Insights from a quantum-mechanical and molecular dynamics study of N -methylacetamide. 1. Geometries, infrared, and ultraviolet spectra in water. J Phys Chem B 109:9818–9829

Villani V, Alagona G, Ghio C (1999) Ab initio studies on N -methylacetamide. Molec Eng B 8:135–153

Katz L, Post B (1960) The crystal structure and polymorphism of N -methylacetamide. Acta Crystallogr 13:624–628

Rauscher S et al (2015) Structural ensembles of intrinsically disordered proteins depend strongly on force field: a comparison to experiment. J Chem Theory Comput 11:5513–5524

Hagler AT (1977) Relation between spatial electron-density and conformational properties of molecular systems. Isr J Chem 16:202–212

Hagler AT, Lapiccirella A (1976) Spatial electron distribution and population analysis of amides, carboxylic acid, and peptides and their relation to empirical poteintial functions. Biopolymers 15:1167–1200

Hagler AT, Lapiccirella A (1978) Basis set dependence of spatial electron-distribution—implications for calculated conformational equilibria. J Am Chem Soc 100:4026–4029

Nyburg SC, Faerman CH (1985) A revision of van der Waals atomic radii for molecular crystals: N, O, F, S, CI, Se, Br and I bonded to carbon. Acta Cryst B 41:274–279

Stone AJ, Price SL (1988) Some new ideas in the theory of intermolecular forces: anisotropic atom-atom. J Phys Chem 92:3325–3335

Price SL, Leslie M, Welch GWA, Habgood M, Price LS, Karamertzanis PG, Day GM (2010) Modelling organic crystal structures using distributed multipole and polarizability-based model intermolecular potentials. Phys Chem Chem Phys 12:8478–8490

Eramian H, Tian Y-H, Fox Z, Beneberu HZ, Kertesz M, Se S (2013) On the anisotropy of van der Waals atomic radii of O. F, Cl, Br, and I. J Phys Chem A 117:14184–14190

Hagler AT (2015) Quantum derivative fitting and biomolecular force fields—functional form, coupling terms, charge flux, nonbond anharmonicity and individual dihedral potentials. J Chem Theory Comput 11:5555–5572

Hagler AT (1977) On the relation between the spatial electron density and the conformational properties of molecular systems. Isr J Chem 16:202–212

Rowlinson JS (1951) The lattice energy of ice and the second virial coefficient of water vapour. Trans Faraday Soc 47:120–129

Berendsen HJC, Postma JPM, van Gunsteren WF, Hermans J (1981) In: Pullman B (ed) Intermolecular forces. Reidel Publishing Company, Dordrecht, pp 331–342

Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML (1983) Comparison of simple potential functions for simulating liquid water. J Chem Phys 79:926–935

Berkovitch-Yellin Z, Leiserowitz L (1982) Atom–atom potential analysis of the packing characteristics of carboxylic acids. A study based on experimental electron density distributions. J Am Chem Soc 104:4052–4064

Spackman MA (2012) Charge densities and crystal engineering in modern charge-density analysis. Springer, Dordrecht

Mannfors B, Palmo K, Krimm S (2000) A new electrostatic model for molecular mechanics force fields. J Mol Struct 556:1–21

Palmo K, Mannfors B, Mirkin NG, Krimm S (2003) Potential energy functions: from consistent force fields to spectroscopically determined polarizable force fields. Biopolymers 68:383–394

Qian W, Krimm S (2005) Limitations of the molecular multipole expansion treatment of electrostatic interactions for C–H⋯O and O–H⋯O hydrogen bonds and application of a general charge density approach. J Phys Chem A 109:5608–5618

Shi Y et al (2013) Polarizable atomic multipole-based AMOEBA force field for proteins. J Chem Theory Comput 9:4046–4063

Weiner PK, Kollman PA, AMBER (1981) Assisted model building with energy refinement. A general program for modeling molecules and their interactions. J Comput Chem 2:287–303

McCammon JA, Gelin BR, Karplus M (1977) Dynamics of folded proteins. Nature 267:585–590

Engler EM, Andose JD, Schleyer PVR (1973) Critical evaluation of molecular mechanics. J Am Chem Soc 95:8005

Weiner SJ, Kollman PA, Nguyen DT, Case DA (1986) An all atom force field for simulations of proteins and nucleic acids. J Comput Chem 7:230–252

van Gunsteren WF, Berendsen HJC (1987) GROningen MOlecular Simulation (GROMOS) Library Manual. Biomos, Groningen

Hermans J, Berendsen HJC, van Gunsteren WF, Postma JPM (1984) A consistent empirical potential for water–protein interactions. Biopolymers 23:1513–1518

Jorgensen WL (1981) Transferable intermolecular potential functions for water, alcohols, and ethers. Application to liquid water. J Am Chem Soc 103:335–340

Jorgensen WL (1982) Revised TIPS for simulations of liquid water and aqueous solutions. J Chem Phys 77:4156

Jorgensen WL, Swenson CJ (1985) Optimized intermolecular potential functions for amides and peptides. Structure and properties of liquid amides. J Am Chem Soc 107:569–578

Jorgensen WL, Tirado-Rives J (1988) The OPLS potential functions for proteins. Energy minimizations for crystals of cyclic peptides and crambin. J Am Chem Soc 110:1657–1666

Dauber P, Osguthorpe DJ, Hagler AT, Structure (1982) Energetics and dynamics of ligand binding to dihydrofolate-reductase. Biochem 10:312–318

Dauber-Osguthorpe P et al (1988) Structure and energetics of ligand binding to proteins: Escherichia coli dihydrofolate reductase-trimethoprim, a drug-receptor system. Proteins Struct Funct Bioinform 4:31–47

Hagler AT, Osguthorpe DJ, Dauber-Osguthorpe P, Hempel JC (1985) Dynamics and conformational energetics of a peptide hormone: vasopressin. Science 227:1309–1315

Struthers RS, Rivier J, Hagler AT (1985) Molecular dynamics and minimum energy conformations of GnRH and analogs. Ann N Y Acad Sci 439:81–96

Struthers RS et al (1990) Design of biologically-active, conformationally constrained Gnrh antagonists. Proteins-Struct Funct Genet 8:295–304

Dauber-Osguthorpe P et al (1988) Structure and energetics of ligand binding to proteins: Escherichia coli dihydrofolate reductase-trimethoprim, a drug-receptor system. Proteins-Struct Funct Genet 4:31–47

Rick SW, Stuart SJ (2002) Potentials and algorithms for incorporating polarizability in computer simulations. Rev Comput Chem 18:89–146

Vega C, Abascal JLF (2011) Simulating water with rigid non-polarizable models: a general perspective. Phys Chem Chem Phys 13:19663–19688

Wang L-P et al (2013) Systematic improvement of a classical molecular model of water. J Phys Chem B 117:9956–9972

Cardamone S, Hughes TJ, Popelier PLA (2014) Multipolar electrostatics. PhysChemChemPhys 16:10367–10387

Clark GNI, Cappa CD, Smith JD, Saykally RJ, Head-Gordon T (2010) The structure of ambient water. Mol Phys 108:1415–1433

Demerdash O, Yap E-H, Head-Gordon T (2014) Advanced potential energy surfaces for condensed phase simulation. Annu Rev Phys Chem 65:149–174

Ben-Naim A, Stillinger FH (1972) Aspects of the statistical-mechanical theory of water. Wiley, Hoboken

Rahman A, Stillinger FH (1971) Molecular dynamics study of liquid water. J Chem Phys 55:3336–3359

Alder BJ, Wainwright TE (1959) Studies in molecular dynamics. I. General method. J Chem Phys 31:459

Alder BJ, Wainwright TE (1960) Studies in molecular dynamics. II. Behavior of a small number of elastic spheres. J Chem Phys 33:1439

Stillinger FH, Rahman A (1974) Improved simulation of liquid water by molecular dynamics. J Chem Phys 60:1545–1557

Mahoney MW, Jorgensen WL (2000) A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions. J Chem Phys 112:8910

Rick SW, Stuart SJ, Berne BJ (1994) Dynamical fluctuating charge force fields: application to liquid water. J Chem Phys 101:6141–6156

Chen B, Xing J, Siepmann JI (2000) Development of polarizable water force fields for phase equilibrium calculations. J Phys Chem B 104:2391–2401

Horn HW et al (2004) Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew. J Chem Phys 120:9665–9678

Abascal JLF, Vega C (2005) A general purpose model for the condensed phases of water: TIP4P/2005. J Chem Phys 123:234505

Bauer BA, Warren GL, Patel S (2009) Incorporating phase-dependent polarizability in nonadditive electrostatic models for molecular dynamics simulations of the aqueous liquid—vapor interface. J Chem Theory Comput 359–373

Bauer BA, Patel S (2009) Properties of water along the liquid-vapor coexistence curve via molecular dynamics simulations using the polarizable TIP4P-QDP-LJ water model. J Chem Phys 131:84709

Cisneros GA et al (2016) Modeling molecular interactions in water: from pairwise to many-body potential energy functions. Chem Rev 116:7501 – 7528

Lipkowitz KB, Allinger NL, Lipkowitz KB, Allinger NL (1987) QCPE Bull 7

Allinger NL, Lii J-H (1987) Benzene, aromatic rings, van der Waals molecules, and crystals of aromatic molecules in molecular mechanics (MM3). J Comput Chem 8:1146–1153

Lii JH, Allinger NL (1989) Molecular mechanics. The MM3 force field for hydrocarbons. 2. Vibrational frequencies and thermodynamics. J Am Chem Soc 111:8566–8575

Lii JH, Allinger NL (1989) Molecular mechanics. The MM3 force field for hydrocarbons. 3. The van der Waals’ potentials and crystal data for aliphatic and aromatic hydrocarbons. J Am Chem Soc 111:8576–8582

Ermer O (1976) Calculation of molecular properties using force fields. Applications in organic chemistry. Bond Forces Struct Bond 27:161–211

Lifson S, Stern PS (1982) Born–Oppenheimer energy surfaces of similar molecules: interrelations between bond lengths, bond angles, and frequencies of normal vibrations in alkanes. J Chem Phys 77:4542

Allinger NL, Rahman M, Lii JA (1990) Molecular mechanics force field (MM3) for alcohols and ethers. J Am Chem Soc 112:8293–8307

Schmitz LR, Allinger NL (1990) Molecular mechanics calculations (MM3) on aliphatic amines. J Am Chem Soc 112:8307–8315

Allinger NL, Quinn M, Rahman M, Chen K (1991) Molecular mechanics (MM3) calculations on sulfides. J Phys Org Chem 4:647–658

Chen K, Allinger NL (1991) Molecular mechanics (MM3) calculations on disulfides. J Phys Org Chem 4:659–666

Allinger NL, Zhu ZQS, Chen K (1992) Molecular mechanics (MM3) studies of carboxylic acids and esters. J Am Chem Soc 114:6120–6133

Tai JC, Yang L, Allinger NL (1993) Molecular mechanics (MM3). Calculations on nitrogen-containing aromatic heterocycles. J Am Chem Soc 115:11906–11917

Lii J, Allinger NL (1991) The MM3 force field for amides, polypeptides and proteins. J Comput Chem 12:186–199

Lii J, Allinger NL (1994) Directional hydrogen bonding in the MM3 force field. I. J Phys Org Chem 7:591–609

Dauber P, Hagler AT (1980) Crystal packing hydrogen bonding, and the effect of crystal forces on molecular conformation. Acct Chem Res 13:105–112

Nevins N, Lii JH, Allinger NL (1996) Molecular mechanics (MM4) calculations on conjugated hydrocarbons. J Comput Chem 17:695–729

Nevins N, Allinger NL (1996) Molecular mechanics (MM4) vibrational frequency calculations for alkenes and conjugated hydrocarbons. J Comput Chem 17:730–746

Langley CH, Allinger NL (2002) Molecular mechanics (MM4) calculations on amides. J Phys Chem A 106:5638–5652

Chen K, Lii J, Fan Y, Allinger NL (2007) Molecular mechanics (MM4) study of amines. J Comput Chem 28:2391–2412

Allinger NL, Chen K, Lii J, Durkin KA (2003) Alcohols, ethers, carbohydrates, and related compounds. I. The MM4 force field for simple compounds. J Comput Chem 24:1447–1472

Nevins N, Chen K, Allinger NL (1996) Molecular mechanics (MM4) calculations on alkenes. J Comput Chem 17:669–694

Langley CH, Allinger NL (2003) Molecular mechanics (MM4) and ab initio study of amide-amide and amide-water dimers. J Chem Phys A 107:5208–5216

Cornell WD et al (1995) A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J Am Chem Soc 117:5179–5197

Wang J, Cieplak P, Kollman PA (2000) How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J Comput Chem 21:1049–1074

Brüesch P (1966) X-ray and infrared studies of bicyclo(2.2.2) octane, triethylenediamine and quinuclidine—II. Normal co-ordinate calculations of bicyclo(2.2.2.)octane, triethylenediamine and quinuclidine. Spectrochim Acta 22:867–875

Bayly CI, Cieplak P, Cornell WD, Kollman PA (1993) A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J Phys Chem 97:10269–10280

MacKerell AD et al (1998) All-atom empirical potential for molecular modeling and dynamics studies of proteins. J Phys Chem B 102:3586–3616

Gelin BR, Karplus M (1975) Sidechain torsional potentials and motion of amino acids in proteins: bovine pancreatic trypsin inhibitor. Proc Natl Acad Sci USA 72:2002–2006

Maple JR, Hwang MJ, Stockfisch TP, Hagler AT (1994) Derivation of class II force fields.3. Characterization of a quantum force field for alkanes. Isr J Chem 34:195–231

Debiec KT et al (2016) Further along the road less traveled: AMBER ff15ipq, an original protein force field built on a self-consistent physical model. J Chem Theory Comput 12:3926 – 3947

Huang J et al (2017) Charmm36M: an improved force field for folded and intrinsically disordered proteins. Nat Methods 14:71–73

Jorgensen WL, Maxwell DS, Tirado-Rives J (1996) Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J Am Chem Soc 118:11225–11236

Kaminski GA, Friesner RA, Tirado-rives J, Jorgensen WL (2001) Evaluation and reparametrization of the OPLS-AA force field for proteins via comparison with accurate quantum chemical calculations on peptides. J Phys Chem B 105:6474–6487

Beachy MD, Chasman D, Murphy RB, Halgren TA, Friesner RA (1997) Accurate ab initio quantum chemical determination of the relative energetics of peptide conformations and assessment of empirical force fields. J Am Chem Soc 119:5908–5920

Daura X, Mark AE, van Gunsteren WF (1998) Parametrization of aliphatic CH United Atoms of GROMOS96 force field. J Comp Chem 19:535–547

Schuler LD, Daura X, van Gunsteren WF (2001) An improved GROMOS96 force field for aliphatic hydrocarbons in the condensed phase. J Comp Chem 22:1205–1218

Halgren TA (1996) Merck molecular force field.I. Basis, form, scope, parameterization, and performance of MMFF94. J Comput Chem 17:490–519

Schuler LD, van Gunsteren WF (2000) On the choice of dihedral angle potential energy functions for n-alkanes. Mol Simul 25:301–319

Oostenbrink C, Villa A, Mark AE, van Gunsteren WF (2004) A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6. J Comput Chem 25:1656–1676

Schmid N et al (2011) Definition and testing of the GROMOS force-field versions 54A7 and 54B7. Eur Biophys J 40:843–856

Horta BAC, Fuchs PFJ, van Gunsteren WF, Hunenberger Philippe H (2011) New interaction parameters for oxygen compounds in the GROMOS force field: improved pure-liquid and solvation properties for alcohols, ethers, aldehydes, ketones, carboxylic acids and esters. J Chem Theory Comput 7:1016–1031

Horta BAC et al (2012) Reoptimized interaction parameters for the peptide-backbone model compound N -methylacetamide in the GROMOS force field: influence on the folding properties of two beta-peptides in methanol. J Comput Chem 33:1907–1917

Arnautova Ya, Jagielska A, Pillardy J, Scheraga H (2003) Derivation of a new force field for crystal-structure prediction using global optimization: nonbonded potential parameters for hydrocarbons and alcohols. J Phys Chem B 107:7143–7154

Jagielska A, Arnautova YA, Scheraga HA (2004) Derivation of a new force field for crystal-structure prediction using global optimization: nonbonded potential parameters for amines, imidazoles, amides, and carboxylic acids. J Phys Chem B 108:12181–12196

Pillardy J, Arnautova Y, Czaplewski C, Gibson KD, Scheraga HA (2001) Conformation-family Monte Carlo: a new method for crystal structure prediction. Proc Natl Acad Sci USA 98:12351–12356

Williams DE, Cox SR (1984) Nonbonded potentials for azahydrocarbons: the importance of the coulombic interaction. Acta Cryst B40:404–417

Berkovitch-Yellin Z (1985) Toward an ab initio derivation of crystal morphology. J Am Chem Soc 107:8239–8253

Mooij WTM, van Eijck BP, Price SL, Verwer P, Kroon J (1998) Crystal structure predictions for acetic acid. J Comput Chem 19:459–474

Hehre WJ, Pople JA (1968) Atomic electron populations for some simple molecules. Chem Phys Lett 2:379–380

Cox SR, Williams DE (1981) Representation of the molecular electrostatic potential by a net atomic charge model. J Comput Chem 2:304–323

Scrocco E, Tomasi T (1978) Electronic molecular structure, reactivity and intermolecular forces: an euristic interpretation by means of electrostatic molecular potentials. Adv Quant Chem I1:116–193

Besler BH, Merz KM, Kollman PA (1990) Atomic charges derived from semiempirical methods. J Comput Chem 11:431–439

Clementi E (1985) Ab initio computational chemistry. J Phys Chem 89:4426–4436

Zhou T, Huang D, Caflisch A (2010) Quantum mechanical methods for drug design. Curr Top Med Chem 10:33–45

Lopes PEM, Guvench O, MacKerell AD (2015) Molecular modeling of proteins. In: Kukol A (ed) Methods in molecular biology, vol. 1215. Springer, New York, pp 47–71

Moore GE, Fellow L (1965) Cramming more components onto integrated circuits. Electronics 114–117

Becke AD, Perspective (2014) Fifty years of density-functional theory in chemical physics. J Chem Phys 140:18A301

St-Amant A, Salahub DR (1990) New algorithm for the optimization of geometries in local density functional theory. Chem Phys Lett 169:387–392

Bajorath J, Kraut J, Li ZQ, Kitson DH, Hagler AT (1991) Theoretical-studies on the dihydrofolate-reductase mechanism—electronic polarization of bound substrates. Proc Natl Acad Sci USA 88:6423–6426

Andzelm JW, Nguyen DT, Eggenberger R, Salahub DR, Hagler AT (1995) Applications of the adiabatic connection method to conformational equilibria and reactions involving formic-acid. Comput Chem 19:145–154

Halgren TA (1996) Merck molecular force field. II. MMFF94 van der Waals and electrostatic parameters for intermolecular interactions. J Comput Chem 17:520–552

Halgren TA (1996) Merck molecular force field. III. Molecular geometries and vibrational frequencies for MMFF94. J Comput Chem 17:553–586

Halgren TA, Nachbar RB (1996) Merck molecular force field. IV.Conformational energies and geometries for MMFF94. J Comput Chem 17:587–615

Halgren TA (1996) Merck molecular force field. V. Extension of MMFF94 using experimental data, additional computational data, and empirical rules. J Comput Chem 17:616–641

Maple JR, Dinur U, Hagler AT (1988) Derivation of force fields for molecular mechanics and dynamics from ab initio energy surfaces. Proc Natl Acad Sci USA 85:5350–5354

Maple JR et al (1994) Derivation of class-Ii force-fields.1. Methodology and quantum force-field for the alkyl functional-group and alkane molecules. J Comput Chem 15:162–182

Hagler AT, Ewig CS (1994) On the use of quantum energy surfaces in the derivation of molecular-force fields. Comput Phys Commun 84:131–155

Dinur U, Hagler AT (1991) Lipkowitz KB, Boyd DB (eds) Reviews in computational chemistry, vol. 15. VCH, New York, pp 99–164

Probe Manual (1989) Biosym Technologies, Inc

Ewig CS et al (2001) Derivation of class II force fields. VIII. Derivation of a general quantum mechanical force field for organic compounds. J Comput Chem 22:1782–1800

Palca J (1986) Computer models: cooperation on new molecules. Nature 322:586

Halgren TA (1992) The representation of van der Waals (vdW) interactions in molecular mechanics force fields: potential form, combination rules, and vdW parameters. J Am Chem Soc 114:7827–7843

Ewig CS, Thacher TS, Hagler AT (1999) Derivation of class II force fields. 7. Nonbonded force field parameters for organic compounds. J Phys Chem B 103:6998–7014

Halgren TA, MMFF VII (1999) Characterization of MMFF94, MMFF94s, and other widely available force fields for conformational energies and for intermolecular interaction energies and geometries. J Comput Chem 20:730–748

Bordner AJ, Cavasotto CN, Abagyan RA (2003) Direct derivation of van der Waals force field parameters from quantum mechanical interaction energies. J Phys Chem B 107:9601–9609

Hagler AT, Maple JR, Thacher TS, Fitzgerald GB, Dinur U Weiner PK (1989) Potential energy functions for organic and biomolecular systems. Comput Simul Biomol Syst ESCOM Leiden 1:149–167

Hwang MJ, Stockfisch TP, Hagler AT (1994) Derivation of class II force fields. 2. Derivation and characterization of a class-II force field, CFF93, for the alkyl functional group and alkane molecules. J Am Chem Soc 116:2515–2525

Ercolessi F, Adams JB (1994) Interatomic potentials from first-principles calculations: the force-matching method. Eur Lett 26:583–588

Dasgupta S, Goddard WA (1989) Hessian-biased force fields from combining theory and experiment. J Chem Phys 90:7207–7215

Maple JR, Hwang MJ, Jalkanen KJ, Stockfisch TP, Hagler AT (1998) Derivation of class II force fields: V. Quantum force field for amides, peptides, and related compounds. J Comput Chem 19:430–458

Palmo K, Mirkin NG, Pietila: L, Krimm S (1993) Spectroscopically determined force fields for macromolecules. 1 n-alkane chains. Macromolecules 26:6831–6840

Stouch TR (2012) The errors of our ways: taking account of error in computer-aided drug design to build confidence intervals for our next 25 years. J Comput Aided Mol Des 26:125–134

Ermer O, Lex J (1987) Shortened C–C, bondsand antiplanar O=C–0–H torsion angles in 1,4-cubanedicarboxyiic acid. Anyen Chem Int Ed Engl 41:447–449

Wilson E, Decius J, Cross PC (1980) Molecular vibrations. Dover, New York

Hwang MJ, Ni X, Waldman M, Ewig CS, Hagler AT (1998) Derivation of class II force fields. VI. Carbohydrate compounds and anomeric effects. Biopolymers 45:435–468

Chen K-H, Allinger NL (2002) Molecular mechanics (MM4) study of saturated four-membered ring hydrocarbons. J Mol Struct Theochem 581:215–237

Irwin JJ, Sterling T, Mysinger MM, Bolstad ES, Coleman RG (2012) ZINC: a free tool to discover chemistry for biology. J Chem Inf Model 52:1757–1768

MDDR Accelrys. http://accelrys.com/products/collaborative-science/databases/bioactivity-databases/mddr.html

Sun H, Mumby SJ, Maple JR, Hagler AT (1994) An ab initio Cff93 all-atom force field for polycarbonates. J Am Chem Soc 116:2978–2987

Sun H, Rigby D (1997) Polysiloxanes: ab initio force field and structural, conformational and thermophysical properties title. Spectrochim Acta A 53:1301–1323

Hill JR, Sauer J (1994) Molecular mechanics potential for silica and zeolite catalysts based on ab initio calculations. 1. Dense and microporous silica. J Phys Chem 98:1238–1244

Zhu W et al (2009) Molecular dynamics simulations of AP/HMX composite with a modified force field. J Hazard Mater 167:810–816

McQuaid MJ, Sun H, Rigby D (2004) Development and validation of COMPASS force field parameters for molecules with aliphatic azide chains. J Comput Chem 25:61–71

Sun HCOMPASS (1998) An ab initio force-field optimized for condensed-phase applications—overview with details on alkane and benzene compounds. J Phys Chem B 102:7338–7364

Dinur U, Hagler AT (1989) Direct evaluation of nonbonding interactions from ab initio calculations. J Am Chem Soc 111:5149–5151

Dinur U, Hagler AT (1989) Determination of atomic point charges and point dipoles from the Cartesian derivatives of the molecular dipole moment and second moments, and from energy second derivatives of planar dimers II. Applications to model systems. J Chem Phys 91:2959–2970

Dinur U, Hagler AT (1995) Geometry-dependent atomic charges—methodology and application to alkanes, aldehydes, ketones, and amides. J Comput Chem 16:154–170

Galimberti D, Milani A, Castiglioni C (2013) Charge mobility in molecules: charge fluxes from second derivatives of the molecular dipole. J Chem Phys 138:164115

Cruz-Cabeza AJ, Bernstein J (2014) Conformational polymorphism. Chem Rev 114:2170–2191

Dauber P, Hagler AT, Crystal, Packing (1980) Hydrogen-bonding, and the effect of crystal forces on molecular-conformation. Acc Chem Res 13:105–112

Liang CX, Ewig CS, Stouch TR, Hagler AT (1993) Abinitio studies of lipid model species.1. Dimethyl-phosphate and methyl propyl phosphate anions. J Am Chem Soc 115:1537–1545

Liang CX, Ewig CS, Stouch TR, Hagler AT (1994) Ab initio studies of lipid model species. 2. Conformational-analysis of inositols. J Am Chem Soc 116:3904–3911

Liang C et al (1995) Force field studies of cholesterol and cholesteryl acetate crystals and cholesterol–cholesterol intermolecular interactions. J Comp Chem 16:883–897

Lorentz HA (1881) Ueber die anwendung des satzes vom virial in der kinetischen theorie der gase. Ann Phys 248:127–136

Berthelot D (1898) Sur le mélange des gaz. Comptes Rendus Hebd Séances l’Acad Sci 126:1703–1855

Diaz Peña M, Pando C, Renuncio JA (1982) R. Combination rules for intermolecular potential parameters. I. Rules based on approximations for the long-range dispersion energy. J Chem Phys 76:325

Tang KT, Toennies JP (1986) New combining rules for well parameters and shapes of the van der Waals potential of mixed rare gas systems. Z Phys D 1:91–101

Kestin J et al (1984) Equilibrium and transport properties of the noble gases and their mixtures at low density. J Phys Chem Ref Data 13:229–303

Al-Matar AK, Rockstraw DW (2004) A generating equation for mixing rules and two new mixing rules for interatomic potential energy parameters. J Comp Chem 25:660–668

Desgranges C, Delhommelle J (2014) Evaluation of the grand-canonical partition function using expanded Wang-Landau simulations. III. Impact of combining rules on mixtures properties. J Chem Phys 140:104109

Song W, Rossky PJ, Maroncelli M (2003) Modeling alkane + perfluoroalkane interactions using all-atom potentials: failure of the usual combining rules. J Chem Phys 119:9145

Reinhold J et al (2014) Molecular dynamics simulations on scattering of single Ar, N 2 , and CO 2 molecules on realistic surfaces. Comput Fluids 97:31–39

Nemethy G, Scheraga HA (1965) Theoretical determination of sterically allowed conformations of a polypeptide chain by a computer method. Biopolymers 3:155–184

Allinger NL (2011) Understanding molecular structure from molecular mechanics. J Comput Aided Mol Des 25:295–316

Clementi E, Ranghino G, Scordamaglia R (1977) Intermolecular pontentials: interaction of water with lysozyme. Chem Phys Lett 49:218–224

Dinur U, Hagler AT (1990) A novel decomposition of torsional potentials into pairwise interactions—a study of energy 2nd derivatives. J Comput Chem 11:1234–1246

Dinur U, Hagler AT (1989) Determination of atomic point charges and point dipoles from the cartesian derivatives of the molecular dipole-moment and 2nd moments, and from energy 2nd derivatives of planar dimers. I. Theory. J Chem Phys 91:2949–2958

Ooi T, Scott RA, Vanderkooi G, Scheraga HA (1967) Conformational analysis of macromolecules. IV. Helical structures of poly- l -alanine, poly- l -valine, poly-β-methyl- l -aspartate, poly-γ-methyl- l -glutamate, and poly- l -tyrosine. J Chem Phys 46:4410–4426

Download references

Acknowledgements

We would like to thank Dr. Mike Gilson, for reading parts of the manuscript and helpful discussions and Dr. Ruth Sharon for reading, discussing and valuable help with editing. We also thank Eitan Hagler for help with the figures. Special thanks to the editor, Dr. Terry Stouch for his invitation to write this perspective, encouragement, and endless patience.

Author information

A. T. Hagler

Present address: Valis Pharma, San Diego, CA, USA

Authors and Affiliations

Department of Chemistry, University of Massachusetts, Amherst, MA, 01003, USA

Art Pearl Studios, Newlyn, Cornwall, UK

Pnina Dauber-Osguthorpe

You can also search for this author in PubMed   Google Scholar

Corresponding author

Correspondence to A. T. Hagler .

Rights and permissions

Reprints and permissions

About this article

Dauber-Osguthorpe, P., Hagler, A.T. Biomolecular force fields: where have we been, where are we now, where do we need to go and how do we get there?. J Comput Aided Mol Des 33 , 133–203 (2019). https://doi.org/10.1007/s10822-018-0111-4

Download citation

Received : 17 November 2017

Accepted : 09 March 2018

Published : 30 November 2018

Issue Date : 15 February 2019

DOI : https://doi.org/10.1007/s10822-018-0111-4

Share this article

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

• Force fields: force field derivation
• Potential functions
• Hydrogen bond: drug discovery
• Protein simulation
• Molecular simulation
• Nonbond interactions
• Combination rules
• Polarizability
• Charge flux
• Nonbond flux
• Polarizability flux
• Free energy
• Coupling terms
• Cross terms
• Electrostatics
• Multipole moments
• Find a journal
• Publish with us
• Track your research

An official website of the United States government

Official websites use .gov A .gov website belongs to an official government organization in the United States.

Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites.

• Publications
• Account settings
• Advanced Search
• Journal List

Recent Developments and Applications of the CHARMM force fields

Pedro em lopes, alexander d mackerell jr.

• Author information
• Article notes
• Copyright and License information

Correspondence: [email protected]

Issue date 2012 January/February.

Empirical force fields commonly used to describe the condensed phase properties of complex systems such as biological macromolecules are continuously being updated. Improvements in quantum mechanical (QM) methods used to generate target data, availability of new experimental target data, incorporation of new classes of compounds and new theoretical developments (eg. polarizable methods) make force-field development a dynamic domain of research. Accordingly, a number of improvements and extensions of the CHARMM force fields have occurred over the years. The objective of the present review is to provide an up-to-date overview of the CHARMM force fields. A limited presentation on the historical aspects of force fields will be given, including underlying methodologies and principles, along with a brief description of the strategies used for parameter development. This is followed by information on the CHARMM additive and polarizable force fields, including examples of recent applications of those force fields.

Introduction

Empirical force fields (FF) have become integral tools in the study of complex phenomena in numerous fields such as materials science and biomolecular simulations. The first protein molecular dynamics (MD) simulation was performed in 1970 on a system of 458 atoms over a length of 9.2 ps on an IBM 370, a top-of-the-line computer at the time, using precursor versions of the program CHARMM and CHARMM force field.( 1 ) From its initial conception in the late 1970s to the present time, CHARMM and the associated force field have evolved to allow for MD simulation studies to be performed on a host of systems of different sizes, functions, and complexity. The current version of the CHARMM program, which has been extensively reviewed by Brooks et al ( 2 ), is compatible with and has been optimized for multiple types of computer platforms capable of running either in serial or in parallel ranging from single processors PCs to clusters of multi-core processor machines to shared-memory supercomputer installations.

The CHARMM additive force field, one of the most successful force fields in the study of biomolecules, was first described in the seminal paper by Karplus and co-workers in 1983( 3 ). Despite being an integral part of the CHARMM program, it can also be used with a number of MD programs such as AMBER ( 4 ), NAMD ( 5 ) and GROMACS ( 6 ). It has been the force field of choice in many landmark MD simulations, having been used, for example, by Shulten and co-workers in 2006 for the first MD simulation that broke the one million atom MD simulation barrier. Their work on the satellite tobacco mosaic virus, which marked the first all-atom computational model of a full life form, was published two years later.( 7 ) The CHARMM force field was also used in MD simulations performed with ANTON, a specialized supercomputer that executes MD simulations orders of magnitude faster than was previously possible.( 8 ) These advances put additional demands on the force field with respect to accuracy requiring continuous evolution of the parameters comprising the force fields. Indeed, the ongoing evolution of the additive CHARMM force field is what allows it to continue to be of utility for the study of molecular systems of more complex nature and under a longer simulation time frame.

The purpose of the current review is to highlight recent developments in the CHARMM force fields and with emphasis on work originating from our laboratory. Topics covered include strategies employed in parameter optimization, limitations of the force field, and their potential application in the discovery of new therapeutic agents. While the earliest additive CHARMM force field focused on proteins and nucleic acids it has been since extended to lipids ( 9 , 10 ), carbohydrates ( 11 – 15 ), and to small, typically drug-like, molecules ( 16 ) as well as more specialized chemicals, such as quartz ( 17 ). In addition, recent enhancements in the protein ( 18 , 19 ), lipid ( 9 ) and nucleic acid ( 20 ) force fields have been presented. Current progress in the treatment of induced polarization in the context of the Drude-based polarizable CHARMM force field will also be discussed.

Molecular Mechanics and the CHARMM potential energy function

The molecular mechanics approach was originally pioneered by Alder and Wainwright in the late 1950’s( 21 ). It consists of using a “balls-on-springs” scheme to represent molecules where each ball corresponds to an atom and each spring a covalent bond. This representation allows molecules to be modeled through Newtonian mechanics. Calculation of the forces and integration of the Newton’s equation of motion allows the time evolution of a given system to be calculated, thus giving birth to the field of MD simulations. In molecular mechanics (MM) the potential energy of a system of particles is determined using a defined functional form and associated parameter set, the combination of which is called a force field. In MM force fields the properties of a system of particles are defined by the elasticity of the covalent bonds, the suppleness of two or three adjacent bonds (valence and dihedral angles), as well as the interactions that arise between non-bonded atoms (van der Waals “vdW” and electrostatic interactions). These elements represent the individual contributions to the energy of the system, which are ultimately combined allowing for the energy and forces of a system to be calculated as a function of geometry.

The potential energy function comprising the CHARMM additive force field is

where R⃗ is the ensemble of coordinates of atoms in a system. Internal or intramolecular parameters include bonds ( b ), angles ( θ ), Urey-Bradley ( UB, S ), dihedrals ( χ ), and impropers (□) where K corresponds to their associated force constants and the symbols having subscript “0” corresponding to the equilibrium values. Non-bonded terms include Lennard-Jones (LJ) parameters ε ij (well-depth) and R min (distance of minimum interaction energy), which are used to define the vdW interactions, and partial charges q , used in Coloumb’s law to describe electrostatic contributions. Force fields based on this type of energy function are referred to as “additive” as they don’t allow the electrostatic parameters to change as a function of environment, such that the total electrostatic energy is simply a sum of all the individual atom-atom terms. Though a large number of force fields have been developed based on the additive model and successfully applied to the study of biomolecular systems, this representation does not account for the phenomenon of electronic polarization and may be considered a significant approximation in these types of force fields.

Explicit inclusion of electronic polarization is expected to represent an important improvement in empirical force fields. For example, in proteins, it was shown that non-polarizable force fields share a systematic error in the calculated electrostatic potential due to a lack of treatment of electronic polarizability.( 22 ) To compensate for the lack of polarizability the additive CHARMM FF, as well as most other macromolecular force fields, deliberately overestimates molecular dipoles of compounds in the gas phase in order to be able to reproduce electrostatic interactions that occur in condensed phase environments. In other words, the parameters are optimized to better reproduce condensed-phase properties at the expense of agreement with gas-phase data. As a basic illustration of the importance of polarization, a water molecule has a gas-phase dipole moment of 1.85 D ( 23 ) but an average dipole moment of 2.1 D in a dimer and even higher value in large water clusters ( 24 ). The value for the dipole moment in condensed or liquid phase is indicated to be as high as 2.9 D in X-ray and neutron diffraction experiments ( 25 ) and ab initio MD simulations have predicted values around 3.0 D ( 26 ). Although systematic overestimation of molecular dipoles has worked well in calculating structure and energetics of molecules in the condensed phase, it fundamentally lacks the ability to describe locally induced polarization events, an essential element in intermolecular interactions.

Several of the major energy functions for biological macromolecules, including CHARMM, have been extended to include polarization effects with varying degrees of success. The efforts have been centered around using three methods including induced dipoles, fluctuating point charges, and classical Drude oscillators to describe polarizability.( 27 ) Details on the various implementations have been previously discussed in other reviews, ( 27 , 28 ) an issue of the Journal of Chemical Theory and Computation that specially addressed polarization ( 29 ) and in a more recently published review ( 30 ). The classical Drude oscillator is currently being implemented in the context of CHARMM by the MacKerell and Roux groups.( 31 ).

The implementation of molecular polarizability in CHARMM, based on the classical Drude oscillator formalism, involves a massless point charge (ie. the Drude oscillator) being attached to each non-hydrogen atom via a harmonic spring ( 32 )( Figure 1 ). To model polarization, the position of the additional point charge reacts to the surrounding electrostatic environment, thereby effectively shifting the electron density of the system leading to a change in the atomic and molecular dipoles. To include this effect the electrostatic part of the potential energy function ( Eq. 1 ) is extended by the following equation

Schematic representation of how atomic polarizability is treated in the CHARMM polarizable force field using methanol as an example. Drude oscillators (or particles)(blue, “D”) are attached to non-hydrogen parent atoms through harmonic bonds (dashed lines). Oxygen lone-pairs (green, “LP”) are connected with constrained bonds, angles, and dihedrals relative to the C-O-H plane. Hydrogens are not considered as polarizable entities in this model.

for a molecule of N non-hydrogen atoms, where q D and q C are the charges on the Drude particle and its parent atom, respectively; r D and r C are their respective positions, and k D is the force constant of the harmonic spring between the Drude oscillator and its parent atom. In this model, the atomic polarizability (α) of an atom A is related to the charge of the attached Drude particle and the associated harmonic spring through the following equation:

Electrostatic interactions between the Drudes and other charged centers for the 1,4 atom pairs (i.e. atoms separated by 3 covalent bonds) and beyond are treated through Coulomb’s Law the same way as in the additive force field ( 33 ). However, 1,2 and 1,3 electrostatic interactions (i.e. atoms separated by 1 or 2 covalent bonds, respectively), which are ignored in the CHARMM additive force field, are reintroduced for Drude particles. These interactions, corresponding to the 1,2 q D (A) - q C (A) , or “Drude-parent”, and the 1,3 q D (A)-q C (B) , or Drude-adjacent parent, interactions described in Equation 2 , are required to create the correct dipole response in a molecule. However, since these electrostatic interactions occur at close proximities, the point charge approximation is unsuitable. To circumvent this issue, the CHARMM Drude force field introduces a damping function initially developed by Thole ( 34 ). The Thole screening function

is applied to interactions between two atom centers with charges q i and q j , where r ij is the distance between the atom centers and the normalized distance is defined as

such that α i and α j are the polarizabilities of respective charge centers and α, the Thole parameter.

In addition to the Drude model two other polarizable models have been introduced into CHARMM. The polarizable fluctuating charge model has been implemented by Patel, Brooks and co-workers.( 35 , 36 ) This method allows the values of the partial charges to respond to the electric field of their environment effectively altering the molecular dipole using an electronegativity or chemical potential equalization scheme. The CHARMM implementation of the fluctuating charge model has been applied to study small-molecule liquid-vapor interfaces, solvated proteins/peptides, and physiological membrane systems, in which calculation of free energies of solvation was made possible.( 37 ) In addition, Gao and coworkers have introduced the induced dipole model into CHARMM ( 38 ). In this model a tensor is applied to each atom in the system to define polarization. Implementation of multiple methods to treat electronic polarization in CHARMM is a feature that will facilitate future force field development.

Parameter Optimization: Strategies and Considerations

Independent of the form of the potential energy function, a force field is only as good as the quality of the associated parameters. To obtain parameters of sufficient accuracy a procedure is followed in which the parameters in Equation 1 are optimized to reproduce a variety of target properties such as molecular geometries and vibrations, pure solvent properties, and free energies of solvation, among others. The general outline of the parametrization process has been described for the CHARMM additive force field in earlier publications ( 33 , 39 ) and is depicted in Figure 2 . Note that parameter optimization is an iterative approach and several rounds of parametrization are typically required until convergence is achieved. A common strategy in parameter optimization of biological macromolecules is decomposition of the large molecules into smaller model compounds. The advantages are twofold. First, it breaks complex systems into smaller models that are easier to treat using both MM and QM methods. Second, for smaller systems more experimental data are available for smaller systems, including thermodynamic properties of condensed phases, such as heats of vaporization or free energies of aqueous solvation, much of which is essential for accurate optimization of the nonbond portion of the force field. This approach assumes that parameters from the smaller models are transferable to the macromolecules, an assumption that is valid to a limited extent. Figure 3 illustrates how proteins can be decomposed into smaller model compounds; N-methyl acetamide describes the peptide bond, the alanine dipeptide models the protein backbone, and each sidechain is modeled by a corresponding valence-saturated model compound. Similar decompositions are employed for nucleic acids, carbohydrate, lipids, as well as in the construction of larger drug-like molecules from fragments. A well-thought optimization strategy is not only the key to obtain accurate reproduction of QM and experimental target data but also to reduce the labor required to parametrize a molecule.

A simplified parametrization flow chart.

Illustration of model compounds used in the parametrization of the CHARMM additive protein force field. (a) N-methyl acetamide is used to model the peptide bond; (b) sidechains, such as in PHE, are modeled by analogous compounds that include terminating methyl or ethyl groups; (c) alanine dipeptide is the model compound for optimization of the ϕ/ψ torsional parameters including CMAP corrections.

Selection of Target Data

A very important step in parameter optimization is the selection of target data. Force field parameter determination is an under-determined and nontrivial problem difficult to automate and very much dependent on human intervention. For the CHARMM force fields a combination of QM data and experimentally determined properties are considered. Although increases in computer power have allowed bigger systems to be studied with more accurate methods, sole reliance on QM methods is inappropriate. This is particularly important for optimization of the nonbond parameters, where dispersion interactions play an important role. However, it should be emphasized that the availability of high quality experimentally determined structural properties (e.g. geometries, vibrations, and conformations) is limited to small molecules, including constituents of nucleic acids and proteins, leading to a reliance on QM data. Macroscopic condensed phase properties such as heats of vaporization or sublimation and free energies of hydration are available for many of the building blocks from which the force field is hierarchically constructed. Although many of these data do not provide an atomic level picture of molecular interactions, when combined with the atomic details provided by QM and structure-based experimental methods (ex. X-ray crystallography and neutron diffraction) they are ideally able to produce a balanced force field, i.e. one that performs well at the atomic level but is also able to reproduce condensed phase experimental properties. The difficulty of attaining this was noted in a 1995 review paper by Halgren ( 40 ) in which he commented on the future of force fields:

“The past two years have seen a veritable beehive of activity in force field development. Remarkably, this period has seen the (…) completion of significant re-parameterizations of two currently widely applied force fields - CHARMM and AMBER. That these force fields differ substantially in form and in manner of derivation serves to emphasize that force field development is still as much a matter of art as of science. Someday, consensus on the form and manner of parameterization of molecular force fields may exist, but for now much remains to be learned. This is as it should be, for the problem being addressed is a hard one: to capture faithfully in a computationally tractable model enough of the real, quantum-mechanical physics to insure that a bio-molecular simulation, properly carried out, will yield a correct answer to useful precision.”

Fifteen years later the vision of Halgren has yet to materialize. Force field development is still as much of an art as it is a science. Improvements in QM methods alone have yet to solve the problem and human intervention continues to be extremely important to properly weight the various QM and experimental data being considered during parameters optimization. Finally, problems associated with parameter correlation continue to be a problem.

Optimization procedure

Parametrization of the CHARMM force field consists of obtaining appropriate intramolecular (bond, angle, dihedral, Urey-Bradley and improper, or “out-of-plane”, terms), vdW, and electrostatic parameters that adequately reproduce the selected target data. Determination of the electrostatic parameters differs between the additive and the Drude polarizable force fields in that the polarizabilities ( Eq. 3 ) and Thole factors ( Eq. 4 ) need to be optimized in addition to the point charges in the additive force field. Here, we describe the strategies behind the various steps depicted in Figure 2 that are involved in parameter optimization for the additive and Drude polarizable force fields. As was suggested earlier, parametrization is a tedious process that directly affects the outcome of subsequent usage. It is therefore a requirement that this task be performed in a careful and consistent fashion.

The development of an accurate electrostatic model is one of the more tedious steps in the parametrization workflow, laying the foundation for the model to reproduce the proper electronic response as a function of environment. In the additive CHARMM force field optimization of point charges is based on a supramolecular approach where the charges are adjusted to reproduce QM HF/6-31G(d) interaction energies and geometries of the model compound with, typically, individual water molecules as illustrated in Figure 4 . Placement of water molecules at different orientations around the molecule has the distinct advantage of implicitly accounting for local electronic polarization, a desirable feature for accurate reproduction of condensed phase properties. In addition, interactions with water are often supplemented with dimers of the model compounds and with dipole moments of the models, with the dipole moments overestimated to account for condensed phase effects.( 39 , 41 ). For a number of biomolecular force fields the most common optimization methods are based on QM electrostatic potential (ESP) approaches. ESP methods consist of determining atomic charges based on reproduction of QM ESP evaluated on a grid surrounding the model compound.( 42 – 44 ) The advantage of these methods is of convenience since charges can be determined quickly for any compound for which the QM ESP can be determined. An extension of ESP methods is the inclusion of restraints during fitting, referred to as the restrained ESP (RESP) approach,( 45 ) which overcomes limitations on the determination of charges on buried atoms.( 46 ) It is important to take note that solutions from both supramolecular and ESP approaches are conformation dependent, requiring care in selecting the appropriate conformation or conformations when performing the charge optimization.

Orientation of test water molecules around the polarizable atom of interest exemplified through ethanethiol. Carbon is in cyan, sulfur in yellow and LP in blue. “120” is probing for interactions along the S-LP axis. “180” is oriented towards the C–S bond. “BIS” water points at the bisection of the C-S-H valence angle. “PR” is pointed directly at the sulfur proton and along the S–H bond.

Electrostatic parameters in the Drude polarizable force field are obtained from restrained fitting to perturbed QM ESP maps calculated on grid points located on concentric Connolly surfaces around the molecule. In some cases fitting is supplemented with reproduction of the molecular dipole moment and diagonal elements of the polarization tensor.( 31 ) The determination of the atomic polarizabilities and Thole factors ( 34 ) requires multiple perturbed ESPs, typically calculated at the B3LYP/aug-cc-pVDZ level, each giving the electronic response of the molecule to a point charge. Perturbing ions are placed mainly along chemical bonds and lone pairs. To more accurately reproduce hydrogen-bond interactions as a function of orientation lone pairs (LPs) are added to hydrogen bond acceptors. LPs typically carry the charge of the atom to which they are attached (N, O or S) while the associated polarizability and Thole factor are assigned to the parent atom. In addition, anisotropic polarizability of hydrogen bond acceptors was found to be required to reproduce interactions with ions as a function of orientation (see reference ( 47 ) for an example on sulfurs and Figure 5 for an illustration). Although gas-phase dipoles are readily reproducible with full atomic polarizabilities, scaling of the polarizabilities was shown to be necessary to reproduce condensed-phase properties as seen in the SWM4-DP and SWM4-NDP water models.( 48 ) A scaling factor of approximately 0.7 was found appropriate for the water models but for other classes of molecules scaling factors ranging from 0.6 to 1.0 were shown to be appropriate with 1.0 being full polarizability. For instance, a scaling factor of 0.7 was used in thiols and a factor of 0.85 was needed for dimethyl disulfide and 0.6 for ethylmethyl sulfide ( 47 ). Other scaling factors are 0.7 for primary and secondary alcohols ( 49 ), 0.85 for aromatics ( 50 ), N-containing heterocycles ( 51 ), nucleic acid bases and ethers ( 52 ) and 1.0 for alkanes ( 53 ). A value of 0.724, identical to the water model, was recently used in ion parameters( 54 ). Final optimization of the polarizabilities consists of testing the model for reproduction of the pure solvent dielectric and adjusting the polarizability scaling if necessary. A generalization can be made here such that, although scaling factors will differ from one class of model compounds to another, they will remain in the aforementioned range.

The anisotropic nature of polarizability for a sulfur atom in ethylmethylsulfide is probed using +0.5e point charges placed on two perpendicular arcs (inset). Differences between QM perturbed and unperturbed electrostatic potentials are used to determine polarization response as a function of orientation. Polarizability anisotropy parameters are fitted to reproduce relative response curves obtained through quantum mechanical calculations.

Determination of the vdW terms is one of the most important aspects in the development of a force field for condensed phase simulations and is indeed the most tedious step in the process. Jorgensen pioneered the use of condensed phase simulations, usually pure liquids, as the basis for optimization of the LJ parameters.( 55 , 56 ) Typically, once partial atomic charges are assigned via the supramolecular approach, the LJ parameters for a model compound can be adjusted to reproduce the experimentally determined pure solvent heat of vaporization and density, as well as isothermal compressibility and heat capacity, crystal heat of sublimation and lattice geometry and/or free energy of aqueous solvation, as available. Although this is an effective method for the fine-tuning of the parameters, the solution is by far trivial. A major issue is the parameter correlation paradigm, such that LJ parameters for different atoms in a molecule (e.g., H and C in ethane) as well as between the ε ij and R min terms on an individual atom can compensate for changes in each other, making it difficult to pinpoint the “correct” LJ parameters of a molecule based on reproduction of condensed phase properties alone.( 57 ) To overcome this problem a method has been developed to determine the relative value of the LJ parameters based on high level QM data ( 58 ) and the absolute values based on localized ε ij and R min scans for the reproduction of experimental data.( 59 , 60 ) This approach requires supramolecular interactions between rare gases and the model compound. However, once satisfactory LJ parameters are optimized for atoms in a class of functional groups they can often be directly transferred to other molecules with those functional groups without further optimization.

Reproduction of experimental hydration free energies reflects how well the electrostatic and vdW parameters reproduce interactions with bulk water. This property is especially important in the determination of the free energy of ligand binding, where the solvation free energy of the ligand plays a major role in the overall binding affinity. Recently, in the context of the polarizable Drude force field, it has been shown that atom-pair specific LJ parameters need to be used in order to minimize discrepancies between calculated and experimental hydration free energies while simultaneously reproducing pure solvent heats of vaporization and molecular volumes. This can be achieved through the NBFIX facility in CHARMM, and consists of overwriting well-depths and vdW radii based on the standard combination rules with atom-pair specific Lennard-Jones parameters (e.g. atom -O w ) that are optimized to reproduce the aqueous free energies ( 19 ). As shown in Table I this has led to considerably better agreement with experimental values for free energies of hydration calculated with the Drude polarizable force field versus those calculated with the additive force field. The need for atom-pair specific LJ parameters is due to limitations in the combining rules currently use; ( 19 ) efforts are ongoing in our laboratory to identify new LJ combining rules such that atom-pair specific parameters will not be required.

Optimization of the internal parameters in the CHARMM force field relies heavily on QM calculations and data from the Cambridge Crystallographic Database ( 61 ). The parametrization effort consists of determining appropriate force constants and equilibrium values for the bond, angle, improper and Urey-Bradley terms and force constants, multiplicities and phase angles for the dihedral term. Target data includes geometrical information, supplemented by vibrational spectra, including assignments of the normal mode contributions in terms of primitive modes.( 62 , 63 ) In compounds with rotatable bonds, conformational energies as a function of rotation around those bonds are used as additional target data for the optimization of the associated dihedral parameters ( 64 , 65 ). Examples include the peptide bond and the glycosidic linkage ( 11 , 12 ). An automated fitting program is currently available ( 66 ) and is especially useful for large multidimensional target data sets such as those generated for glycosidic linkages between monosaccharides ( 11 ).

It should again be noted that the parameter optimization process follows an iterative approach ( Figure 2 ). Thus, if changes are made to one aspect of the parameters for a model compound, other parameters have to be re-evaluated with respect to reproduction of the target data and optimized as required. While this may be considered tedious, given the now extensive nature of the CHARMM force fields, this procedure typically only requires two or three cycles to reach an adequate level of convergence.

Recent developments and applications of the CHARMM additive force field

While a range of parameters in the context of the additive force field are available, continuous optimization and development efforts have improved the accuracy and expanded the coverage of the CHARMM additive force field allowing more accurate simulations on heterogeneous molecular systems. Current versions of the additive force field include CHARMM22/CMAP ( 39 , 67 ) for proteins, CHARMM27( 68 ) for nucleic acids, CHARMM36( 9 , 10 ) for lipids, CHARMM37( 11 – 15 ) for carbohydrates, and CGenFF( 16 ) for small molecules. As discussed below, updates of the CHARMM22/CMAP protein and CHARMM27 nucleic acid parameters are nearing completion as of the writing of this article. The parameters are available at the MacKerell group web page at ( http://mackerell.umaryland.edu/CHARMM_ff_params.html ) as well as being included with the CHARMM package. It should be emphasized that these different parameter sets are subgroups of the unified CHARMM additive force field as all the contained parameters are compatible, a prerequisite for calculations of complex heterogeneous systems. Recent adjustments to the organization of the individual files will facilitate the application of the CHARMM force field to such systems. This section will provide information complementary to previous published reviews ( 33 , 39 ) on recent improvements and extension in the additive force field in CHARMM. Discussions about the current progress in the CHARMM polarizable force field will ensue.

Polypeptides and Proteins (CMAP)

The protein φ/ψ parameters are crucial for accurate simulation and sampling of the conformational properties of polypeptides and protein. Studies from around the turn of the millennium indicated that the intrinsic form of the dihedral term in the potential energy function limits the ability of force fields in reproducing QM data and experimental structures of proteins ( 69 ) as well as the ability to fold small peptides. To address this issue, improvements in the backbone energetics of the CHARMM22 force field were undertaken in the form of a 2D dihedral energy correction map (CMAP) that could be applied to the φ/ψ surface.( 67 ) Since alanine is the common denominator for all natural amino acids with the exception of glycine and proline, CMAP was calculated based on QM φ/ψ energy surfaces of the alanine dipeptide. The correction grid applies to all residues in the CHARMM22 force field except glycine and proline for which separate CMAPs were implemented. This alanine dipeptide-based CMAP was further adjusted empirically through MD simulations of protein crystals and small peptides, with the latter using an implicit solvent model, bringing significant improvements in reproducing crystallographic φ/ψ distributions.( 70 ) Furthermore, a study on hen lysozyme has shown that CMAP corrections significantly reduce the discrepancies in the backbone dynamics between CHARMM22 and NMR observations.( 71 ) The force field containing this grid-based energy correction scheme is referred to as CHARMM22/CMAP.

The CHARMM22/CMAP force field has been used to study the dynamics of various systems in contexts ranging from structure-function analysis of proteins as well as drug development. Guvench et al ( 72 ) performed MD simulations to elucidated the structural basis for substrate binding in SHP-2, a non-receptor tyrosine phosphatase. It was shown that substrate pY-peptide binding is modulated by Tyr66, which, through conformational flexibilities in its backbone, acts as a gate for the binding cleft. CHARMM was also used to derive the mechanisms for binding of the anticancer drug imatinib to its receptor.( 73 ) In this study, a combination of CHARMM22/CMAP and newly developed parameters for imatinib were used to conduct free energy calculations yielding good agreement with experimental binding affinities and allowed the researchers to suggest imatinib modifications with improved binding. Long MD simulations on the order of hundreds of nanoseconds are often required to trace subtle structural changes in proteins that could play pivotal roles in activity. Such is the case for the recent structure-activity relationship (SAR) assessment of the naturally occurring small protein PcFK1 that displays antimalarial properties.( 74 ) Correlated movements of the loop, turns, and sheet regions on the protein were identified through 150 ns simulations of wild-type and mutant forms of PcFK1 using the CHARMM22/CMAP force field, from which a SAR model was developed. This diverse set of examples demonstrates the breadth of applications possible with and the fidelity of simulation results offered by the CHARMM22/CMAP force field.

The ability to accurately reproduce experimental structures and energies is essential for successful applications of force fields in MD simulations. One major issue with CHARMM22/CMAP is that although CMAP corrections gave better agreement with NMR measurements for folded proteins, simulations of unstructured peptides such as Ala 5 is biased towards α-helical conformations.( 75 ) Efforts to address the problem including revisions of the CMAP and amino acid side-chain torsions are ongoing (R. Best, M. Feig, X. Zhu and A. MacKerell, Personal communication). Tests of the new parameters are currently pending completion and release of these parameters is anticipated by the time this article is published. Improvements of secondary structure and sidechain torsion sampling brought to the CHARMM protein force field will likely improve the quality of MD simulations for the study of both the structural and dynamical properties of proteins.

Nucleic Acids

The CHARMM27 nucleic acid force field( 68 ) was developed to fully support computational analysis of DNA and RNA. It consists of a full reoptimization of an earlier version of the CHARMM additive force field for nucleic acids that was shown to favor the A form of DNA.( 76 ) CHARMM27 has been used in a study by Priyakumar and MacKerell ( 77 ) to investigate DNA base-flipping, which has provided a detailed understanding of the mechanisms and free energy profile for the breaking of Watson-Crick base-pairing and base stacking interactions associated with a C:G base flip. In another study, DNA:RNA hybrid systems were simulated and shown to yield good agreement with experimental observations ( 78 ). A separate study examined the structural basis of DNA triplex formation ( 79 ) and the impact of base modification on triplex sequence specificity. Results from the simulation will aid the design of efficient triple helix-forming oligonucleotides useful for gene targeting. One issue that remains to be addressed is the overestimation of B I /B II ratio for DNA duplexes such that unrestrained simulations fail to reproduce the sequence-dependent B II populations, a problem also seen in other force fields.( 80 ) This offers room for further improvements to the CHARMM27 force field so to expand its utility in the study of the conformation, energy, and sequence-dependent behavior of nucleic acids; efforts to address the B I /B II issue are ongoing (L. Nilsson and A. MacKerell, Personal communication)

CHARMM27 has also been shown to be of utility in the context of RNA. One recent example is a study of the adenine riboswitch ( 81 ) as well as the DNA:RNA hybrid study mentioned above. Another study noted the presence of local opening of WC base pairs in duplex RNA, an phenomena indicated to be consistent with experimental data for GC base pairs but not for AU pairs ( 82 ). Other studies include MD simulations of a protein-RNA complex ( 83 ), structural characterization of the locked nucleic acid-RNA duplexes ( 84 ), and free energy profile of base-flipping in dsRNA ( 85 ). Two recent studies in which both the CHARMM27 and AMBER (param99) force fields were used are of particular interest. In a study investigating folding of two RNA hairpins based on potential of mean force calculations ( 86 ) it was observed that the AMBER force field significantly overestimated the experimental folding energy while CHARMM27 was in significantly better agreement with experiment, though one of the hairpins was indicated to be slightly unstable at room temperature. A second study involved a RNA 18-mer duplex ( 87 ) that showed the RMS fluctuations in MD simulations using AMBER to be significantly lower than that of CHARMM27, with the AMBER model staying closer to the canonical starting structure. To understand these phenomena, as well as local base opening seen in these CHARMM27 RNA simulations, systematic analysis was undertaken from which the sampling of the 2′OH moiety was identified as a possible cause of the local opening. Adjustments to the corresponding parameters were made and a number of models were obtained in which base opening was decreased ( 20 ). Notably, the new models lead to significantly improved reproduction of a range of experimental data, including better agreement with crystal and NMR structures, while still yielding free energies of folding of the RNA hairpins studied by Deng and Cieplak ( 86 ) similar to those calculated with CHARMM27, indicating the adjustments to not over stabilize the folded state. It is anticipated that the optimized parameter set will be referred to as CHARMM36.

The CHARMM27r( 88 ) lipid force field has been used for the study of biological systems involving the lipid bilayer. Recent publications include, for example, elucidation of the free energy profile and barriers for cholesterol flip-flops in three different types of lipid environments( 89 ), simulation of a peptide in lipid bilayers( 90 ) to better understand the dynamics of peptides as a function of membrane type and thickness, and exploration of the mechanism behind which a serine protease is anchored to a bilayer( 91 ). Success of these studies is associated, in part, with revisions of the original CHARMM27 force field. In the revision parameters for the aliphatic tails were optimized using high level QM calculations as target data to yield better agreements with NMR order parameters for the aliphatic tails. However, one major flaw exists in the CHARMM27r version where positive surface tension of the bilayer results in shrinking of the membrane to a gel-phase state well above the transition temperature.( 92 ) This flaw has been recently addressed by Klauda et al. ( 92 ) where headgroup torsions parameters and LJ and partial atomic charges of the ester linkage were re-optimized to yield the CHARMM36 lipid force field. This version reproduces experimental surface areas for saturated and unsaturated lipid bilayers, which prevents unwarranted shrinking of the lipid bilayer. One issue that remains to be addressed is that of the treatment of long-range LJ interactions, which, as discussed by the authors of the CHARMM36 article, is different for monolayers than with bilayers, with the latter being correctly treated with standard truncation methods common to all simulation programs. Accordingly, the CHARMM36 version of the lipid force field is suitable for the majority of the current research interests involving lipids including heterogeneous systems involving bilayers.

Carbohydrates

The additive CHARMM force field has historically contained a limited number of parameters for carbohydrates( 93 ). It has recently been extended to include a range of carbohydrates, an omnipresent entity in cells that is essential in many cellular processes. Modeling of carbohydrates is further driven by scientific interest in understanding their biological and physiological functions in complex systems involving glycoconjugates and carbohydrate-protein complexes. In contrast to proteins, carbohydrates lack clear structural classifications.( 94 ) While some polysaccharides have a regular repeating unit that favors structural periodicities, many are widely diverse, varying in the context of multiple factors such as the types of constituent monosaccharides and the type of the glycosidic linkages, variations that effect their ability to form regular structures and repeating networks of hydrogen bonds formed by the constituent hydroxyl groups.( 95 – 97 ) Such structural diversity presents a challenge for the development of reliable parameters.

In the past few years, development of a carbohydrate force field compatible with the additive CHARMM force field has taken place in our group. The work has been systematic, covering major classes of mono- and disaccharides. As a result parameters are currently available for hexopyranose ( 15 ), aldopentofuranose ( 14 ), and fructofuranose ( 14 ) monosaccharides and as well as an extensive set of glycosidic linkage types ( 11 , 12 ) as required to study di- and polysaccharides. In addition, optimization of parameters for non-hydroxyl moieties present in eukaryotic carbohydrates is nearing completion.

Development of CHARMM force field parameters for cyclic hexopyranoses and furanoses ( Figure 6 ) and acyclic polyalcohols, acyclic carbohydrates, and inositol followed a common approach, including the usual array of target data: QM optimized molecular geometries, experimental geometries derived from small crystals, vibrational frequencies, solute-water interaction energies, molecular volumes, and heats of vaporization or sublimation. One of the more time consuming stages is the extensive QM sampling of conformational space of the furanose monosaccharides was for the two-dimensional ring pucker energy surfaces. Parametrization of acyclic polyalcohols, acyclic carbohydrates, and inositol had been complicated by the intramolecular hydrogen bonds and by the increased flexibility of the compounds versus their cyclic counterparts. The strategy employed to solve this problem was to perform constrained optimizations, with a total of ~1,500 conformations considered. In all cases, simultaneous consideration of a vast amount of QM data presents a challenge. To deal with that an automated scheme( 66 ), based on Monte Carlo simulated annealing (MCSA), was used to accommodate the large amount of data. The performance of the newly developed force field parameters was evaluated by comparing experimental and calculated condensed-phase properties from MD simulations. Good agreement was found for crystal lattice parameters and densities, diffusion coefficients as well as NMR proton-proton couplings in aqueous solution in a subsequent study( 98 ) and is indicative of the quality of the CHARMM force field parameters for carbohydrates.

Cyclic hexopyranose and furanose compounds parametrized for the CHARMM carbohydrate force field.

Parameters for Small Molecules: Applications in Drug Design

An important utility of empirical force fields is to aid in the design and development of new therapeutic agents.( 99 , 100 ) Because of the extent of structural and chemical diversity in drug molecules, parameters geared towards transferability across a wide range of compounds is required and sought after by many force field development communities. Fulfilling this purpose are specialized force field parameters including MMFF ( 101 ), AMBER GAFF ( 102 ), CVFF ( 103 ), COMPASS ( 104 ) and the commercial version of CHARMm ( 105 ). These force fields are appropriate, for example, for methods that identify biologically active or “lead” compounds from a large database of small molecules as well as computer-based de novo molecular design of drug-like molecules. Using the aforementioned force fields new drug entities can be readily assembled, analyzed, and subjected to computational predictions of activity. However, it should be emphasized that reliable calculations such as those required for target-based structure optimization of lead compounds require the use of a well-balanced force field that properly treats both the drug and target (e.g. protein) molecules. In addition, to accurately capture the subtle differences across a series of test compounds it is generally necessary that at least some level of parameter optimization of the molecules of interest be performed, due to the inherent limited transferability of parameters within a force field. With sufficient optimization, a force field can provide superior accuracy in conformational sampling and energy evaluations, leading to improved results from a variety of techniques such as 4D-QSAR( 106 ), conformationally-sampled pharmacophore (CSP) ( 107 ) and related methods, as well as free energy perturbation calculations ( 108 ).

The increasing need for computational analysis of drug-like molecules was the driving force behind the development of the CHARMM General Force field (CGenFF) ( 16 ). Over 100 common and medicinally useful molecular scaffolds and linkers have been formally added to the existing collection of model compounds in the CHARMM force field. From this pool of model compounds, researchers can intuitively assemble compounds of interest, followed by evaluation of the quality of the resulting model and optimization of selected parameters to improve the accuracy of the treatment of the molecule of interest. Typically, one would only need to optimize the linker regions between ring systems. A summary of the steps involved in creating a new molecule is shown in Figure 7 . CGenFF was conceived using a methodology consistent with the parametrization philosophy behind the biological CHARMM additive force field. As a result, CGenFF is fully compatible with the biological portions of the force field, which makes possible simulations of heterogeneous systems such as protein-ligand interactions. This concept is similar to the philosophy used in OPLS ( 109 , 110 ). For complex systems or cases where a large number of molecules need to be studied, scripts that generate initial guesses of the atom type and parameters are available on the MacKerell group web page ( http://mackerell.umaryland.edu/ ) and efforts are ongoing with respect to the development of a web server that automatically assigns atom types, charges and parameters ( https://www.paramchem.org/ ). The ultimate goal of that server will be the ability to validate and optimize selected parameters, a capability that is anticipated to facilitate the expansion of CGenFF to cover a larger range of chemical space. The availability of CGenFF will likely spur additional interest in computer-aided drug design using the CHARMM force field.

Flowchart for the determination of parameters for drug-like molecules as described in reference 16 . 3-phenoxymethylpyrrolidine is assembled from two of its constituents, ethoxybenzene and 3-hydroxymethyltetrahydrofuran, available from the CGenFF. Parameters were identified for the pyrrolidine group by analogy. Optimization of the dihedrals I, II, and III is required to produce an accurate computational model for this molecule.

Recent advances in the CHARMM polarizable force field

The initial step in developing the polarizable force field based on the Drude oscillator was the parametrization of water. A four-site polarizable water model SWM4-NDP was developed( 48 ). In this model, Drude particles have negative charges rather than positive, as was in the previous SWM4-DP model ( 111 ). Such representation more accurately reflects the shifting of electron distributions associated with changes in the environment. The convention of using negatively charged Drude particles has been maintained in subsequent development efforts. SWM4-NDP yields a clear improvement over the additive TIP3P ( 112 ) model in that both gas-phase dipole moments and condensed-phase dielectric constant are simultaneously reproduced with high accuracy. The ability of reproducing the orientation and extent of quadrupoles is also a clear indication of the superiority of the polarizable water model over that of the additive.

From the time the 4-point polarizable water model SWM4-NDP was published ( 48 ) until now, significant progress has been made in the development of CHARMM polarizable force fields for proteins, lipids and nucleic acids. Current available model compounds for proteins include the backbone N-methylacetamide ( 113 ) and amino acid sidechain models such as alkanes ( 53 ), alcohols ( 49 ), aromatics ( 50 , 51 ), amides ( 113 ), and sulfur-containing compounds ( 47 ). Compared to the additive force field, polarizable models display increased accuracy in gas-phase dipole moments, dielectric constants of bulk solvents, and free energies of solvation. Finalization of both the peptide backbone and sidechain parameters is ongoing and preliminary simulations of small polypeptides as well as full proteins have been performed. A full polarizable model of DPPC lipids has been created and successfully applied to simulate a bilayer ( 114 ); additional optimization of that model is ongoing. In addition, parameters for nucleic acid bases have recently been reported ( 115 ) and preliminary simulations of DNA and RNA have been performed. Finally, optimization of Drude parameters for carbohydrates has been initiated, such that the availability of a comprehensive Drude polarizable force field for biological molecules in the next several years is anticipated.

Development efforts have also extended to ions, which are not only an integral part of the cellular environment but also an important component of metalloproteins ( 54 , 116 ). One major issue encountered during this optimization was that of polarization catastrophe, as seen in other polarizable models ( 38 , 117 , 118 ). This primarily results from strong Coulombic interactions that overcome LJ repulsion and lead to uncontrolled drifting of the Drude particle away from the atomic center when atoms approach to closely, a problem enhanced most notably in the case of dications. To address this issue the Drude energy function was extended to include an anharmonic hyperpolarizability damping term

where a force constant K hyp is applied once the Drude displacement, Δ r, is a distance greater than Δ r cut . This term limits the maximum displacement of Drude particles with negligible effects on the overall electrostatic properties. The force constant of 40,000 kcal/mol/Å 4 was selected based on QM-estimated maximum induced dipole of halide anions such as Br − and I − ( 54 ). A Δ r cut of 0.2 Å is currently used for simulations of biomolecules using the Drude polarizable force field in CHARMM. Implementation of the damping function provides an effective way to minimize the polarizability catastrophe issue commonly seen in polarizable force fields.

Conclusion and Perspectives

Molecular mechanics as a tool to investigate experimentally difficult to access properties of molecular systems has gained wide acceptance within the scientific community. After years of effort, the CHARMM additive force field has evolved into a full featured force field capable of treating biologically relevant macromolecules: proteins, lipids nucleic acids and carbohydrates as well as a range of small drug-like molecules. Improvements in this force field, which relies on fine-tuning of parameters to reproduce an ever increasing collection of QM and experimental data, are ongoing.

Aside from fundamental algorithmic changes and advances in computer power, improvements in empirical force field can be made in three tiers. These include modifications to the underlying form of the potential energy function on which the force field is based, expansion of the number of explicitly parametrized model compounds, and additional refinements of available force fields. In this review all three of the issues have been addressed in the context of the CHARMM force field development efforts. Moreover, as computational power continues to increase and conformational sampling methodologies improves, it will be possible to both improve present force fields as well as extend empirical force fields to more complex energy functions, as evidenced by ongoing efforts in the development of polarizable force fields. It is anticipated that the inherent advantage of polarizable force fields over current additive models will allow them to be of significant utility in the application of molecular mechanics approaches to solving biological questions.

As a final word, we must bear in mind that even the most accurate force field is only as good as the computational design with which it is applied. Carefully crafted calculations require an understanding of the strengths and weaknesses of the applied force field. These considerations are also important for interpreting results from the calculations. With this review we hope to provide a basis towards an expanded understanding of force fields with the purpose of facilitating proper implementation of the CHARMM force fields.

Comparison of experimental/QM and empirical gas phase dipole moments, μ, pure solvent dielectric constants, ε, and free energies of hydration (kcal/mol) for various model compounds using the additive and Drude polarizable CHARMM force fields.

ΔG hydr free energies of aqueous hydration in kcal/mol and dipole moments in debeye. QM (versus experimental) dipole moments are shown in italic.

Ref. ( 53 )

Ref. ( 49 )

Ref. ( 47 )

Ref. ( 52 )

Acknowledgments

We would like to thank Dr. Sairam S. Mallajosyula, Dr. Christopher Baker, as well as everyone in the MacKerell group for help with various references, active discussions and moral support. This work was supported by a departmental fellowship from the School of Pharmacy, University of Maryland, Baltimore to XZ, the University of Maryland Computer-Aided Drug Design Center and NIH grants GM070855, GM015101 and GM0772558 and NSF grant CHE-0823198.

• 1. McCammon JA, Gelin BR, Karplus M. Dynamics of folded proteins. Nature. 1977;267(5612):585–590. doi: 10.1038/267585a0. [ DOI ] [ PubMed ] [ Google Scholar ]
• 2. Brooks BR, Brooks CL, 3rd, Mackerell AD, Jr, Nilsson L, Petrella RJ, et al. CHARMM: the biomolecular simulation program. J Comput Chem. 2009;30(10):1545–1614. doi: 10.1002/jcc.21287. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 3. Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M. CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J Comput Chem. 1983;4:187–217. [ Google Scholar ]
• 4. Case DA, TAD, Cheatham TE, III, Simmerling CL, Wang J, Duke RE, Luo R, RCW, Zhang W, Merz KM, Roberts B, Wang B, Hayik S, Roitberg A, Seabra GIK, Wong KF, Paesani F, Vanicek J, Liu J, Wu X, Brozell SR, Steinbrecher T, HG, Cai Q, Ye X, Wang J, Hsieh M-J, Cui G, Roe DR, Mathews DH, MGS, Sagui C, Babin V, Luchko T, Gusarov S, Kovalenko A, Kollman PA. AMBER. University of California; San Francisco: 2010. p. 11. [ Google Scholar ]
• 5. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, et al. Scalable molecular dynamics with NAMD. Journal of Computational Chemistry. 2005;26(16):1781–1802. doi: 10.1002/jcc.20289. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 6. Bjelkmar Pr, Larsson P, Cuendet MA, Hess B, Lindahl E. Implementation of the CHARMM Force Field in GROMACS: Analysis of Protein Stability Effects from Correction Maps, Virtual Interaction Sites, and Water Models. Journal of Chemical Theory and Computation. 2010;6(2):459–466. doi: 10.1021/ct900549r. [ DOI ] [ PubMed ] [ Google Scholar ]
• 7. Chandler DE, Hsin J, Harrison CB, Gumbart J, Schulten K. Intrinsic curvature properties of photosynthetic proteins in chromatophores. Biophys J. 2008;95(6):2822–2836. doi: 10.1529/biophysj.108.132852. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 8. Shaw DE, Deneroff MM, Dror RO, Kuskin JS, Larson RH, et al. Anton, a special-purpose machine for molecular dynamics simulation. SIGARCH Comput Archit News. 2007;35(2):1–12. [ Google Scholar ]
• 9. Klauda JB, Venable RM, Freites JA, O’Connor JW, Tobias DJ, et al. Update of the CHARMM All-Atom Additive Force Field for Lipids: Validation on Six Lipid Types. The Journal of Physical Chemistry B. 2010;114(23):7830–7843. doi: 10.1021/jp101759q. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 10. Feller SE, Yin D, Pastor RW, MacKerell AD., Jr Molecular Dynamics Simulation of Unsaturated Lipid Bilayers at Low Hydration: Parametrization and Comparison with Diffraction Studies. Biophy J. 1997;73:2269–2279. doi: 10.1016/S0006-3495(97)78259-6. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 11. Raman EP, Guvench O, MacKerell AD., Jr CHARMM additive all-atom force field for glycosidic linkages in carbohydrates involving furanoses. J Phys Chem B. 2010;114(40):12981–12994. doi: 10.1021/jp105758h. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 12. Guvench O, Hatcher ER, Venable RM, Pastor RW, Mackerell AD. CHARMM Additive All-Atom Force Field for Glycosidic Linkages between Hexopyranoses. J Chem Theory Comput. 2009;5(9):2353–2370. doi: 10.1021/ct900242e. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 13. Hatcher E, Guvench O, Mackerell AD. CHARMM Additive All-Atom Force Field for Acyclic Polyalcohols, Acyclic Carbohydrates and Inositol. J Chem Theory Comput. 2009;5(5):1315–1327. doi: 10.1021/ct9000608. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 14. Hatcher E, Guvench O, Mackerell AD. CHARMM additive all-atom force field for aldopentofuranoses, methyl-aldopentofuranosides, and fructofuranose. J Phys Chem B. 2009;113(37):12466–12476. doi: 10.1021/jp905496e. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 15. Guvench O, Greene SN, Kamath G, Brady JW, Venable RM, et al. Additive empirical force field for hexopyranose monosaccharides. J Comput Chem. 2008;29(15):2543–2564. doi: 10.1002/jcc.21004. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 16. Vanommeslaeghe K, Hatcher E, Acharya C, Kundu S, Zhong S, et al. CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J Comput Chem. 2010;31(4):671–690. doi: 10.1002/jcc.21367. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 17. Lopes PE, Murashov V, Tazi M, Demchuk E, Mackerell AD., Jr Development of an empirical force field for silica. Application to the quartz-water interface. J Phys Chem B. 2006;110(6):2782–2792. doi: 10.1021/jp055341j. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 18. MacKerell AD, Jr, Feig M, Brooks CL., 3rd Improved treatment of the protein backbone in empirical force fields. J Am Chem Soc. 2004;126(3):698–699. doi: 10.1021/ja036959e. [ DOI ] [ PubMed ] [ Google Scholar ]
• 19. Baker C, Lopes P, Zhu X, Roux Bt, MacKerell A. Accurate Calculation of Hydration Free Energies using Pair-Specific Lennard-Jones Parameters in the CHARMM Drude Polarizable Force Field. Journal of Chemical Theory and Computation. 2010;6(4):1181–1198. doi: 10.1021/ct9005773. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 20. Denning EJ, Priyakumar DU, Nilsson L, Mackerell AD., Jr Impact of hydroxyl sampling on the conformational properties of RNA: Update of the CHARMM all-atom additive force field for RNA. 2011 doi: 10.1002/jcc.21777. In Press. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 21. Alder BJWTE. Phase Transition for a Hard Sphere System. J Chem Phys. 1957;27(1208) [ Google Scholar ]
• 22. Hernandez G, Anderson JS, LeMaster DM. Polarization and polarizability assessed by protein amide acidity. Biochemistry. 2009;48(27):6482–6494. doi: 10.1021/bi900526z. [ DOI ] [ PubMed ] [ Google Scholar ]
• 23. David R, Lide E. Handbook of Chemistry and Physics. 84. CRC Press; Boca Raton, Florida: 2003. [ Google Scholar ]
• 24. Gregory JK, Clary DC, Liu K, Brown MG, Saykally RJ. The Water Dipole Moment in Water Clusters. Science. 1997;275(5301):814–817. doi: 10.1126/science.275.5301.814. [ DOI ] [ PubMed ] [ Google Scholar ]
• 25. Badyal YS, Saboungi ML, Price DL, Shastri SD, Haeffner DR, Soper AK. Electron distribution in water. The Journal of Chemical Physics. 2000;112(21):9206–9208. [ Google Scholar ]
• 26. Silvestrelli PL, Parrinello M. Water Molecule Dipole in the Gas and in the Liquid Phase. Physical Review Letters. 1999;82(16):3308. [ Google Scholar ]
• 27. Rick SW, Stuart SJ. Potentials and algorithms for incorporating polarizability in computer simulations. In: Lipkowitz KB, Boyd DB, editors. Reviews in Computational Chemistry. Wiley-VCH; Hoboken, NJ: 2002. pp. 89–146. [ Google Scholar ]
• 28. Halgren TA, Damm W. Polarizable force fields. Curr Opin Struct Biol. 2001;11(2):236–242. doi: 10.1016/s0959-440x(00)00196-2. [ DOI ] [ PubMed ] [ Google Scholar ]
• 29. Jorgensen WL. Special Issue on Polarization. Journal of Chemical Theory and Computation. 2007;3(6):1877–1877. doi: 10.1021/ct700252g. [ DOI ] [ PubMed ] [ Google Scholar ]
• 30. Lopes P, Roux B, MacKerell A. Molecular modeling and dynamics studies with explicit inclusion of electronic polarizability: theory and applications. Theoretical Chemistry Accounts: Theory, Computation, and Modeling (Theoretica Chimica Acta) 2009;124(1):11–28. doi: 10.1007/s00214-009-0617-x. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 31. Anisimov VM, Lamoureux G, Vorobyov IV, Huang N, Roux B, MacKerell AD. Determination of Electrostatic Parameters for a Polarizable Force Field Based on the Classical Drude Oscillator. J Chem Theory Comput. 2005;1(1):153–168. doi: 10.1021/ct049930p. [ DOI ] [ PubMed ] [ Google Scholar ]
• 32. Lamoureux G, Roux B. Modeling induced polarization with classical Drude oscillators: Theory and molecular dynamics simulation algorithm. The Journal of Chemical Physics. 2003;119(6):3025–3039. [ Google Scholar ]
• 33. MacKerell AD., Jr Empirical force fields for biological macromolecules: overview and issues. J Comput Chem. 2004;25(13):1584–1604. doi: 10.1002/jcc.20082. [ DOI ] [ PubMed ] [ Google Scholar ]
• 34. Thole BT. Molecular polarizabilities calculated with a modified dipole interaction. Chem Phys. 1981;59:341–350. [ Google Scholar ]
• 35. Patel S, Brooks CL. CHARMM fluctuating charge force field for proteins: I parameterization and application to bulk organic liquid simulations. J Comp Chem. 2004;25(1):1–15. doi: 10.1002/jcc.10355. [ DOI ] [ PubMed ] [ Google Scholar ]
• 36. Patel S, Mackerell AD, Brooks CL. CHARMM fluctuating charge force field for proteins: II -Protein/solvent properties from molecular dynamics simulations using a nonadditive electrostatic model. J Comp Chem. 2004;25(12):1504–1514. doi: 10.1002/jcc.20077. [ DOI ] [ PubMed ] [ Google Scholar ]
• 37. Patel S, Brooks CL. Fluctuating charge force fields: Recent developments and applications fromsmall molecules to macromolecular biological systems. Mol Simul. 2006;32(3–4):231–249. [ Google Scholar ]
• 38. Xie W, Pu J, Mackerell AD, Gao J. Development of a polarizable intermolecular potential function (PIPF) for liquid amides and alkanes. J Chem Theory Comput. 2007;3(6):1878–1889. doi: 10.1021/ct700146x. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 39. MacKerell AD, Bashford D, Bellott M, Dunbrack RL, Evanseck JD, et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J Phys Chem B. 1998;102(18):3586–3616. doi: 10.1021/jp973084f. [ DOI ] [ PubMed ] [ Google Scholar ]
• 40. Halgren TA. POTENTIAL-ENERGY FUNCTIONS. Current Opinion in Structural Biology. 1995;5(2):205–210. doi: 10.1016/0959-440x(95)80077-8. [ DOI ] [ PubMed ] [ Google Scholar ]
• 41. Jorgensen WL, Tirado-Rives J. The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin. Journalof the American Chemical Society. 1988;110(6):1657–1666. doi: 10.1021/ja00214a001. [ DOI ] [ PubMed ] [ Google Scholar ]
• 42. Singh UC, Kollman PA. AN APPROACH TO COMPUTING ELECTROSTATIC CHARGES FOR MOLECULES. Journal of Computational Chemistry. 1984;5(2):129–145. [ Google Scholar ]
• 43. Chirlian LE, Francl MM. ATOMIC CHARGES DERIVED FROMELECTROSTATIC POTENTIALS -A DETAILED STUDY. Journal of Computational Chemistry. 1987;8(6):894–905. [ Google Scholar ]
• 44. Merz KM. ANALYSIS OF A LARGE DATA-BASE OF ELECTROSTATIC POTENTIAL DERIVED ATOMIC CHARGES. Journal of Computational Chemistry. 1992;13(6):749–767. [ Google Scholar ]
• 45. Bayly CI, Cieplak P, Cornell W, Kollman PA. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. The Journal of Physical Chemistry. 1993;97(40):10269–10280. [ Google Scholar ]
• 46. Francl MM, Carey C, Chirlian LE, Gange DM. Charges fit to electrostatic potentials .2. Can atomic charges be unambiguously fit to electrostatic potentials? Journal of Computational Chemistry. 1996;17(3):367–383. [ Google Scholar ]
• 47. Zhu X, MacKerell AD., Jr Polarizable empirical force field for sulfur-containing compounds based on the classical Drude oscillator model. J Comput Chem. 2010;31(12):2330–2341. doi: 10.1002/jcc.21527. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 48. Lamoureux G, Harder E, Vorobyov IV, Roux B, MacKerell AD., Jr A polarizable model of water for molecular dynamics simulations of biomolecules. Chem Phys Lett. 2006;418(1–3):245–249. [ Google Scholar ]
• 49. Anisimov VM, Vorobyov IV, Roux B, MacKerell JAD. Polarizable empirical force field for the primary and secondary alcohol series based on the classical Drude model. J Chem Theory Comp. 2007;3:1927–1946. doi: 10.1021/ct700100a. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 50. Lopes PE, Lamoureux G, Roux B, Mackerell AD., Jr Polarizable empirical force field for aromatic compounds based on the classical drude oscillator. J Phys Chem B. 2007;111(11):2873–2885. doi: 10.1021/jp0663614. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 51. Lopes PE, Lamoureux G, Mackerell AD., Jr Polarizable empirical forcefield for nitrogen-containing heteroaromatic compounds based on the classical Drude oscillator. J Comput Chem. 2009;30(12):1821–1838. doi: 10.1002/jcc.21183. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 52. Baker CM, Mackerell AD., Jr Polarizability rescaling and atom-based Thole scaling in the CHARMM Drude polarizable force field for ethers. J Mol Model. 2009 doi: 10.1007/s00894-009-0572-4. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 53. Vorobyov IV, Anisimov VM, MacKerell AD., Jr Polarizable empirical force field for alkanes based on the classical Drude oscillator model. J Phys Chem B. 2005;109(40):18988–18999. doi: 10.1021/jp053182y. [ DOI ] [ PubMed ] [ Google Scholar ]
• 54. Yu H, Whitfield TW, Harder E, Lamoureux G, Vorobyov I, et al. Simulating Monovalent and Divalent Ions in Aqueous Solution Using a Drude Polarizable Force Field. J Chem Theory Comput. 2010;6(3):774–786. doi: 10.1021/ct900576a. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 55. Jorgensen WL. Optimized intermolecular potential functions for liquid alcohols. The Journal of Physical Chemistry. 1986;90(7):1276–1284. [ Google Scholar ]
• 56. Jorgensen WL, Madura JD, Swenson CJ. Optimized intermolecular potential functions for liquid hydrocarbons. Journal of the American Chemical Society. 1984;106(22):6638–6646. [ Google Scholar ]
• 57. MacKerell AD., Jr . Atomistic Models and Force Fields. In: Becker OM, MacKerell AD Jr, Roux B, Watanabe M, editors. Computational Biochemistry and Biophysics. Marcel Dekker, Inc; New York: 2001. pp. 7–38. [ Google Scholar ]
• 58. Yin DX, MacKerell AD. Ab initio calculations on the use of heliumand neon as probes of the van der Waals surfaces of molecules. Journal of Physical Chemistry. 1996;100(7):2588–2596. [ Google Scholar ]
• 59. Yin DX, Mackerell AD. Combined ab initio empirical approach for optimization of Lennard-Jones parameters. Journal of Computational Chemistry. 1998;19(3):334–348. doi: 10.1002/jcc.1166. [ DOI ] [ PubMed ] [ Google Scholar ]
• 60. Chen IJ, Yin DX, MacKerell AD. Combined ab initio/empirical approach for optimization of Lennard-Jones parameters for polar-neutral compounds. Journal of Computational Chemistry. 2002;23(2):199–213. doi: 10.1002/jcc.1166. [ DOI ] [ PubMed ] [ Google Scholar ]
• 61. Allen F. The Cambridge Structural Database: a quarter of a million crystal structures and rising. Acta Crystallographica Section B. 2002;58(3 Part 1):380–388. doi: 10.1107/s0108768102003890. [ DOI ] [ PubMed ] [ Google Scholar ]
• 62. Bright Wilson JE, Decius JC, Cross PC. Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra. Dover Publications; 1980. [ Google Scholar ]
• 63. Pulay P, Fogarasi G, Pang F, Boggs JE. Systematic ab initio gradient calculation of molecular geometries, force constants, and dipole moment derivatives. Journal Of The American Chemical Society. 1979;101(10):2550–2560. [ Google Scholar ]
• 64. Beachy MD, Chasman D, Murphy RB, Halgren TA, Friesner RA. Accurate ab Initio Quantum Chemical Determination of the Relative Energetics of Peptide Conformations and Assessment of Empirical Force Fields. Journal of the American Chemical Society. 1997;119(25):5908–5920. [ Google Scholar ]
• 65. Brooks C, Case DA. Simulations of peptide conformational dynamics and thermodynamics. Chemical Reviews. 1993;93(7):2487–2502. [ Google Scholar ]
• 66. Guvench O, MacKerell AD., Jr Automated conformational energy fitting for force-field development. J Mol Model. 2008;14(8):667–679. doi: 10.1007/s00894-008-0305-0. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 67. Mackerell AD, Jr, Feig M, Brooks CL., 3rd Extending the treatment of backbone energetics in protein force fields: limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J Comput Chem. 2004;25(11):1400–1415. doi: 10.1002/jcc.20065. [ DOI ] [ PubMed ] [ Google Scholar ]
• 68. MacKerell AD, Jr, Banavali N, Foloppe N. Development and current status of the CHARMM force field for nucleic acids. Biopolymers. 2000;56(4):257–265. doi: 10.1002/1097-0282(2000)56:4<257::AID-BIP10029>3.0.CO;2-W. [ DOI ] [ PubMed ] [ Google Scholar ]
• 69. Feig M, MacKerell AD, Brooks CL. Force Field Influence on the Observation of π-Helical Protein Structures in Molecular Dynamics Simulations. The Journal of Physical Chemistry B. 2003;107(12):2831–2836. [ Google Scholar ]
• 70. Price DJ, Brooks CL., 3rd Modern protein force fields behave comparably in molecular dynamics simulations. J Comput Chem. 2002;23(11):1045–1057. doi: 10.1002/jcc.10083. [ DOI ] [ PubMed ] [ Google Scholar ]
• 71. Buck M, Bouguet-Bonnet S, Pastor RW, MacKerell AD., Jr Importance of the CMAP correction to the CHARMM22 protein force field: dynamics of hen lysozyme. Biophys J. 2006;90(4):L36–38. doi: 10.1529/biophysj.105.078154. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 72. Guvench O, Qu CK, MacKerell AD., Jr Tyr66 acts as a conformational switch in the closed-to-open transition of the SHP-2 N-SH2-domain phosphotyrosine-peptide binding cleft. BMC Struct Biol. 2007;7:14. doi: 10.1186/1472-6807-7-14. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 73. Aleksandrov A, Simonson T. A molecular mechanics model for imatinib and imatinib:kinase binding. J Comput Chem. 2010;31(7):1550–1560. doi: 10.1002/jcc.21442. [ DOI ] [ PubMed ] [ Google Scholar ]
• 74. Gleeson MP, Deechongkit S, Ruchirawat S. Molecular dynamics investigation of psalmopeotoxin I. Probing the relationship between 3D structure, anti-malarial activity and thermal stability. J Mol Model. 2010 doi: 10.1007/s00894-010-0732-6. [ DOI ] [ PubMed ] [ Google Scholar ]
• 75. Best RB, Buchete NV, Hummer G. Are current molecular dynamics force fields too helical? Biophys J. 2008;95(1):L07–09. doi: 10.1529/biophysj.108.132696. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 76. MacKerell AD, Wiorkiewicz-Kuczera J, Karplus M. An all-atom empirical energy function for the simulation of nucleic acids. Journal of the American Chemical Society. 1995;117(48):11946–11975. [ Google Scholar ]
• 77. Priyakumar UD, MacKerell AD. Base flipping in a GCGC containing DNA dodecamer: A comparative study of the performance of the nucleic acid force fields, CHARMM, AMBER, and BMS. Journal of Chemical Theory and Computation. 2006;2(1):187–200. doi: 10.1021/ct0501957. [ DOI ] [ PubMed ] [ Google Scholar ]
• 78. Priyakumar UD, Mackerell AD., Jr Atomic detail investigation of the structure and dynamics of DNA.RNA hybrids: a molecular dynamics study. J Phys Chem B. 2008;112(5):1515–1524. doi: 10.1021/jp709827m. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 79. Semenyuk A, Darian E, Liu J, Majumdar A, Cuenoud B, et al. Targeting of an interrupted polypurine:polypyrimidine sequence in Mammalian cells by a triplex-forming oligonucleotide containing a novel base analogue. Biochemistry. 2010;49(36):7867–7878. doi: 10.1021/bi100797z. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 80. Heddi B, Foloppe N, Oguey C, Hartmann B. Importance of accurate DNA structures in solution: the Jun-Fos model. J Mol Biol. 2008;382(4):956–970. doi: 10.1016/j.jmb.2008.07.047. [ DOI ] [ PubMed ] [ Google Scholar ]
• 81. Priyakumar UD, MacKerell AD., Jr Role of the Adenine Ligand on the Stabilization of the Secondary and Tertiary Interactions in the Adenine Riboswitch. Journal of Molecular Biology. 2010;396(5):1422–1438. doi: 10.1016/j.jmb.2009.12.024. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 82. Pan Y, MacKerell AD., Jr Altered structural fluctuations in duplex RNA versus DNA: a conformational switch involving base pair opening. Nucleic Acids Research. 2003;31(24):7131–7140. doi: 10.1093/nar/gkg941. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 83. MacKerell AD, Jr, Nilsson L. Molecular dynamics simulations of nucleic acid-protein complexes. Current Opinion in Structural Biology. 2008;18(2):194–199. doi: 10.1016/j.sbi.2007.12.012. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 84. Pande V, Nilsson L. Insights into structure, dynamics and hydration of locked nucleic acid (LNA) strand-based duplexes from molecular dynamics simulations. Nucleic Acids Res. 2008;36(5):1508–1516. doi: 10.1093/nar/gkm1182. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 85. Hart K, Nystrom B, Ohman M, Nilsson L. Molecular dynamics simulations and free energy calculations of base flipping in dsRNA. RNA. 2005;11(5):609–618. doi: 10.1261/rna.7147805. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 86. Deng N-J, Cieplak P. Free Energy Profile of RNA Hairpins: A Molecular Dynamics Simulation Study. Biophysical Journal. 2010;98(4):627–636. doi: 10.1016/j.bpj.2009.10.040. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 87. Faustino I, Pérez A, Orozco M. Toward a Consensus View of Duplex RNA Flexibility. Biophysical Journal. 2010;99(6):1876–1885. doi: 10.1016/j.bpj.2010.06.061. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 88. Klauda JB, Brooks BR, MacKerell AD, Jr, Venable RM, Pastor RW. An ab initio study on the torsional surface of alkanes and its effect on molecular simulations of alkanes and a DPPC bilayer. J Phys Chem B. 2005;109(11):5300–5311. doi: 10.1021/jp0468096. [ DOI ] [ PubMed ] [ Google Scholar ]
• 89. Jo S, Rui HA, Lim JB, Klauda JB, Im W. Cholesterol Flip-Flop: Insights from Free Energy Simulation Studies. Journal of Physical Chemistry B. 2010;114(42):13342–13348. doi: 10.1021/jp108166k. [ DOI ] [ PubMed ] [ Google Scholar ]
• 90. Rui HA, Im W. Protegrin-1 Orientation and Physicochemical Properties in Membrane Bilayers Studied by Potential of Mean Force Calculations. Journal of Computational Chemistry. 2010;31(16):2859–2867. doi: 10.1002/jcc.21580. [ DOI ] [ PubMed ] [ Google Scholar ]
• 91. Broemstrup T, Reuter N. How does proteinase 3 interact with lipid bilayers? Phys Chem Chem Phys. 2010;12(27):7487–7496. doi: 10.1039/b924117e. [ DOI ] [ PubMed ] [ Google Scholar ]
• 92. Klauda JB, Venable RM, Freites JA, O’Connor JW, Tobias DJ, et al. Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types. J Phys Chem B. 2010;114(23):7830–7843. doi: 10.1021/jp101759q. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 93. Ha SN, Giammona A, Field M, Brady JW. A revised potential-energy surface for molecular mechanics studies of carbohydrates. Carbohydr Res. 1988;180(2):207–221. doi: 10.1016/0008-6215(88)80078-8. [ DOI ] [ PubMed ] [ Google Scholar ]
• 94. Woods RJ. Three-dimensional structures of oligosaccharides. Curr Opin Struct Biol. 1995;5(5):591–598. doi: 10.1016/0959-440x(95)80049-2. [ DOI ] [ PubMed ] [ Google Scholar ]
• 95. Jeffrey GA, Saenger W. Hydrogen Bonding in Biological Structures. Springer-Verlag; Berlin: 1991. [ Google Scholar ]
• 96. Schmidt RK, Tasaki K, Brady JW. Computer Modeling Studies of the Interacton of Water with Carbohydrates. Journal of Food Engineering. 1994;22:43–57. [ Google Scholar ]
• 97. Kirschner KN, Woods RJ. Solvent interactions determine carbohydrate conformation. PNAS. 2001;98:10541–10545. doi: 10.1073/pnas.191362798. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 98. Hatcher E, Säwén E, Widmalm G, MacKerell JAD. Conformational properties of methyl β-maltoside and methyl α-and β-cellobioside disaccharides. J Phys Chem B. 2010 doi: 10.1021/jp109475p. In press. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 99. Young DC. Molecular Mechanics. John Wiley & Sons, Inc; 2009. [ Google Scholar ]
• 100. Wong CF, McCammon JA. Protein simulation and drug design. Adv Protein Chem. 2003;66:87–121. doi: 10.1016/s0065-3233(03)66003-1. [ DOI ] [ PubMed ] [ Google Scholar ]
• 101. Halgren TA. MMFF VII. Characterization of MMFF94, MMFF94s, and Other Widely Available Force Fields for Conformational Energies and for Intermolecular-Interaction Energies and Geometries. J Comp Chem. 1999;20:730–748. doi: 10.1002/(SICI)1096-987X(199905)20:7<730::AID-JCC8>3.0.CO;2-T. [ DOI ] [ PubMed ] [ Google Scholar ]
• 102. Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA. Development and testing ofa general amber force field. J Comput Chem. 2004;25(9):1157–1174. doi: 10.1002/jcc.20035. [ DOI ] [ PubMed ] [ Google Scholar ]
• 103. Ewig CS, Berry R, Dinur U, Hill J-R, Hwang M-J, et al. Derivation of Class II Force Fields. VIII. Derivation of a General Quantum Mechanical Force Field for Organic Compounds. J Comp Chem. 2001;22:1782–1800. doi: 10.1002/jcc.1131. [ DOI ] [ PubMed ] [ Google Scholar ]
• 104. Sun H. COMPASS: An ab Inito Force-Field Optimized for Condensed-Phase Applications-Overview with Details on Alkane and Benzene Compounds. J Phys Chem B. 1998;102:7338–7364. [ Google Scholar ]
• 105. Momany FA, Rone R. Validation of the General Purpose QUANTA 3.2/CHARMm Force Field. Journal of Computational Chemistry. 1992;13(7):888–900. [ Google Scholar ]
• 106. Hopfinger AJ, Wang S, Tokarski JS, Jin BQ, Albuquerque M, et al. Construction of 3D-QSAR models using the 4D-QSAR analysis formalism. Journal of the American Chemical Society. 1997;119(43):10509–10524. [ Google Scholar ]
• 107. Bernard D, Coop A, MacKerell AD., Jr 2D conformationally sampled pharmacophore: a ligand-based pharmacophore to differentiate delta opioid agonists from antagonists. J Am Chem Soc. 2003;125(10):3101–3107. doi: 10.1021/ja027644m. [ DOI ] [ PubMed ] [ Google Scholar ]
• 108. Reddy MR, Erion MD. Free Energy Calculations in Rational Drug Design. Springer -Verlag; [ Google Scholar ]
• 109. Price MLP, Ostrovsky D, Jorgensen WL. Gas-phase and liquid-state properties of esters, nitriles, and nitro compounds with the OPLS-AA force field. Journal of Computational Chemistry. 2001;22(13):1340–1352. [ Google Scholar ]
• 110. McDonald NA, Jorgensen WL. Development of an All-Atom Force Field for Heterocycles. Properties of Liquid Pyrrole, Furan, Diazoles, and Oxazoles. The Journal of Physical Chemistry B. 1998;102(41):8049–8059. [ Google Scholar ]
• 111. Lamoureux G, Jr, Roux B. A simple polarizable model of water based on classical Drude oscillators. Journal of Chemical Physics. 2003;119(10):5185–5197. [ Google Scholar ]
• 112. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of Simple Potential Functions for Simulating Liquid Water. Journal of Chemical Physics. 1983;79:926–935. [ Google Scholar ]
• 113. Harder E, Anisimov VM, Whitfield T, MacKerell AD, Jr, Roux B. Understanding the dielectric properties of liquid amides from a polarizable force field. J Phys Chem B. 2008;112(11):3509–3521. doi: 10.1021/jp709729d. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 114. Harder E, Mackerell AD, Jr, Roux B. Many-body polarization effects and the membrane dipole potential. J Am Chem Soc. 2009;131(8):2760–2761. doi: 10.1021/ja806825g. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 115. Baker CM, Anisimov VM, MacKerell AD., Jr Development of CHARMM polarizableforce field for nucleic acid bases based on the classical Drude oscillator model. J Phys Chem B. 2011;115(3):580–596. doi: 10.1021/jp1092338. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ]
• 116. Lamoureux G, Roux B. Absolute Hydration Free Energy Scale for Alkali and Halide Ions Established from Simulations with a Polarizable Force Field. J Phys Chem B. 2006;110:3308–3322. doi: 10.1021/jp056043p. [ DOI ] [ PubMed ] [ Google Scholar ]
• 117. Yu H, van Gunsteren WF. Accounting for polarization in molecular simulation. Computer Physics Communications. 2005;172(2):69–85. [ Google Scholar ]
• 118. Darley MG, Handley CM, Popelier PLA. Beyond Point Charges: Dynamic Polarization from Neural Net Predicted Multipole Moments. Journal of Chemical Theory and Computation. 2008;4(9):1435–1448. doi: 10.1021/ct800166r. [ DOI ] [ PubMed ] [ Google Scholar ]
• View on publisher site
• PDF (1.9 MB)
• Collections

Similar articles

Cited by other articles, links to ncbi databases.

• Download .nbib .nbib
• Format: AMA APA MLA NLM

Add to Collections