Evaluating the London Dispersion Coefficients of Protein Force Fields

Subscriber access provided by Kaohsiung Medical University

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

Evaluating the London Dispersion Coefficients of Protein Force Fields Using the Exchange-Hole Dipole Moment Model Evan Walters, Mohamad Mohebifar, Erin R Johnson, and Christopher N. Rowley J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b02814 • Publication Date (Web): 07 Jun 2018 Downloaded from http://pubs.acs.org on June 7, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

TOC Graphic 88x48mm (300 x 300 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Evaluating the London Dispersion Coefficients of Protein Force Fields Using the Exchange-Hole Dipole Moment Model Evan T. Walters,† Mohamad Mohebifar,† Erin R. Johnson,‡ and Christopher N. Rowley∗,† †Department of Chemistry, Memorial University of Newfoundland, St. John’s, Newfoundland and Labrador, Canada ‡Department of Chemistry, Dalhousie University, Halifax, Nova Scotia, Canada E-mail: [email protected] Phone: +1 (709) 864-2229

1


Page 2 of 43

Page 3 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Abstract London dispersion is one of the fundamental interactions involved in protein folding and dynamics. The popular CHARMM36, Amber ff14sb, and OPLS-AA force fields represent these interactions through the C6 /r6 term of the Lennard-Jones potential, where the C6 parameters are assigned empirically. Here, dispersion coefficients of these three force fields are shown to be roughly 50% larger than values calculated using the quantum-mechanically-derived exchange-hole dipole moment (XDM) model. The CHARMM36 and Amber OL15 force fields for nucleic acids also exhibit this trend. The hydration energies of the side-chain models were calculated using REMD-TI for the CHARMM36, Amber ff14sb, and OPLS-AA force fields. These force fields predict side-chain hydration energies that are in generally good agreement with the experimental values, which suggests the total strength of aqueous dispersion interactions is correct, despite C6 coefficients that are considerably larger than XDM predicts. An analytical expression for the dispersion hydration energy using XDM coefficients shows that higher-order dispersion terms (i.e., C8 and C10 ) account for roughly 37.5% of the hydration energy of methane. This suggests that the C6 dispersion coefficients used in contemporary force fields are elevated to account for the neglected higher-order terms.

1

Introduction

London dispersion is an attractive force between atoms and molecules that results from fluctuations in the electron density that create instantaneous electric moments. These electric moments induce electric moments in neighboring atoms, resulting in an attractive Coulombic interaction. The potential energy of these interactions can be well-approximated by

Edisp,ij (rij ) = −

Cn,ij , n r ij n=6,8,10... X

(1)

where rij is the distance between atomic pairs and the Cn,ij are coefficients that depend on the identity of the interacting pair. 1 In general, terms of order greater than 10 are negligible. 2



Page 4 of 43

Although the dispersion interaction between a pair of atoms is generally weak and short range, these interactions amount to a significant cohesive force in condensed phases. They are notably significant in determining the structure and function of biomolecules like proteins, nucleic acids, and lipids. 2,3 Dispersion is also an important factor in the stability of crystal phases and interfacial complexes. 4–6 Atomistic simulations of these materials require a molecular-mechanical force field that provides a realistic description of dispersion interactions. Most popular force fields include dispersion interactions through the Lennard-Jones potential, 7 C6,ij Aij − 6 . 12 rij rij

ELJ,ij (rij ) =

(2)

The A/r12 term of this equation is intended to represent Pauli repulsion. Dispersion interactions are represented only by the C6 /r6 term, neglecting the higher-order dispersion terms in Eqn 1. More commonly, this potential is defined in terms of the atomic radius (σ) and well-depth (ε), " ELJ,ij (rij ) = 4εij

σij rij

12

−

σij rij

6 # ,

(3)

where C6,ij = 4εσ 6 . In principle, parameters can be defined for each unique pair of atoms, 8–10 but the most common practice is to define parameters for like-pairs and calculate the parameters for unlike pairs using the Lorentz–Berthelot combination rule, 11

σij =

σii + σjj √ , εij = εii εjj . 2

(4)

The Lennard-Jones parameters are typically assigned empirically such that simulations using these parameters accurately predict the properties of representative bulk liquids (e.g., density, enthalpy of vaporization, dielectric constant, etc.), 10,12–17 although some QM-based approaches have also been applied. 18–22 A drawback of parameterizing Lennard-Jones param-

3


Page 5 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


eters empirically is that they are generally underdetermined when molecules contain multiple types of atoms, allowing widely different parameters to be defined for them. For example, the Amber ff14sb force field uses Lennard-Jones parameters for the thiol sulfur atom type of σ = 3.64 ˚ A and ε = 1.046 kJ/mol, yielding a dispersion coefficient of 149 a.u. 23 For the same atom type, the CHARMM36 force field assigns parameters of σ = 3.52 ˚ A and ε = 1.59 kJ/mol, yielding a dispersion coefficient of 267 a.u. Further, because the underlying force field is approximate, the dispersion parameters can be assigned unphysical values in order to compensate for other issues. In the past, there was no general method to validate the parameters that come out of this process, so they are commonly treated as purely empirical, unrestrained parameters. The exchange-hole dipole moment (XDM) model provides a method for calculating atomic and molecular dispersion coefficients from a simple density-functional theory (DFT) calculation. 24–26 XDM is based on second-order perturbation theory 27,28 and uses the multipole moments for a reference electron and its exchange hole as the source of the instantaneous moments that give rise to dispersion. It allows non-empirical calculation of atomic and molecular dispersion coefficients, to arbitrarily high order, 29 based on properties of the DFT electron density. Because the atomic dispersion coefficients are dependent on the electron density via the exchange hole, they vary with local chemical environment 30 and implicitly include many-body electronic effects 31,32 (or Type A and B non-additive effects in Dobson’s classification scheme 33 ). In general, the computed molecular C6 coefficients are accurate to within ca. 10% when compared to experimental reference data. 34 XDM calculations provide a straightforward way to evaluate force-field dispersion coefficients using quantum chemistry. Mohebifar et al. used XDM to assess the dispersion coefficients of the GAFF, CGenFF, and OPLS-AA force fields for a set of small organic molecules. 35 This report showed that the molecular C6 coefficients for GAFF, CGenFF, and OPLS-AA force fields are systematically 50% higher than the XDM values. This trend is not universal; the C6 coefficient of the TIP3P water model is close to the XDM value, suggesting

4



Page 6 of 43

that the strength of dispersion interactions may be unbalanced in simulations of molecules in solution. In this paper, we extend this analysis to evaluate the dispersion coefficients of molecularmechanical force fields for proteins (Amber ff14sb, 23 CHARMM36, 36 and OPLS-AA 37 ) using XDM. The Amber OL15 and CHARMM36 nucleic acid force fields are also evaluated. We note that many earlier versions of these force fields use the same Lennard-Jones parameters. The relationship between the molecular dispersion coefficient and the computed hydration energy is investigated. Finally, we explore the possible origins and potential resolutions to issues with the treatment of dispersion in molecular-mechanical force fields.

2

Computational Methods

2.1

Side-Chain Models NH2

CH4

N H

methane Ala

O

O

NH2

NH2

1-propylguanidine Arg

acetamide Asn

O

H 3C

O acetic acid Asp

N

propanoic acid Glu O

N-methylacetamide NMA

H 3C toluene Phe

propanamide Gln

NH3

4-methyl-3 H-imidazole Hsd

N H

NH2

methanethiol Cys

HN

O

O

SH

butane Ile

OH

methanol Ser

OH ethanol Thr

2-methylpropane Leu

butan-1-amine Lys

N H 3-methyl-1 H-indole Trp

S methylthioethane Met

OH 4-methylphenol Tyr

propane Val

Figure 1: Structures of compounds used to model amino acid side chains.

Molecular models for the amino acid side chains were defined by deleting the amine and 5


Page 7 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


carboxylic acid moieties and replacing the α-carbon with a hydrogen atom. The hydration energy of N-methylacetamide (NMA), a model for the protein backbone, was also considered. The structures of these species are presented in Figure 1.

2.2

XDM calculations

To calculate the XDM dispersion coefficients, the geometries of all side-chain and nucleobase structures were optimized using DFT calculations with Gaussian 09. 38 The PBE0 exchangecorrelation functional 39 and the aug-cc-pVTZ basis set 40 were used for all calculations. This functional and basis set combination is generally reliable for the calculation of molecular dipole moments and polarizabilities. 41 The XDM dispersion coefficients of each side chain and nucleobase were calculated from the Gaussian wavefunction file using the postg code. 42,43

2.3

Free Energy Simulations

The hydration energies of the side-chain models were calculated with CHARMM 41b1 using the staged Weeks, Chandler, and Andersen protocol of Deng and Roux. 44 The electrostatic component of the Gibbs energy of solvation was calculated from thermodynamic integration (TI) simulations with 11 values of λ spaced evenly between 0 and 1, which corresponds to the solute–solvate electrostatic interaction. The dispersion component of the hydration energy is then calculated by turning off the solute–solvent dispersion interactions. Lastly, the repulsive component was calculated from a 9-stage free energy perturbation (FEP) calculation. A correction for the truncation of long-range Lennard-Jones interactions was obtained by calculating the average of the difference of the solute–solvent Lennard-Jones interactions for a 12 ˚ A cutoff vs a 40 ˚ A cutoff, for a 1 ns simulation of the solute in a 72 ˚ A ×72 ˚ A ×72 ˚ A cubic cell of water. Replica exchange was used to sample the configurational space of the states in the TI/FEP process more effectively by allowing exchanges between neighboring replicas every 1000 time steps, following the implementation of Jiang et al. 45 Uncertainties were estimated by dividing the sampled data into three equal sets and taking the standard 6



deviation of the energies calculated from these values. Water molecules were represented using the TIP3P model. 46 A water–vacuum interfacial potential of -520 mV was included in the calculation of the Gibbs energy of hydration of charged residues. 47

2.4

Molecular Dynamics Calculations

To provide quantitative examples of dispersion interactions, several molecular dynamics simulations were conducted using NAMD 2.12. 48 In these simulations, the protein and lipid components were described using the CHARMM36 force field, while the Tamoxifen and LSD molecules were described using the CGenFF force field. 49 CHARMM-GUI was used to construct a pure POPC bilayer containing 48 lipids. 50 The coordinates for the folded state of the protein barnase were taken from the crystallographic structure (PDB ID: 1BRS). 51 Unfolded states were generated by 10 ns simulations at 600 K, followed by a 1 ns re-equilibration at 298 K with the protein backbone atoms restrained to their unfolded coordinates. The simulations of the LSD–GPCR complex were generated from the crystallographic structure of 5-HT2B reported by Wacker et al (PDB ID: 4IB4). 52 In all cases, the average LennardJones interactions were calculated using the NAMDenergy plugin of VMD from three 2 ns simulations, where coordinates were saved every 20 ps.

2.5

Ab Initio Molecular Dynamics

The ab initio molecular dynamics (AIMD) simulation of aqueous methylthiol was performed using CP2K 2.5.1. 53 The simulation cell contained methylthiol and 162 water molecules in a cubic cell with an edge length of 16.7 ˚ A. The revPBE exchange-correlation functional 54 with the D3 dispersion correction 55 was used, as it has been reported to be effective in describing the structural properties of liquid water. 56 The electron density of the system was represented using a TZVP basis set with a multigrid cutoff of 300 Ry and Goedecker-TeterHutter pseudopotentials. 57 The methylthiol radial distribution function was calculated from a 100 ps simulation that followed a 20 ps equilibration. Temperature was regulated using 7


Page 8 of 43

Page 9 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Langevin dynamics with a friction coefficient of 5 ps−1 . Bonds containing hydrogen were constrained to a distance of 0.96 ˚ A using the SHAKE algorithm. 58

3

Results and Discussion

3.1

The Significance of Dispersion Interactions in Biophysics

When presenting this work, we have been asked if force field dispersion parameters have any significant impact on biophysical processes, so we will begin by highlighting several processes where dispersion interactions are significant. The hydrophobic effect is most commonly discussed in terms of changes in electrostatic interactions and entropy, but changes in dispersion interactions can also play a significant, albeit secondary, role. Water is a polar, hydrogen-bonding liquid, but water molecules are only moderately polarizable, so the dispersion interaction between a molecule and an aqueous solution is relatively weak. Conversely, biomolecules like lipids, proteins, and nucleic acids are largely comprised of more polarizable carbon atoms, which results in stronger dispersion interactions. As a result, there can be a large change in the strength of dispersion interactions when molecules transfer between different cellular environments or a biomolecule changes conformation such that its exposure to aqueous solution changes. For example, a significant change in the strength of dispersion interactions can occur when a protein folds. The protein–solvent dispersion interactions decrease when a protein folds because a high proportion of side-chains are solvent-exposed in unfolded states, but many of these residues are buried in the interior of the protein in the folded states. Conversely, protein–protein dispersion interactions become much stronger when a protein folds, largely because aliphatic hydrophobic residues move into contact in the folded state. Early analysis of protein folding energetics by Honig and coworkers 59 using a continuum model for the solvent concluded that changes in dispersion interactions were not significant in protein folding, although more recent models have shown that more complete descriptions of 8



A

B

Page 10 of 43

C

Figure 2: Examples of biophysical processes that involve changes in dispersion interactions. (a) Hydrocarbon side-chains (purple) of the protein barnase have dispersion interactions with water when unfolded, but form stronger dispersion interactions with each other in the folded state. (b) LSD (green) has stronger dispersion interactions with the amino acid residues of the binding site of a GPCR membrane protein (magenta) than with water molecules in solution. (c) Tamoxifen (green) has stronger dispersion interactions with the acyl chains of a lipid bilayer than with water molecules in solution.

9


Page 11 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


dispersion interactions are needed to achieve quantitative accuracy. 2,60–63 The change in the strength of dispersion interactions can be roughly quantified by calculating the changes in the averages of the protein–protein and protein–water LennardJones interactions in an MD simulation. For example, in a simulation of the protein barnase, the water–protein Lennard-Jones interactions decrease from −533.3 ± 10.0 kcal/mol to −235.4 ± 1.1 kcal/mol when the protein folds, but the protein–protein Lennard-Jones interactions increase from −29.7 ± 9.3 kcal/mol to −339.2 ± 1.6 kcal/mol. In this example, the net result is that the Lennard-Jones interactions are roughly the same in the folded state and unfolded states, although the protein–water and protein–protein dispersion interactions must be carefully balanced to achieve this and proteins will vary in the relative strength of dispersion interactions in their various conformations. Another example where changes in dispersion interactions are significant is the binding of hydrophobic ligands to proteins. In the proteinaceous environment, the ligand will have stronger interactions with the polarizable carbon atoms of the protein than with the aqueous solution from which they originate. For example, the Lennard-Jones interactions of lysergic acid diethylamide (LSD) become stronger when this drug moves from solution into the binding site of a GPCR membrane protein (−31.0 ± 0.1 kcal/mol in solution vs −44.2 ± 0.6 kcal/mol in the binding site). This is in keeping with the analysis of Roux and coworkers, who have shown that increases in dispersion interactions are often the largest component of the Gibbs energy of protein–ligand binding. 64,65 A final example of a biophysical process that involves a large change in dispersion interactions is the permeation of molecules through lipid bilayers. Hydrophobic solutes will generally have stronger dispersion interactions in the membrane interior than they have with the aqueous solutions above and below the bilayer. 66,67 For example, in a simulation of tamoxifen permeating through a POPC bilayer, the total dispersion interactions of tamoxifen become stronger when the molecule moves from solution to the center of the bilayer (−36.0 ± 0.1 kcal/mol in solution vs −46 ± 0.8 kcal/mol in the bilayer).

10



Page 12 of 43

Exptl. XDM SWM6 SWM4-NDP TIP4P-D TIP4P/2005 TIP4P-FB TIP4P-Ew TIP4P TIP3P-FB TIP3P 0

10

20

30

40

50

60

70

C6 (a.u.)

Figure 3: Dispersion coefficients of water models. The reference XDM value is indicated by the dotted line. Values are taken from Ref. 35. The experimental value is taken from Ref. 68. The strength of the dispersion interactions in these simulations is sensitive to the parameters used in the force field. Figure 3 shows the broad range of C6 coefficients used in popular water models that represent dispersion interactions through the single Lennard-Jones term centered on the oxygen atom. For example, the C6 coefficient of the TIP3P 46 model is 43.2 a.u., while the C6 coefficient of the TIP4P/2005 69 model is 53.4 a.u. The large variations in C6 affect the strength of the dispersion interactions between water and the rest of the system. This is apparent in the tamoxifen–bilayer permeation models; when the TIP3P water model is used, the average water–ligand Lennard-Jones interaction of tamoxifen in solution is −36.4 ± 0.1 kcal/mol, but this interaction is strengthened to −38.7 ± 0.1 kcal/mol when the TIP4P/2005 model is used. The stronger drug–water interaction reduces the degree to which dispersion facilitates the transfer of the drug from solution into the bilayer interior. The wide variation of dispersion coefficients in water models that have similar physical properties suggests the potential for inconsistency in these simulations. Each of these processes involves the hydrophobic effect, where there are competing changes in dispersion interactions, electrostatic interactions, and entropy. For simulations to describe these processes rigorously, each component, including dispersion, must be represented accurately. In the most widely-used force fields, the only term that accounts for

11


Page 13 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


dispersion interactions is the C6 /r6 term of the Lennard-Jones potential, so this term and its C6 parameters must describe the dispersion interactions in these systems realistically.

3.2

Molecular Dispersion Coefficients of Amino Acid Side-Chains

One model for understanding the folding and membrane partitioning of proteins is the hydration energies of side-chain models. Molecular models for protein side chains are devised by truncating the side chain with a hydrogen atom that replaces the alpha carbon, as shown in Figure 1. This provides a rough estimate of the thermodynamic cost for the side chain to be exposed to the solvent. Together with N-methyl acetamide (NMA), a model for the amide moiety of the protein backbone, these molecules serve as test systems to evaluate the non-bonded force-field parameters. The atomic and molecular dispersion coefficients of the molecules in Figure 1 were calculated using the CHARMM36, 36 OPLS-AA, 37 and Amber ff14sb, 23 force fields using the relation, C6 = 4εσ 6

(5)

where ε and σ are the atomic Lennard-Jones parameters defined for a given atom in this force field. The molecular dispersion coefficient (C6,mol ) is calculated from the sum of all atomic pairs (C6,ij ) C6,mol =

X

C6,ij .

(6)

ij

The molecular dispersion coefficients of NMA and the side-chain models calculated using XDM are correlated to the force-field dispersion coefficients in Figure 4. All three force fields have molecular C6 dispersion coefficients for the side chain models that are systematically 40–50% higher than the QM values. In our previous study, the molecular C6 coefficients for a set of organic molecules from the Generalized Amber Force Field (GAFF), OPLS-AA, and CGenFF force fields were also found to be roughly 50% larger than the XDM values. 35

12



6000

5000

CHARMM36 (y=1.4x) Amber (y=1.5x) OPLS-AA (y=1.5x)

4000

∑ C6 FF (a.u.)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 43

Arg

3000

2000

1000

0 0

1000

2000

3000

4000

∑ C6 XDM (a.u.)

Figure 4: The correlation between molecular C6 coefficients of the protein side chain models. The dispersion coefficients for arginine are highlighted. The gray line indicates where a 1:1 correlation would lie. The points corresponding to the arginine side chain are highlighted because of their broad range.

13


Page 15 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.3

Atomic Dispersion Coefficients

To identify why the force-field molecular C6 coefficients are so much larger than the XDM values, we compare the individual atomic dispersion coefficients. The correlation between the force field and XDM atomic dispersion coefficients is presented in Figure 5. One immediate observation is that the variation in atomic dispersion coefficients for a given element is much larger for the force fields than XDM predicts. The XDM analysis suggests that the range of atomic dispersion coefficients in molecular systems should be modest.1 Specific comparisons of atomic dispersion coefficients are presented in the following sections.

3.3.1

Carbon

As carbon is the most common non-hydrogen element in the side-chain models, its parameters affect the molecular dispersion coefficients of the side chain-models most significantly. Broadly, the force-field dispersion coefficients for carbon are systematically larger and are distributed over a larger range than the XDM values; the XDM values range from 19 a.u. to 24 a.u.,2 but the CHARMM36 values range from 19 to 65 a.u. This trend is also present in the Amber ff14sb and OPLS-AA force fields, where the ranges for C-atom dispersion coefficients are 38–39 a.u. and 35–84 a.u., respectively. This overestimation is partially countered by the small dispersion coefficients assigned to hydrogen atoms in the force field; the force field C6 coefficients of hydrogen atoms are generally less than 1 a.u., but XDM predicts these values to lie in the 2–3 a.u. range. Nevertheless, the carbon atoms in all residues tend to have dispersion coefficients that are significantly larger than the XDM values. 1 Iterative Hirshfeld approaches have been proposed as an alternative method of performing electron density partitioning, which might lead to a larger variation in atomic dispersion coefficients. 70 2 The carbon atom of the carboxylic acid in Asp and Glu have C6 coefficients of 32 a.u. due to the electron-richness of the anion.

14



C6 CHARMM36 (a.u.)

(a) 300 250 200 150

C H O N S

100 50 0 0

25

50

75

100

75

100

75

100

C6 XDM (a.u.) (b) 300

C6 AMBER (a.u.)

250 200 150

C H O N S

100 50 0 0

25

50

C6 XDM (a.u.) (c) 300 250

C6 OPLS (a.u.)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 43

200 150

C H O N S

100 50 0 0

25

50

C6 XDM (a.u.)

Figure 5: The atomic C6 coefficients of the amino acid side chain models. The gray diagonal line indicates where the force field would match the XDM value.

15


Page 17 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.3.2

Hydrogen-Bonding Atom Types

Atoms that participate in hydrogen-bonding interactions tend to have dispersion coefficients that are much larger than the XDM value. This is pronounced in the CHARMM36 force field, where the dispersion coefficients of all N atoms are assigned values of 74.5 a.u., while the XDM values range between 10.5 a.u. to 19.6 a.u. Similarly, the oxygen atoms of amides have values of 26.9 a.u., while the XDM value is only 15.1 a.u. Hydroxyl oxygens have dispersion coefficients of 43.4 a.u., but the XDM value is only 12.4 a.u. The force field atomic dispersion coefficients of these oxygen and nitrogen atoms are even larger than those of the carbon atoms, although physically, these atoms are significantly less polarizable than carbon atoms, so they should have weaker dispersion interactions. This occurs in all three force fields, but to different degrees. One possible explanation for this trend is the empirical procedure used to parameterize these coefficients, where the Lennard-Jones parameters were adjusted to predict the liquid properties. In some cases, the hydrogen atom is assigned no Lennard-Jones interaction term at all, so the dispersion interactions stemming from the hydrogen atoms must be effectively absorbed into the dispersion-interaction term of their parent atom. 22 Further, because these models neglect induced polarization and charge transfer, the ε and σ terms are heavily adjusted to compensate for these missing associative interactions. This is an unattractive solution, because it could have unintended consequences, such as making these groups spuriously soluble in high-dispersion environments (e.g., the interior of lipid bilayers). To resolve this issue, it would likely be necessary to introduce additional terms to describe the unique electrostatic and repulsive components of hydrogen bonding rigorously. Including induced polarization alone does not appear to resolve this issue because the parameters in the Drude polarizable force field also show this trend. 35 A recent reparameterization of the non-bonded parameters of the GROMOS-compatible parameter set for small organic molecules (2016H66) yielded atomic C6 dispersion coefficients for nitrogen atoms in the 40– 42 a.u. range, indicating that effective force fields can employ significantly smaller dispersion 16



coefficients for nitrogen atoms. 71 3.3.3

Sulfur

The parameters for sulfur atom types found in cysteine and methionine vary significantly between the force fields. The dispersion coefficient for sulfur atom types in the CHARMM36 force field is 267 a.u., while the XDM values are 96 a.u. and 89 a.u., respectively. Similarly, the OPLS-AA force field assigns values of 268 a.u. and 224.3 a.u. to C6 coefficients of the Cys and Met sulfur atoms, respectively. In contrast, the dispersion coefficient assigned to these atom types in the Amber ff14sb force field is only 149 a.u.

3.4

Side-Chain Hydration Energies 5

ΔGcalc (kcal/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 43

0

−5

CHARMM36 Amber ff14sb OPLS-AA

−10 −10

−5

0

5

ΔGexptl (kcal/mol)

Figure 6: Correlation between experimental and calculated hydration energies of the sidechain models for the neutral amino acids. Although it is apparent that the force fields have dispersion coefficients that are higher than the XDM values, the effect of this difference on the calculated properties is less apparent. To explore this, we calculated the absolute Gibbs energy of hydration of the side-chain models in an explicit solvent using the Weeks–Chandler–Andersen (WCA) decomposition scheme 17


Page 19 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of Deng and Roux. 44 This allows the hydration energy to be divided into electrostatic, dispersion, and repulsive components.3 The dispersion component is of particular interest here, as this term results from the interactions between the solute and the solvent through the dispersion component of the Lennard-Jones interactions of the force fields. The correlation between the calculated and experimental hydration energies is presented in Figure 6. Tables of the hydration energies and their decomposition into electrostatic, dispersion, and repulsive components are included in the Supporting Information. In general, the agreement between the side-chain hydration energies and the experimental values is reasonably good for all three force fields, which is consistent with previous reports. 44,72–76 All three force fields exhibit a small bias to underestimate the solubility of the side-chains and the performance of the CHARMM36 force field is incrementally better on the whole. The dispersion component of the hydration energy is generally larger for the CHARMM36 force field than for the Amber ff14sb or OPLS-AA models. Significantly, all three force fields predict reasonably-accurate hydration energies for the aliphatic side-chains (e.g., Ala, Val, Ile, Leu), where the hydration energy is almost exclusively dependent on the balance between the dispersion and repulsive components. This suggests that, although the force field C6 dispersion coefficients are larger than the XDM values, the total solute–water dispersion interactions are probably of approximately the correct strength. The origin and potential solution of this apparent paradox is discussed in Section 3.6.

3.4.1

Side Chain Hydration Energy Prediction

The solvation energy of a hard sphere of radius R that experiences C6 -dispersion interactions with a surrounding solvent can be estimated analytically to be4

∆Gdisp = −4πρw

C6,w−mol . 3R3

3

(7)

This is an inherently path-dependent procedure, so other decomposition schemes will assign different values to the terms of the hydration energy. 4 There are also more complex numerical methods for non-spherical structures. See Ref. 77.

18



The C6,mol terms are for like pairs of side chains, so the dispersion coefficient for the interp action between a water molecule and the side chain would be C6,w−mol = C6,mol C6,w . ρw is

(a)

0

ΔGdisp (kcal/mol)

the number density of liquid water.

−5

CHARMM36 Amber OPLS-AA

−10 −15 −20 0

20

40

60

80

√C6,FF (a.u.)

(b)

0

ΔGdisp (kcal/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 43

−5

CHARMM36 Amber OPLS-AA

C

−10 M

NMA

−15 R

−20 0

20

40

60

80

√C6,XDM (a.u.)

Figure 7: (a) Correlation between the force field C6 molecular dispersion coefficients and the dispersion component of the hydration energy calculated from the free-energy simulations. A strong correlation is observed for all three force fields. (b) Correlation between the XDM C6 molecular dispersion coefficients and the dispersion component of the hydration energy calculated using the free-energy simulations. This correlation illustrates where the forcefield dispersion parameters yield a dispersion solvation energy that is inconsistent with the XDM-predicted dispersion coefficient. This suggests a linear relationship between the square root of the molecular C6 coefficient and the calculated dispersion component of the hydration energy,

∆Gdisp ≈ a

p

19

C6,mol + b.


(8)

Page 21 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


This linear relationship is apparent in Figure 7(a), which correlates the dispersion component of the hydration energy calculated using free-energy simulations with the molecular C6 dispersion coefficient of the force fields. This relation is strongly linear for all three force fields, with coefficients of determination of about 0.99. Section 3.2 showed that there is a roughly linear correlation between the XDM and force-field molecular dispersion coefficients (i.e. C6,FF ∝ C6,XDM ) for the side-chain models. Combining this with Eqn. 8, which relates the molecular dispersion coefficient and dispersion component of the hydration energy, the XDM dispersion coefficients can be used to highlight where the calculated dispersion hydration energy is inconsistent with the XDM dispersion coefficient. The correlations between the XDM molecular dispersion coefficient and the calculated dispersion component of the hydration energy are presented in Figure 7 (b). Points lying far from the line of regression for each force field indicate that the XDM molecular dispersion coefficient would predict a significantly different dispersion component of the Gibbs energy of hydration for that residue. The CHARMM36 models for Arg and NMA are significant outliers, as are the Amber ff14sb models for Cys and Met. The hydration energies of the aliphatic residues (Ala, Leu, Ile, and Val) are generally in good agreement with the experimental values. As these residues have minimal electrostatic interactions with the solvent, the attractive component of the hydration energy is almost entirely due to dispersion interactions. Noting that the force field molecular dispersion coefficient is strongly correlated to the hydration energy and this molecular dispersion coefficient is proportional to the XDM dispersion coefficient, Eqn. 8 can be used to estimate dispersion hydration free energies directly from the XDM dispersion coefficients. Linear regression was used to fit the dispersion component of the CHARMM36 hydration energies of the aliphatic side chains using the XDM molecular dispersion coefficients (a = −0.279, b = −1.44, R2 = 0.9973). This linear relation was used to predict the dispersion contribution to the hydration energy for the remaining side-chain models. This analysis is intended to highlight where

20



the force field dispersion component of the hydration energy is inconsistent with the firstprinciples, XDM prediction. The predicted values are compared with those calculated using free energy simulation in Table 1. Table 1: Dispersion component of side-chain hydration energies extrapolated from XDM and those calculated for each force field using REMD-TI. Energies are in kcal/mol. residue Ile Leu Ala Val Asn Tyr Thr Cys Ser Met Trp Hid Gln Phe Asp Lys Arg Glu NMA

XDM -10.94 -10.98 -4.52 -8.80 -8.28 -15.16 -7.54 -7.32 -5.39 -11.56 -18.45 -11.61 -10.39 -14.34 -8.51 -11.86 -14.09 -10.58 -10.37

CHARMM36 Amber ff14-SB OPLS-AA -11.00 -10.30 -10.28 -10.80 -10.00 -10.09 -4.46 -4.00 -4.15 -8.98 -8.40 -8.43 -9.47 -8.80 -9.18 -15.75 -13.90 -14.19 -8.01 -7.70 -7.63 -8.52 -6.40 -7.42 -5.94 -5.50 -5.51 -12.75 -10.50 -11.24 -19.10 -16.40 -16.81 -11.95 -11.70 -11.92 -11.24 -10.60 -10.93 -14.60 -12.60 -13.16 -8.42 -8.50 -8.89 -12.83 -11.60 -11.02 -16.75 -14.40 -12.50 -11.36 -10.40 -10.64 -12.36 -10.90 -11.51

Generally, the dispersion components of the hydration energies calculated using the free energy simulations are within 1 kcal/mol of the XDM-predicted values. There are some systematic deviations, such as the relatively small dispersion component for the aromatic residues (Phe, Trp, and Tyr) for the Amber ff14sb force field. For the CHARMM36 force field, the dispersion component of the hydration energies of sulfur-containing and nitrogencontaining residues show the greatest deviations, which are discussed in the following sections.

21


Page 22 of 43

Page 23 of 43

ΔG (kcal/mol)

10 CHARMM36 Amber ff14sb OPLS-AA exptl.

5

Cys

0 −5 −10

elec

disp

repul

total

15

ΔG (kcal/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Met

10 5 0 −5 −10 −15

elec

disp

repul

total

Figure 8: Components of the hydration energy of the Cys and Met side chain models. The CHARMM36 force field has the strongest dispersion interactions, but this is countered by a larger repulsive component.

22



3.4.2

Sulfur Containing Residues

The sulfur-containing side chains, Cys and Met, show a large variation in hydration energies (Figure 8). The difference in the dispersion component of the hydration energy is particularly significant; for the Cys side chain, this component ranges from -8.5 kcal/mol for the CHARMM36 force field to -6.4 kcal/mol for the Amber ff14sb force field. The OPLS-AA force field value of -7.4 kcal/mol is closest to the XDM-predicted value of -7.3 kcal/mol. This variation in dispersion energy stems from the atomic dispersion coefficients of the sulfur atoms (see Section 3.3); the Amber ff14sb sulfur C6 coefficients are less than half the OPLS-AA and CHARMM36 values. These disparate values stem from the variation in the Lennard-Jones parameters for sulfur atoms in the force fields (Table 2). 1.75 1.5 1.25

gS-O(r)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 43

1 0.75

AIMD CHARMM36 Amber OPLS

0.5 0.25 0 2.5

3

3.5

4

4.5

5

r (Å)

Figure 9: S–O radial distribution functions for methyl thiol (model for the Cys side chain) in an aqueous solution. The Amber ff14sb model underestimates the radius predicted by the AIMD simulation while the OPLS-AA and CHARMM36 models overestimate this difference. Table 2: Sulfur–oxygen radii (rS−O ) calculated from the methylthiol radial distribution functions, the force field Lennard-Jones parameters (i.e., σ and ε), and the dispersion component of the methylthiol solvation energy (∆Gdisp ). Method CHARMM36 Amber ff14sb OPLS-AA AIMD

rS−O (˚ A) σ (˚ A) 3.25 3.56 3.17 3.56 3.26 3.6 3.21

ε (kcal/mol) ∆Gdisp (kcal/mol) -0.45 -8.27 -0.25 -6.40 -0.425 -7.42

The hydration energy of the Cys side chain is underestimated by all three models, al23


Page 25 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


though there is a significant variation in the components. The XDM-predicted dispersion component of the hydration energy is approximately −8 kcal/mol, most similar to value predicted by the CHARMM36 force field. The repulsive energy predicted by the CHARMM36 and OPLS-AA models are higher than those of the Amber ff14sb model. Analysis of the S–O(H2 ) radial distribution function of aqueous methylthiol (Figure 9) shows that the CHARMM36 and OPLS-AA models yield an atomic radius (i.e., g(rS−O ) = 1) for sulfur that is approximately 0.05 ˚ A larger than the ab initio value, suggesting that the repulsive component is too large. The same trend is repeated for the methionine side chain. This analysis suggests that the parameters for the sulfur-containing residues could be improved in all three force fields.

3.5

Nitrogen-Containing Residues

For the CHARMM36 force field, residues containing one or more nitrogen atom (e.g., Lys, Arg, NMA, Asn, and Gln) tend to have larger water–side-chain dispersion interactions than predicted from the XDM results. For each of these molecules, the XDM-predicted dispersion component of the hydration energy is 1–3 kcal/mol smaller than the value calculated using the CHARMM36 free-energy simulations. The spuriously large contribution of nitrogen atoms to the non-polar component of the hydration energy has also been noted by Mobley et al. 78 This is consistent with the large atomic dispersion coefficients assigned to nitrogen atoms in the CHARMM36 force field (see Section 3.3), which results in a spurious increase in the strength of water–side-chain dispersion interactions. Arginine is an important example because the dispersion interactions between arginine and the interior of lipid bilayers are significant in the gating of ion channels 79,80 and arginine residues are often critical to binding nucleic acids. 81 As no experimental hydration energy of the arginine side chain is available, the XDM prediction is particularly informative. The dispersion component of the hydration energy calculated using the CHARMM36 force field is significantly larger than the XDM-predicted value (-16.75 kcal/mol vs -14.09 kcal/mol, 24



Page 26 of 43

respectively). The XDM analysis predicts that the Amber ff14-SB parameters are more consistent (-14.40 kcal/mol).

3.6

Improving the Description of Dispersion in Force Fields

The neglect of higher-order dispersion terms in the Lennard-Jones equation may explain why the force fields have high C6 dispersion coefficients but yield reasonably accurate hydration energies. Quantum chemical calculations have found that C8 /r8 and C10 /r10 terms account for approximately 30% of the dispersion interactions in crystals. 82 When these terms are neglected, the empirical parameterization process used to define Lennard-Jones parameters will increase the strength of the C6 /r6 interactions to compensate. To test this possibility, we calculated the XDM C6 , C8 , and C10 dispersion coefficients for methane and water. These data can be used to estimate the dispersion component of the hydration energy using the expression, ∆Gdisp = −4πρw

C6,w−mol C8,w−mol C10,w−mol + + , 3R3 5R5 7R7

(9)

where ρw is the number density of water, Cn,w−mol are the water–solute dispersion coefficients calculated using XDM, and R is the solute–water radius (taken from the arithmetic mean of methane and water Lennard-Jones radii, σ). The coefficients are calculated from the XDMcalculated molecular dispersion coefficients of the individual molecules, combined using the p relation Cn,w−mol = C8,w C8,w−mol . Based on this analytical expression, the C6 dispersion interaction accounts for the largest portion of the hydration energy, but the C8 and C10 terms contribute 25% and 12% of the energy, respectively (Table 3). The total interaction calculated from all the dispersion terms is −4.0 kcal/mol, which is close to the dispersion component of the Gibbs energy of hydration of methane (e.g., Ala) calculated with the CHARMM36 force field. Previous analysis by Floris et al. concluded that C8 and C10 interactions account for 12% of the

25


Page 27 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: The dispersion component of the hydration energy of methane calculated using the free-energy simulations (CHARMM36 alanine side-chain parameters) and those calculated using Eqn. 9 using the XDM C6 , C8 , and C10 coefficients. method n simulation analytical 6 8 10 total

∆G -4.4 -2.5 -1.0 -0.5 -4.0

dispersion component of the hydration energy of methane. 83 This supports the hypothesis that the dispersion interactions stemming from higher-order terms are included effectively in the C6 coefficients of the Lennard-Jones potential. There are several significant drawbacks to the practice of including higher-order dispersion interactions through large “effective” C6 coefficients. The distance dependence of this dispersion interaction will be spuriously strong at long range, but weak at short range. More generally, non-physical descriptions of dispersion interactions limit the accuracy of calculations where solutes move between environments, such as in membrane permeation and protein–ligand binding. Replacing the Lennard-Jones potential with a non-bonded potential that explicitly includes higher-order dispersion terms would resolve this issue in a more rigorous way. Nonbonded potentials that include higher-order dispersion terms have been used by chemical physicists for decades and there was discussion of including them in early biomolecular force fields. 84 For example, in 1938, Buckingham devised an equation of state for gaseous helium, neon, and argon using a potential that included C8 dispersion: 85

Enb (rij ) = A · exp (−b · r) −

C8,ij C6,ij − 8 . 6 rij rij

(10)

The primary challenge associated with the inclusion of higher order dispersion terms is that two to three dispersion coefficients would have to be defined for each atom type, although

26



methods like XDM could provide reasonable estimates. Separately, many of the anomalously large atomic dispersion coefficients are those in hydrogen-bonding groups (e.g., amides and hydroxyls). Adjustments to the Lennard-Jones parameters of these terms to capture the strength of hydrogen bonding despite the neglect of induced polarization and charge transfer may be the origin of this effect. The adoption of force fields that describe induced polarization explicitly, 86,87 or include explicit terms to represent hydrogen bonds, 88–91 may resolve this issue.

3.7

Evaluation of Nucleic Acid Force Fields

To test if trends in force-field dispersion coefficients identified for proteins also exist in nucleic-acid force fields, the dispersion coefficients of the nucleobases from the CHARMM36 92 and Amber OL15 93 force fields were compared to the XDM values (Figure 10). These models overestimate the molecular dispersion coefficients to an even larger degree, with dispersion coefficients that are more than 200% larger than the XDM values. A major cause of these large force field dispersion coefficients is the large atomic dispersion coefficient for the nucleobase nitrogen atoms; the CHARMM and Amber nucleic-acid force fields have dispersion coefficients of 74 a.u. and 58 a.u., respectively, for the nucleobase nitrogens, while XDM C6 coefficients for these atoms range from 13–19 a.u. The propensity of the force fields to attribute stronger dispersion interactions should favor states where the nucleobases are in close contact with each other (e.g., in the base-paired state vs a solvent-exposed state), although, in reality, base-pairing is predominantly an electrostatic interaction. 94 In the absence of experimental nucleobase hydration energies, we cannot perform the same analysis on the nucleic acid force fields as we did on the protein force fields. The relationships between the force field dispersion coefficients and the XDM dispersion coefficients are very similar to those of the protein force fields, so nucleic acid force fields are expected to demonstrate similar trends in the strength of dispersion interactions as the nitrogenous amino acids.

27


Page 28 of 43

Page 29 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


O

NH 2 N N H

N

N

NH

N H

N adenine

guanine O

NH2 N

O NH

NH

O N H cytosine

NH 2

N

N O H thymine

N O H uracil

Gua Ade Thy Cyt XDM CHARMM36 Amber

Ura 0

2000

4000 ∑ C6 (a.u.)

6000

Figure 10: Molecular dispersion coefficients for nucleobases calculated using XDM and the CHARMM36 and Amber OL15 force fields.

28



4

Conclusions

The XDM model was used to evaluate the dispersion coefficients in the CHARMM36, Amber ff14sb, and OPLS-AA force fields for proteins. In keeping with prior reports, we find that the C6 dispersion coefficients of these force fields are systematically larger than the QM-based XDM model. The Amber OL15 and CHARMM nucleic acid force fields also showed this trend. This trend results from large atomic dispersion coefficients, particularly for hydrogenbonding atom types. Interestingly, not all water models exhibit this trend; for example, the C6 dispersion coefficient of the popular TIP3P water model is close to the XDM value. Despite these large dispersion coefficients, hydration energies of the amino acid side chains calculated using the force fields are generally in good agreement with the experimental values. This is true even for aliphatic residues, where dispersion is the dominant attractive intermolecular interaction between the solvent and the solution. This suggests that the forcefield parameters capture the correct strength of water–side-chain dispersion interactions, even though the C6 dispersion coefficients are systematically larger than the XDM values. Internally consistent dispersion components of the solvation energy can also be estimated by this analysis, and the CHARMM force field was found to define anomalously strong dispersion interactions for the Cys, Met, Arg, and amide residues. Analysis of the methane–water solvation energy using an analytical expression shows that higher order terms (i.e., C8 and C10 ) account for a significant fraction of the dispersion energy. This suggests that the C6 coefficients in contemporary force fields are “effective” dispersion coefficients that capture neglected associative interactions, like higher-order dispersion, in addition to dispersion strictly due to instantaneous-dipole–induced-dipole interactions. Although this is clearly effective for describing some simulated properties of proteins, this practice may limit the accuracy of molecular simulations, particularly for processes where there is a net change in the total strength of dispersion interactions. Explicitly including higher-order dispersion terms in molecular-mechanical force fields is a potential solution to this problem. 29


Page 30 of 43

Page 31 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Acknowledgement The authors thank NSERC of Canada for funding through the Discovery Grants program (Application 418505-2012 and RGPIN-05795-2016). MM thanks the School of Graduate Studies at Memorial University for a graduate fellowship and Dr. Liqin Chen for a scholarship. EW thanks the Faculty of Science of Memorial University for a Science Undergraduate Research Award and the NSERC of Canada for an Undergraduate Student Research Award (USRA). Computational resources were provided by Compute Canada (RAPI: djk-615-ab). We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan Xp GPU used for this research.

30



Supporting Information Available Spreadsheets of the XDM and force field dispersion coefficients for each side chain model for each force field are available in the Supporting Information. This material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Stone, A. The Theory of Intermolecular Forces, 2nd ed.; Oxford University Press: Oxford, 2013. (2) Lazaridis, T.; Archontis, G.; Karplus, M. In Enthalpic Contribution to Protein Stability: Insights from Atom-Based Calculations and Statistical Mechanics; Anfinsen, C., Richards, F. M., Edsall, J. T., Eisenberg, D. S., Eds.; Advances in Protein Chemistry; Academic Press, 1995; Vol. 47; pp 231–306. (3) MacKerell, A. D.; Karplus, M. Importance of Attractive van der Waals Contribution in Empirical Energy Function Models for the Heat of Vaporization of Polar Liquids. J. Phys. Chem. 1991, 95, 10559–10560. (4) Moellmann, J.; Grimme, S. Importance of London Dispersion Effects for the Packing of Molecular Crystals: A Case Study for Intramolecular Stacking in a Bis-thiophene Derivative. Phys. Chem. Chem. Phys. 2010, 12, 8500–8504. (5) Fowkes, F. M. Additivity of Intermolecular Forces at Interfaces. I. Determination of the Contribution to Surface and Interfacial Tensions of Dispersion Forces in Various Liquids. J. Phys. Chem. 1963, 67, 2538–2541. (6) DiStasio, R. A.; von Lilienfeld, O. A.; Tkatchenko, A. Collective Many-Body van der Waals Interactions in Molecular Systems. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 14791–14795. 31


Page 32 of 43

Page 33 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(7) Jones, J. E. On the Determination of Molecular Fields. II. From the Equation of State of a Gas. Proc. R. Soc. London, A 1924, 106, 463–477. (8) Song, W.; Rossky, P. J.; Maroncelli, M. Modeling Alkane + Perfluoroalkane Interactions Using All-Atom Potentials: Failure of the Usual Combining Rules. J. Chem. Phys. 2003, 119, 9145–9162. (9) Baker, C. M.; Lopes, P. E. M.; Zhu, X.; Roux, B.; MacKerell Jr., A. D. Accurate Calculation of Hydration Free Energies Using Pair-Specific Lennard-Jones Parameters in the CHARMM Drude Polarizable Force Field. J. Chem. Theory Comput. 2010, 6, 1181–1198. (10) Riahi, S.; Rowley, C. N. Solvation of Hydrogen Sulfide in Liquid Water and at the Water–Vapor Interface Using a Polarizable Force Field. J. Phys. Chem. B 2014, 118, 1373–1380. (11) Berthelot, D. Sur Le Mélange Des Gaz. Compt. Rendus 1898, 126, 1703–1706. (12) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638–6646. (13) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225–11236. (14) Lamoureux, G.; MacKerell Jr., A. D.; Roux, B. A Simple Polarizable Model of Water Based on Classical Drude Oscillators. J. Chem. Phys. 2003, 119, 5185–5197. (15) Riahi, S.; Rowley, C. N. A Drude Polarizable Model for Liquid Hydrogen Sulfide. J. Phys. Chem. B 2013, 117, 5222–5229. (16) Adluri, A. N. S.; Murphy, J. N.; Tozer, T.; Rowley, C. N. Polarizable Force Field

32



with a σ-Hole for Liquid and Aqueous Bromomethane. J. Phys. Chem. B 2015, 119, 13422–13432. (17) Riniker, S. Fixed-Charge Atomistic Force Fields for Molecular Dynamics Simulations in the Condensed Phase: An Overview. J. Chem. Inf. Model. 2018, 58, 565–578. (18) Chen, I. J.; Yin, D.; MacKerell, A. D. Combined ab initio/empirical Approach for Optimization of Lennard-Jones Parameters for Polar-Neutral Compounds. J. Comput. Chem. 2002, 23, 199–213. (19) Adamovic, I.; Gordon, M. S. Dynamic Polarizability, Dispersion Coefficient C6 and Dispersion Energy in the Effective Fragment Potential Method. Mol. Phys. 2005, 103, 379–387. (20) Stone, A.; Misquitta, A. Atom–Atom Potentials from ab initio Calculations. Int. Rev. Phys. Chem. 2007, 26, 193–222. (21) Bajaj, P.; Götz, A. W.; Paesani, F. Toward Chemical Accuracy in the Description of Ion–Water Interactions through Many-Body Representations. I. Halide–Water Dimer Potential Energy Surfaces. J. Chem. Theory Comput. 2016, 12, 2698–2705. (22) Cole, D. J.; Vilseck, J. Z.; Tirado-Rives, J.; Payne, M. C.; Jorgensen, W. L. Biomolecular Force Field Parameterization via Atoms-in-Molecule Electron Density Partitioning. J. Chem. Theory Comput. 2016, 12, 2312–2323. (23) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696–3713. (24) Becke, A. D.; Johnson, E. R. Exchange-Hole Dipole Moment and the Dispersion Interaction Revisited. J. Chem. Phys. 2007, 127, 154108.

33


Page 34 of 43

Page 35 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(25) Otero-de-la Roza, A.; Johnson, E. R. Non-Covalent Interactions and Thermochemistry using XDM-Corrected Hybrid and Range-Separated Hybrid Density Functionals. J. Chem. Phys. 2013, 138, 204109. (26) Johnson, E. R. In Non-covalent Interactions in Quantum Chemistry and Physics; Otero-de-la Roza, A., DiLabio, G. A., Eds.; Elsevier, 2017; Chapter 5, pp 215–248. (27) Salem, L. The Calculation of Dispersion Forces. Mol. Phys. 1960, 3, 441–452. (28) Dalgarno, A.; Davison, W. D. The Calculation of van der Waals Interactions. Adv. At. Mol. Phys. 1966, 2, 1–32. (29) Otero-de-la Roza, A.; Johnson, E. R. Many-Body Dispersion Interactions from the Exchange-Hole Dipole Moment Model. J. Chem. Phys. 2013, 138, 054103. (30) Johnson, E. R. Dependence of Dispersion Coefficients on Atomic Environment. J. Chem. Phys. 2011, 135, 234109. (31) Christian, M. S.; Otero-de-la-Roza, A.; Johnson, E. R. Surface Adsorption from the Exchange-Hole Dipole Moment Dispersion Model. J. Chem. Theory Comput. 2016, 12, 3305–3315. (32) Gould, T.; Johnson, E. R.; Tawfik, S. A. Are Dispersion Corrections Accurate Outside Equilibrium? A Case Study on Benzene. 2018, 14, 1181–1191. (33) Dobson, J. Beyond Pairwise Additivity in London Dispersion Interactions. Int. J. Quantum Chem. 2014, 114, 1157–1161. (34) Becke, A. D.; Johnson, E. R. A Density-functional Model of the Dispersion Interaction. J. Chem. Phys. 2005, 123, 154101. (35) Mohebifar, M.; Johnson, E. R.; Rowley, C. N. Evaluating Force-Field London Dispersion Coefficients Using the Exchange-Hole Dipole Moment Model. J. Chem. Theory Comput. 2017, 13, 6146–6157. 34



(36) Huang, J.; MacKerell, A. D. CHARMM36 All-Atom Additive Protein Force Field: Validation Based on Comparison to NMR Data. J. Comput. Chem. 2013, 34, 2135– 2145. (37) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides. J. Phys. Chem. B 2001, 105, 6474–6487. (38) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian 09 Revision D.01. Gaussian Inc. Wallingford CT 2009. (39) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods Without Adjustable Parameters: The PBE0 Model. J. Chem. Phys. 1999, 110, 6158–6170. (40) Woon, D. E.; Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. III. The Atoms Aluminum Through Argon. J. Chem. Phys. 1993, 98, 1358–1371. (41) Hickey, A. L.; Rowley, C. N. Benchmarking Quantum Chemical Methods for the Calculation of Molecular Dipole Moments and Polarizabilities. J. Phys. Chem. A 2014, 118, 3678–3687. (42) PyPostg. https://github.com/RowleyGroup/postg-parser, (accessed March 17, 2017). (43) Postg. http://schooner.chem.dal.ca, (accessed June 5, 2018). (44) Deng, Y.; Roux, B. Hydration of Amino Acid Side Chains: Nonpolar and Electrostatic Contributions Calculated from Staged Molecular Dynamics Free Energy Simulations with Explicit Water Molecules. J. Phys. Chem. B 2004, 108, 16567–16576.

35


Page 36 of 43

Page 37 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(45) Jiang, W.; Roux, B. Free Energy Perturbation Hamiltonian Replica-Exchange Molecular Dynamics (FEP/H-REMD) for Absolute Ligand Binding Free Energy Calculations. J. Chem. Theory Comput. 2010, 6, 2559–2565. (46) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926–935. (47) Lin, Y.-L.; Aleksandrov, A.; Simonson, T.; Roux, B. An Overview of Electrostatic Free Energy Computations for Solutions and Proteins. J. Chem. Theory Comput. 2014, 10, 2690–2709. (48) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781–1802. (49) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I. 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, 671–690. (50) Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: A Web-Based Graphical User Interface for CHARMM. J. Comput. Chem. 2008, 29, 1859–1865. (51) Buckle, A. M.; Schreiber, G.; Fersht, A. R. Protein-Protein Recognition: Crystal Structural Analysis of a Barnase-Barstar Complex at 2.0 ˚ A Resolution. Biochemistry 1994, 33, 8878–8889. (52) Wacker, D.; Wang, S.; McCorvy, J. D.; Betz, R. M.; Venkatakrishnan, A.; Levit, A.; Lansu, K.; Schools, Z. L.; Che, T.; Nichols, D. E. et al. Crystal Structure of an LSDBound Human Serotonin Receptor. Cell 2017, 168, 377 – 389.e12.

36



(53) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. Comp. Phys. Comm. 2005, 167, 103–128. (54) Zhang, Y.; Yang, W. Comment on “Generalized Gradient Approximation Made Simple”. Phys. Rev. Lett. 1998, 80, 890. (55) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys 2010, 132, 154104. (56) Bankura, A.; Karmakar, A.; Carnevale, V.; Chandra, A.; Klein, M. L. Structure, Dynamics, and Spectral Diffusion of Water from First-Principles Molecular Dynamics. J. Phys. Chem. C 2014, 118, 29401–29411. (57) Goedecker, S.; Teter, M.; Hutter, J. Separable Dual-Space Gaussian Pseudopotentials. Phys. Rev. B 1996, 54, 1703–1710. (58) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327–341. (59) Nicholls, A.; Sharp, K. A.; Honig, B. Protein Folding and Association: Insights from the Interfacial and Thermodynamic Properties of Hydrocarbons. Proteins 1991, 11, 281–296. (60) Levy, R. M.; Zhang, L. Y.; Gallicchio, E.; Felts, A. K. On the Nonpolar Hydration Free Energy of Proteins: Surface Area and Continuum Solvent Models for the Solute–Solvent Interaction Energy. J. Am. Chem. Soc. 2003, 125, 9523–9530. (61) Dzubiella, J.; Swanson, J. M. J.; McCammon, J. A. Coupling Hydrophobicity, Dis-

37


Page 38 of 43

Page 39 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


persion, and Electrostatics in Continuum Solvent Models. Phys. Rev. Lett. 2006, 96, 087802. (62) Tan, C.; Tan, Y.-H.; Luo, R. Implicit Nonpolar Solvent Models. J. Phys. Chem. B 2007, 111, 12263–12274. (63) Fennell, C. J.; Kehoe, C.; Dill, K. A. Oil/Water Transfer Is Partly Driven by Molecular Shape, Not Just Size. J. Am. Chem. Soc. 2010, 132, 234–240. (64) Wang, J.; Deng, Y.; Roux, B. Absolute Binding Free Energy Calculations Using Molecular Dynamics Simulations with Restraining Potentials. Biophys. J. 2006, 91, 2798– 2814. (65) Lin, Y.-L.; Meng, Y.; Jiang, W.; Roux, B. Explaining Why Gleevec Is a Specific and Potent Inhibitor of Abl Kinase. Proc. Natl. Acad. Sci. U.S.A. 2013, 110, 1664–1669. (66) Riahi, S.; Rowley, C. N. Why Can Hydrogen Sulfide Permeate Cell Membranes? J. Am. Chem. Soc. 2014, 136, 15111–15113. (67) Awoonor-Williams, E.; Rowley, C. N. Molecular Simulation of Nonfacilitated Membrane Permeation. Biochim. Biophys. Acta – Biomembranes 2016, 1858, 1672–1687. (68) Zeiss, G.; Meath, W. J. The H2O-H2O Dispersion Energy Constant and the Dispersion of the Specific Refractivity of Dilute Water Vapour. Mol. Phys. 1975, 30, 161–169. (69) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (70) Bultinck, P.; Alsenoy, C. V.; Ayers, P. W.; Carbó-Dorca, R. Critical Analysis and Extension of the Hirshfeld Atoms in Molecules. J. Chem. Phys. 2007, 126, 144111. (71) Horta, B. A. C.; Merz, P. T.; Fuchs, P. F. J.; Dolenc, J.; Riniker, S.; Hnenberger, P. H. A GROMOS-Compatible Force Field for Small Organic Molecules in the Condensed Phase: The 2016H66 Parameter Set. J. Chem. Theory Comput. 2016, 12, 3825–3850. 38



(72) Villa, A.; Mark, A. E. Calculation of the Free Energy of Solvation for Neutral Analogs of Amino Acid Side Chains. J. Comput. Chem. 2002, 23, 548–553. (73) MacCallum, J. L.; Tieleman, D. P. Calculation of the Water–Cyclohexane Transfer Free Energies of Neutral Amino Acid Side-Chain Analogs Using the OPLS All-Atom Force Field. J. Comput. Chem. 2003, 24, 1930–1935. (74) Hess, B.; van der Vegt, N. F. A. Hydration Thermodynamic Properties of Amino Acid Analogues: A Systematic Comparison of Biomolecular Force Fields and Water Models. J. Phys. Chem. B 2006, 110, 17616–17626. (75) Chang, J.; Lenhoff, A. M.; Sandler, S. I. Solvation Free Energy of Amino Acids and Side-Chain Analogues. J. Phys. Chem. B 2007, 111, 2098–2106. (76) Zhang, H.; Yin, C.; Jiang, Y.; van der Spoel, D. Force Field Benchmark of Amino Acids: I. Hydration and Diffusion in Different Water Models. J. Chem. Inf. Model. 2018, 58, 1037–1052. (77) Labute, P. The Generalized Born/Volume Integral Implicit Solvent Model: Estimation of the Free Energy of Hydration Using London Dispersion Instead of Atomic Surface Area. J. Comput. Chem. 2008, 29, 1693–1698. (78) Mobley, D. L.; Bayly, C. I.; Cooper, M. D.; Shirts, M. R.; Dill, K. A. Small Molecule Hydration Free Energies in Explicit Solvent: An Extensive Test of Fixed-Charge Atomistic Simulations. J. Chem. Theory Comput. 2009, 5, 350–358. (79) Long, S. B.; Tao, X.; Campbell, E. B.; MacKinnon, R. Atomic Structure of a VoltageDependent K+ Channel in a Lipid Membrane-like Environment. Nature 2007, 450, 376. (80) Schwaiger, C. S.; Börjesson, S. I.; Hess, B.; Wallner, B.; Elinder, F.; Lindahl, E. The

39


Page 40 of 43

Page 41 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Free Energy Barrier for Arginine Gating Charge Translation Is Altered by Mutations in the Voltage Sensor Domain. PLOS ONE 2012, 7, 1–11. (81) Weiss, M. A.; Narayana, N. RNA Recognition by Arginine-rich Peptide Motifs. Biopolymers 1998, 48, 167–180. (82) Otero-de-la Roza, A.; Johnson, E. R. A Benchmark for Non-Covalent Interactions in Solids. J. Chem. Phys. 2012, 137, 054103. (83) Floris, F. M.; Tomasi, J.; Ahuir, J. L. P. Dispersion and Repulsion Contributions to the Solvation Energy: Refinements to a Simple Computational Model in the Continuum Approximation. J. Comput. Chem. 1991, 12, 784–791. (84) Pastor, R. W. Topics in Stochastic Dynamics of Polymers. Ph.D. thesis, Harvard University, 1984. (85) Buckingham, R. A. The Classical Equation of State of Gaseous Helium, Neon and Argon. Proc. R. Soc. Lond. A 1938, 168, 264–283. (86) Shi, Y.; Xia, Z.; Zhang, J.; Best, R.; Wu, C.; Ponder, J. W.; Ren, P. Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins. J. Chem. Theory Comput. 2013, 9, 4046–4063. (87) Lopes, P. E. M.; Huang, J.; Shim, J.; Luo, Y.; Li, H.; Roux, B.; MacKerell Jr., A. D. Polarizable Force Field for Peptides and Proteins Based on the Classical Drude Oscillator. J. Chem. Theory Comput. 2013, 9, 5430–5449. (88) Morozov, A. V.; Kortemme, T.; Tsemekhman, K.; Baker, D. Close Agreement Between the Orientation Dependence of Hydrogen Bonds Observed in Protein Structures and Quantum Mechanical Calculations. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 6946– 6951.

40



(89) Huang, J.; Meuwly, M. Explicit Hydrogen-Bond Potentials and Their Application to NMR Scalar Couplings in Proteins. J. Chem. Theory Comput. 2010, 6, 467–476. (90) Choi, H.; Kang, H.; Park, H. New Angle-Dependent Potential Energy Function for Backbone–Backbone Hydrogen Bond in Protein–Protein Interactions. J. Comput. Chem. 2010, 31, 897–903. (91) Korth, M. Empirical Hydrogen-Bond Potential Functions – An Old Hat Reconditioned. ChemPhysChem 2011, 12, 3131–3142. (92) Hart, K.; Foloppe, N.; Baker, C. M.; Denning, E. J.; Nilsson, L.; MacKerell, A. D. Optimization of the CHARMM Additive Force Field for DNA: Improved Treatment of the BI/BII Conformational Equilibrium. J. Chem. Theory Comput. 2012, 8, 348–362. ˇ (93) Zgarbová, M.; Sponer, J.; Otyepka, M.; Cheatham, T. E.; Galindo-Murillo, R.; Jureˇcka, P. Refinement of the Sugar–Phosphate Backbone Torsion Beta for AMBER Force Fields Improves the Description of Z- and B-DNA. J. Chem. Theory Comput. 2015, 11, 5723–5736. (94) Su, P.; Li, H. Energy Decomposition Analysis of Covalent Bonds and Intermolecular Interactions. J. Chem. Phys. 2009, 131, 014102.

41


Page 42 of 43

Page 43 of 43

Graphical TOC Entry

6000 4000

CHARMM Amber OPLS-AA

2000 0 0

2000 4000 ∑ C6 XDM (a.u.)

...but predict hydration energies accurately ΔGcalc (kcal/mol)

Dispersion Coefficients in Protein Force Fields are ~50% Higher than QM-XDM ∑ C6 FF (a.u.)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5 0 −5

CHARMM Amber OPLS-AA

−10 −10 −5 0 ΔGexptl (kcal/mol)

5

Neglected Higher-Order Dispersion is Suspected ( )= −

6 6

−

8 8

−

10 10

but only

42

( )=


12

−

6 6

Evaluating the London Dispersion Coefficients of Protein Force Fields

Recommend Documents