A physically grounded damped dispersion model with particle mesh Ewald summation Joshua A. Rackers, Chengwen Liu, Pengyu Ren, and Jay W. Ponder Citation: The Journal of Chemical Physics 149, 084115 (2018); doi: 10.1063/1.5030434 View online: https://doi.Org/10.1063/1.5030434 View Table of Contents: http://aip.scitation.Org/toc/jcp/149/8 Published by the American Institute of Physics Articles you may be interested in Toward quantum-chemical method development for arbitrary basis functions The Journal of Chemical Physics 149, 084106 (2018); 10.1063/1.5044765 Fitting a round peg into a round hole: Asymptotically correcting the generalized gradient approximation for correlation The Journal of Chemical Physics 149, 084116 (2018); 10.1063/1.5021597 Single-root networks for describing the potential energy surface of Lennard-Jones clusters The Journal of Chemical Physics 149, 084102 (2018); 10.1063/1.5043330 The pressure in interfaces having cylindrical geometry The Journal of Chemical Physics 149, 084109 (2018); 10.1063/1.5037054 A novel density functional theory for atoms, molecules, and solids The Journal of Chemical Physics 149, 074104 (2018); 10.1063/1.5038262 Hybrid grid/basis set discretizations of the Schrodinger equation The Journal of Chemical Physics 147, 244102 (2017); 10.1063/1.5007066 Take a closer look at what these 11 <• 1 ,, ,, . PRESENTED BY environmentally friendly adhesive systems can do &MASTERBOND" PHYSICS TODAY WHITEPAPERS THE JOURNAL OF CHEMICAL PHYSICS 149, 084115 (2018) A physically grounded damped dispersion model with particle mesh Ewald summation Joshua A. Rackers,1 Chengwen Liu,2 Pengyu Ren,2 and Jay W. Ponder1'3'a) 1 Program in Computational and Molecular Biophysics, Washington University School of Medicine, Saint Louis, Missouri 63110, USA ^Department of Biomedical Engineering, The University of Texas at Austin, Austin, Texas 78712, USA 3Department of Chemistry, Washington University in Saint Louis, Saint Louis, Missouri 63130, USA (Received 21 March 2018; accepted 16 August 2018; published online 31 August 2018) Accurate modeling of dispersion is critical to the goal of predictive biomolecular simulations. To achieve this accuracy, a model must be able to correctly capture both the short-range and asymptotic behavior of dispersion interactions. We present here a damped dispersion model based on the overlap of charge densities that correctly captures both regimes. The overlap damped dispersion model represents a classical physical interpretation of dispersion: the interaction between the instantaneous induced dipoles of two distinct charge distributions. This model is shown to be an excellent fit with symmetry adapted perturbation theory dispersion energy calculations, yielding an RMS error on the S 101x7 database of 0.5 kcal/mol. Moreover, the damping function used in this model is wholly derived and parameterized from the electrostatic dipole-dipole interaction, making it not only physically grounded but transferable as well. Published by AIP Publishing. https://doi.org/10.1063/L5030434 I. INTRODUCTION The range of possible problems for molecular mechanics models to solve is immense. For problems that are too large to solve with Schrodinger's equation but too small to be observed experimentally, we rely on classical models to make predictions and generate hypotheses. This ability has made molecular mechanics force fields integral to the study of problems from RNA folding1 to new alloy characterization.2 Because they are classical approximations to quantum mechanical reality, the success of these models is entirely dependent on how accurate that approximation is on a wide variety of systems. To achieve this, most force fields split the interaction energies of interacting atoms into physically meaningful components. Among the most significant of these components is the dispersion interaction that arises from the correlation of instantaneous induced dipoles. No force field can provide fully accurate predictions for every component of the total energy of a system. In current models, this has been typically handled by careful cancellation of errors between the various components (electrostatics, polarization, repulsion, dispersion, etc.). More recently, however, a new crop of next-generation force fields is emerging that aim to reduce this dependence on error cancellation by comparing directly to ab initio energy decomposition analysis (EDA) data.3"10 We are working on a model with this same objective. Previously we have shown that it is possible to accurately model electrostatics (to within 1 kcal/mol) in regions where previously error cancellation had long been relied upon, the so-called "charge penetration" error.11 In a'Author to whom correspondence should be addressed: ponder @ dasher. wustl.edu this work, we shall demonstrate that the same is possible for the dispersion interaction component. While this does not represent a complete force field capable of condensed phase simulations, it is an important step toward such a full model. Accurately modeling dispersion in classical force fields is known to be important, particularly for biological systems. On a phenomenological level, dispersion is what causes neutral atoms and molecules to be weakly attracted to each other. This makes it essential to modeling simple Lennard-Jones fluids such as liquid argon, but it is also critically important to more complex systems. Dispersion has been shown to be an essential component of modeling nucleic acid structure,12 where it contributes to the so-called stacking energy of nucleic acid bases. It is known to play a part in halogen bonding, supporting, along with electrostatics, the stabilization energy of the interaction.13'14 Additionally, long-range dispersion is widely recognized to be important for the simulation of lipid bilayers.15 This broad spectrum of applications motivates the necessity of accurate dispersion models. The history of dispersion models dates back to Fritz London, who first established the canonical 1/r6 dependence of the London dispersion energy. This model has been enormously influential. The vast majority of biological force fields in use today still use this simple model (Amber,16 CHARMM,17 etc.) or derivatives thereof such as the attractive part of Halgren's buffered 14-7 potential18 used in the AMOEBA19 force field. It is well known, however, that the l/rn potential expansion breaks down for short-range interactions where charge distributions of interacting molecules overlap.20 There is a long history of attempts to correct this divergence though the use of damping functions. An important early damped dispersion 0021 -9606/2018/149(8)/084115/18/$30.00 149, 084115-1 Published by AIP Publishing. 084115-2 Rackers et al. J. Chem. Phys. 149, 084115 (2018) model was the empirical HFD (Hartree-Fock-Dispersion) scheme proposed by Scoles and co-workers.21 Another notable attempt to describe this phenomenon was undertaken by Tang and Toennies who introduced a damping function parameterized to account for the overlap in charge distributions. A comprehensive review of dispersion damping functions is beyond the scope of this work, but the original Tang and Toennies report provides a thorough overview of dispersion damping functions up to that point. These types of formalisms have seen the widest use as dispersion corrections to density functional theory (DFT) calculations.23 DFT-D schemes have used the Tang-Toennies function, as well as various other damping functions proposed by Wu and Yang,24 Chai and Head-Gordon,25 and Johnson and Becke.26 Despite wide use in the DFT community, damped dispersion functions have been taken up in decidedly fewer molecular mechanics models. Notably, the Effective Fragment Potential (EFP) model employs a dispersion model that utilizes an overlay-based parameter-free modification of the Tang-Toennies damping function.27'28 And recently, Verma et al. proposed using the dispersion part of the DFT-D3 formulation of Grimme as a molecular mechanics model. However, while it has been shown that previous damping functions can effectively account for the change in dispersion upon charge overlap, they do so largely empirically. In the case of the Tang-Toennies damping function, for example, the form is based on a Born-Mayer potential described by an empirically fit width parameter. In this paper, we propose a damped dispersion function similar in spirit to that of Tang and Toennies but rooted in a physical model of charge distribution overlap. In previous work, we have shown that a relatively simple model can capture the physical extent of atomic charge distributions that leads to the so-called charge penetration error in electrostatic interactions between molecules.11 Here we will show that this same model can be used directly and without modification to create a dispersion model that is elegantly unified with the electrostatic model. This unification is possible because both the electrostatic and dispersion terms depend on the density. The electrostatic term is simply the interaction between two static densities, while the dispersion term arises from the interaction of densities associated with instantaneous induced dipoles. In this work, we will show that the same rough description of the density can be used in both cases to great effect. This will be done in five parts. First, we elucidate the theory that starts from dipole-dipole interactions and gives rise to this new damped dispersion function. Second, we describe the methods of the study. Third, we evaluate the performance of this function against benchmark Symmetry Adapted Perturbation Theory (SAPT) calculations. Fourth, we will describe how the model has been implemented with dispersion particle mesh Ewald (DPME) to boost its efficiency. And finally, we will discuss the implications of this work and some general conclusions. II. THEORY To present our new damped dispersion model, we shall first revisit a simple derivation of the original London dispersion model. We do so first and foremost because it forms the basis for our damped model, but also because it is instructive. One of the defining characteristics of a damped dispersion model, as we shall argue later in the paper, is that it has a straightforward physical interpretation. Dispersion is correctly said to be a fundamentally non-classical phenomenon, but the model we use to describe it need not to be so bound. We will show that an interpretable model of dispersion can be constructed from physical models of atomic polarizability and charge density. A. London dispersion For our description of canonical London dispersion energy, we will follow that of Maitland et al.30 The dispersion energy between two atoms arises from the interaction between instantaneous dipoles of these atoms. To model this system, we consider a simplified one-dimensional Drude oscillator model, as illustrated in Fig. 1. In this representation, each atom is represented by a fixed charge +Q bound by a spring with spring constant, k, and an equal and opposite charge -Q with mass, M. This model is crude, but it captures the essential elements of the dispersion interaction. At any point in time, each atom has a dipole moment, [i = Qz (dependent on the atomic polarizability determined by k) and those dipole moments are free to interact with each other. When atom i and atom j are infinitely separated, the Schrodinger equation for each can be written as 1 d2x¥; m dz2 + n^Ei l. fcf m: = o, (1) where the potential energy term is merely the energy of a simple harmonic oscillator. The same can be written for atom j. The solutions to this equation can be found trivially, yielding ground state energies of E{ = —fujQ and Ej = —fujQ, where the frequency, a>o, is (2) (3) In the complete non-interacting limit, the total energy of the system is E(r -> oo) = Et + Ej = hcoo- (4) This limit in itself is not useful, but if we consider what happens when the two atoms get closer, we shall see that it sets a useful reference for our potential energy function. If we bring the two atoms closer so that they do interact, but not so close that their charge distributions overlap, our Schrodinger equation is no longer trivially separable. The wave equation for two -Q +Q •^AAAAAVO -Q +Q FIG. 1. Classical model of dispersion. 084115-3 Rackers et al. J. Chem. Phys. 149, 084115 (2018) interacting atoms now includes the electrostatic interaction between the two dipoles and becomes 1 d2*¥ 1 d2*¥ 2 I 1 2 1 2 \ _ M~dJ+M~dJ+¥\E ~ 2kZi ~ 2kzj ~ u'<"*™***]* " °" ' j (5) One can see that in addition to the simple harmonic oscillator terms, a new potential appears in Eq. (5). This is the potential energy at any given instant between the two interacting instantaneous multipole distributions. For the dipole-dipole interaction of the simple Drude model of Fig. 1, the form of this potential is easily obtained from simple electrostatics Uelectrostatic ~ ^dipole-dipole ~ VV[/( chg—chg ßi ■ fij - 3 (6) If we plug in the Drude dipoles from Fig. I, fi = Qz, Eq. (6) becomes 1 ^dipole-dipole ~ a (7) This dipole-dipole energy is the source, as we shall show, of the canonical 1/r6 leading term dependence of the dispersion energy. Combining Eq. (7) with Eq. (5) yields 1 d2y¥ 1 d2y¥ 2 , ___I____I__\£. M dz2 M dz2 n2' 1. -M .2 ¥ = 0. (8) Following the transformation of variables of Maitland, Rigby, Smith, and Wakeham, we define Zi+Zj V2 Ä2 V2 (9) and rewrite Eq. (8) as 1 d2x¥ 1 <92vP 2v M dz2 where M dz2 h2 2Q2 2Q -3 i2 (ID k\ — k H--—, k2 — k - Equation (10) is simply a transformed version of the original problem of two independent harmonic oscillators. It can be solved in the same manner giving E(r) = ^h(ü>i +w2), (12) where h I 2Q2 Wi — a/ — = Woa/ 1--;—, Cl>2 — a/ — — VAf V r3k y M 1 + 2ß2 r3k (13) One can see that as r becomes large, ojj and a>2 converge to a>o where we recover the independent oscillator solution. For small perturbations, a>] and a>2 can be approximated with a binomial expansion Vl +x = 1 + l/2x- l/8x+ ■ so the total energy becomes E(r) — hdjQ 2r6k2 (14) (15) The final step is to subtract the energy of infinitely separated atoms. This gives the dispersion potential energy QAhu0 Udispersion = E(r) - £(oo) = - + • • • , (16) where the canonical r6 dependence arises from the first nonzero term from the application of the binomial expansion. It should be noted that while the dipole-dipole interaction is the dominant electrostatic term of Eq. (5), there are terms arising from higher-order multipole interactions as well. The dipole-quadrupole and quadrupole-quadrupole interactions giving rise to the 1/r8 and 1/r10 potentials are derived in the Appendix. There are a number of models that use these terms, including EFP, SIBFA (Sum of Interactions Between Fragments Ab initio computed), and Misquitta and Stone's model for small organic molecules.27'31'32 For a perspective on the importance of these higher order terms for the case of the neon dimer, the reader is directed to the work of Bytautas and Ruedenberg.33 The latter reference showed that even for this simple dimer, the 1/r8 and 1/r10 terms are nearly impossible to distinguish at reasonable separations. There are odd-power terms (1/r7, 1/r9, etc.) that can be included in the expansion as well. These arise from the mixing of the even order terms, are highly angularly dependent, and spherically average to zero at long range.34 There has also been recent work on incorporating these terms into dispersion models,35'36 where these higher-order terms give successively better approximations to the exact dispersion energy. As we will show, however, to reach the stated accuracy goal of < 1 kcal/mol, only the leading term will be necessary. For most systems, the perturbation of the dipole-dipole interaction energy is small compared to the energy holding the electrons to their respective atoms. This makes taking only the leading r6 term of Eq. (16) a good approximation for most long-range intermolecular interactions. In practice, this is done by introducing a parameter, C, to capture this dependence 084115-4 Rackers eř a/. J. Chem. Phys. 149, 084115 (2018) jjdispen C6C6 (17) This model will be referred to throughout the remainder of the paper as the London dispersion model. Unfortunately, this method of approximation starts to break down when the charge distributions of interacting atoms start to overlap. We will handle this situation through the introduction of short-range damping, but rather than relying on empiricism for the damping function, we look to the underlying electrostatics to provide a consistent model. B. Short-range electrostatics A long-standing problem in the modeling of electrostatics for molecular mechanics models is the so-called charge penetration error. The error arises when charge distributions of interacting atoms overlap, causing the true electrostatic energy of the interacting densities to diverge from the point charge or point multipole approximation. We have shown in previously published studies11'37'38 that a simple hydrogen-like approximation of the Coulomb potential does a remarkably good job at correcting this error. Why is this germane to a study of dispersion? Dispersion, as shown above, can be modeled as arising from a dipole-dipole interaction. In the context of the multipolar AMOEBA force field, we have shown that the hydrogen-like approximation to the Coulomb interaction can be extended to the interactions between higher-order multipole moments. In fact, including these corrections for charge-dipole, dipole-dipole, dipole-quadrupole, etc. interactions is essential to the transferability and accuracy of the model.11 Here we show that the dipole-dipole interaction arising from this earlier model can be used directly to create a new damped dispersion model. To illustrate where the dipole-dipole damping comes from, we follow a similar derivation to that of Ref. 11. The potential due to the electrons for this model is defined as V(r)=%-(l-e-a'r), (18) r where r is the distance from the center of the charge distribution and a is a parameter describing the width of the distribution. Application of Poisson's equation, V2V = A so yields the corresponding density, p(r) = -l—e a'r. (20) These two quantities can be used to approximate the Coulomb interaction energy between two charge distributions, PiiXdpjiTj) U' chg—chg electrostatic _ J J Pi.(r^(r,) dYidY] (21) Application of the one-center integral method of Coulson39 gives U chg—chg electrostatic mi r a2 j a2 \ (22) Equation (22) gives the charge-charge electrostatic energy. To get the dipole-dipole energy, recall that the full multipole energy of the i-j interaction can be written as jjtotal _ jjchg—chg _j_ jjchg—dipole _j_ jjdipole—chg electrostatic _l_ jj-dipole-dipole _|_ , , , = 'A'/;/'// + 'l^'l'iifij ~ fi^'l'a'lj + HiVVTijUj + ■ ■ ■ . (23) For a point-point interaction, Ty is simply 1/r, but for our model, direct inspection of Eq. (22) yields la2 „2 \ j 1 Tij = -r 1 rdamp j] (24) We can now apply this new relation for Tjj to the definition of the dipole-dipole energy from Eq. (23), Udamp~diP°le = ^'VVr9^ _ ,damp fii ■ fij _ fdamp 3(^' ' (r'J ' M/) where /j and/j are the damping terms that come from the damp ■ (19) derivatives of rhe//^ term of Eq. (24) rdamp _ i h - 1 ~ ■damp --J---(1 + air)e-a'r ---í-- (a2 - a2j (a2 - a2j a2 j W - °f) 1 +atr+ -(air)2)e a c c ra 10 5 0 AMBER ■ CHARMM London Dispersion ■ Overlap Damped Dispersion FIG. 10. Mean unsigned error in dispersion energy for nucleic acid structures. Model error is relative to SAPT for each of the six structural parameters. The overlap damped dispersion model reduces the error in the dispersion across all six degrees of freedom. Amber and CHARMM results from Ref. 49. Fig. 10 that parameterizing a 1/r6 (London dispersion) model directly to SAPT reduces some of the need for the cancellation of error, but not all. The second and more important feature one observes is the agreement throughout the potential energy surface of the overlap damped dispersion model. In addition to relieving itself of the cancellation of errors burden, one can see that the damped model provides a minimum factor of two improvements in the mean unsigned error over the undamped London dispersion model for every degree of freedom. This has little to do with the behavior of the dispersion energy at equilibrium; the divergence occurs primarily for structures where the electron densities of the two base-pairs start to overlap. As an instructive example, take the change in dispersion energy with respect to the tilt angle for the CATG base step shown in Fig. 11. One can see that at equilibrium, both the London and overlap damped dispersion models predict the SAPT dispersion energy with good precision. However, as one changes the tilt angle in either direction, the dispersion energy of the undamped model diverges quickly from the SAPT while the overlap damped dispersion model follows the shape of the SAPT curve with fidelity. This trend holds across all six degrees of freedom and all ten base pair steps. Plots like Fig. for each combination are available in the supplementary material. The divergence observed for non-damped models matters because it is not simply confined to high total energy areas of the DNA potential energy surface. In fact, Parker and Sher-rill showed that for the stacked A-C pair (one half of the CATG base step) at a tilt angle of -15°, the total energy is -5 kcal/mol. This is only 0.5 kcal/mol above the minimum total energy of -5.5 kcal/mol. Figure 11 suggests that in order to accurately model this region of the potential energy surface without large cancellation of errors, a damped dispersion model is necessary. SAPT0 Dispersion Overlap Damped Dispersion London Dispersion Amber Dispersion CHARMM Dispersion -5 0 5 Tilt (degrees) 10 15 FIG. 11. Dispersion energy of CATG interaction vs. tilt. The non-damped dispersion models uniformly overestimate the magnitude of the dispersion energy as the angle varies from equilibrium in either direction. The overlap damped dispersion model predicts the shape of the SAPT curve at both equilibrium and near-equilibrium geometries. V. DISPERSION PARTICLE MESH EWALD SUMMATION Accuracy and efficiency are both important features of a molecular mechanics model. A good dispersion model must not only be accurate, but also fast to compute. While the accuracy of the overlap damped dispersion model has been solidly established in this paper, the exponentials required for its evaluation have the potential to slow potential energy calculations. To make the overlap damped dispersion model computationally efficient and tractable for use in biomolecular simulations, we have implemented the model with particle mesh Ewald (PME) summation in the Tinker molecular mechanics software package. In this section, we present a brief overview of the damped dispersion PME implementation and show how this implementation provides a substantial speed and accuracy improvement over the standard cutoff-based van der Waals implementation. Ewald summation is classically considered to be primarily a solution to the pairwise long-range electrostatics problem. The 21/r electrostatic potential is conditionally 084115-13 Rackers et al. J. Chem. Phys. 149, 084115 (2018) convergent which makes direct computation of the electrostatic energy of a periodic system difficult. To circumvent this problem, Ewald methods split the sum into short-range and long-range parts, with short-range part being computed directly and the long-range via Fourier transformation. This separation not only makes periodic calculations possible, but also increases the speed with which the energy and gradient can be evaluated. The same method can be applied to the dispersion energy calculation. Here we note that the following derivation is by no means original. In fact, Essman and co-workers proposed the possibility of using particle mesh Ewald summation for dispersion in their 1995 paper describing the method of smooth particle mesh Ewald summation.59 We present here a brief summary simply to show that the inclusion of a damping term in this case does not change the ability to use the method. The total dispersion energy, as given by Eq. (33) is cl d 9 jjdamp _ \ 1 6 6 / rdamp \z (36) dispersion / i ^.6 Vdispersion)jj' ^ ' i*j ij This can be split into a short-range part, a long-range part, and a "self term, jjdispersion _ ^dispersion ^ ^dispersion ^ ^dispersion (37) short—range long—range self ' ^ ^ total with ^dispersion _ 6 6 If damp ) M + fl2r? + _fi4r4. le"^ short—range / i 6 Vdispersion) (A ^ ij 9^ ij J u: r.. 9 9/2 dispersion _ ^-^T ^ ^ «2 2 long—range 3 y 111*0 U-M ~ 12 ^ { 6V 1^ ! ---J\ - 2{n\m\l ßf)e^m^)2^e^m^ S(m)S(-m), fl6 ß37T3/2 ^dispersion _ P \ 1 ^2 _j_ P self (38a) (38b) (38c) Equations (38a) and (38b) are commonly known as the direct space sum and reciprocal space sum, respectively. The variable, /3, is the parameter determining the Gaussian width, m is defined by the reciprocal lattice vectors, a, as m = mjai* + rri2a2* + 111333*, and V is the volume of the unit cell. The structure factor, S, is defined for dispersion as f2n(m-Tj) (39) One can see that for cutoffs longer than 6 A, the energy of the PME implementation is constant due to the fact that/jam;, is effectively unity for all atom pairs outside of this radius. For our model, we chose this cutoff of 6 A to balance the direct and reciprocal space computational effort. Comparing the PME and non-PME curves in Fig. 12 illustrates an obvious advantage of using Ewald summation for dispersion interactions. While the non-PME curve converges to the asymptotic total The summation in Eq. (38b) is handled in the same manner as the reciprocal space sum for electrostatics. Tinker uses the FFTW (Fastest Fourier Transform in the West) package to perform the needed Fourier transforms.60 To speed the calculation and because the dispersion energy decreases quickly with distance, Eq. (38a), the direct space sum, is truncated at a fixed distance. For simple dispersion PME, the choice of direct space cutoff matters very little; one simply chooses a cutoff that balances computational effort between direct space and reciprocal space. For overlap damped dispersion PME, however, some care must be taken with the choice. This is because Eqs. (38a)-(38c) as written do not strictly sum to Eq. (36). This imbalance is caused by the presence of the damping function in the direct space sum, without an equivalent component in the reciprocal space. In practice, however, this is easily overcome with a rational choice of cutoff distance. The function f damp g°es to unity very quickly with distance (much faster than 1/r6 goes to zero), so reasonable cutoff distances are easy to obtain. Figure 12 shows dispersion energy as a function of cutoff distance. £ -8200 P -8400 Overlap Damped Dispersion - standard Overlap Damped Dispersion - PME 4 6 8 10 12 14 16 Cutoff Distance (angstroms) FIG. 12. Cutoff distance convergence of the overlap damped dispersion model. The total dispersion energy of a 36 A water box is shown for the standard and particle mesh Ewald (PME) implementations of the overlap damped dispersion model. The cutoff of the PME implementation refers to the cutoff of the real space summation. Computational details are enumerated in Sec. III. 084115-14 Rackers et al. J. Chem. Phys. 149, 084115 (2018) 1PME Overlap Damped Dispersion ■Standard Overlap Damped Dispersion Standard Buffered 14-7 FIG. 13. Computational effort for overlap damped dispersion. The time to complete 100 iterations of overlap damped dispersion and buffered 14-7 is plotted as a function of cutoff distance. cutoff (angstroms) energy very slowly with cutoff distance, the Ewald sum is converged within the 6 A cutoff distance. The slow convergence of the non-Ewald sum is the reason many molecular mechanics models use 12-16 A cutoffs or van der Waals correction terms for their dispersion interactions. Because our model is not forced to use a longer cutoff distance, it can be faster to compute than standard dispersion models. As a point of reference, in the AMOEBA model, the van der Waals calculations currently comprise 10%-15% of the total calculation time. While this is certainly not the bottleneck for efficiency, it is important to keep this relative cost low. In Fig. 13, we compare energy evaluation times for our PME implementation of the model with the standard implementation for various cutoff distances. Figure 13 shows that for standard implementation cutoff distances of greater than 9 A, the PME implementation of the overlap damped dispersion model provides a performance boost relative to the non-PME implementation. As a point of reference, Fig. 13 also shows the speed of the AMOEBA buffered 14-7 van der Waals functional form with its suggested cutoff distance of 12 A. Even compared to this model, which has no required exponential evaluation, the overlap damped dispersion PME model provides a factor of 2.5 speed increase. The use of particle mesh Ewald summation minimizes the work needed in real space, thus enabling the use of our more complicated and accurate functional form without the loss of computational efficiency. Finally, our model benefits from the utilization of simple combination rules. Using a multiplicative combination rule as indicated in Eq. (33) makes our DPME method exact. It is worth noting that many popular force fields, including all three mentioned in this paper (Amber, CHARMM, and AMOEBA), use additive combining rules for van der Waals parameters. Particle mesh Ewald methods can be used with additive combining rules in an approximate manner first proposed by Erik Lindahl and co-workers.61'62 This method prevents what would otherwise be a discontinuity in the forces at the cutoff distance, but does introduce complexity when switching from direct to reciprocal space combining rules.15 The overlap damped dispersion PME model avoids this by explicitly parameterizing to a multiplicative combining rule. This implementation is certainly not unique or novel, but it is important to the future use of the overlap damped dispersion model in a complete force field because it shows that the benefit of a more physics-based short-range model can be realized at no increase in cost. VI. DISCUSSION AND CONCLUSIONS The universe of possible molecular mechanics models is immense. To discriminate between models and decide which models work best with each other, we must evaluate them based upon the goals they wish to achieve. The goal of our proposed model is use in biomolecular modeling, molecular dynamics simulation, and free energy calculations. To this end, it is important for the model to be accurate, transferable, and interpretable. Accuracy has obvious importance for describing the interactions of biological molecules with fidelity, but transferability and interpretability are no less important. In this last section, we summarize how the overlap damped dispersion model measures on each of these attributes. The data in Table I show that when measured against symmetry adapted perturbation theory, damping is necessary to achieve 1 kcal/mol accuracy for the S 101x7 data set. Interestingly, the overlap damped dispersion model is also shown to be more accurate than the attractive part of the slightly more complex buffered 14-7 potential. The root of this behavior seems to lie with the behavior of the dispersion energy at short range. Figure 6 shows that a damping function is necessary to fit both close-contact and large-separation dimer points. In most force fields, this inaccuracy at short range is handled through the cancellation of error. Relying on such cancellation, however, will not work as force fields become more accurate and, more importantly, is not guaranteed to function favorably across the wide variety of intermolecular interactions that occur in biomolecular applications. The transferability of the overlap damped dispersion model is coupled to this idea of eliminating a reliance on the 084115-15 Rackers et al. J. Chem. Phys. 149, 084115 (2018) cancellation of errors. The most notable feature of the model, aside from its accuracy, is the fact that it has no additional adjustable parameters beyond the simple London dispersion model. The damping function, as presented in Sec. II, is entirely determined by the electrostatic model presented in Ref. 11 with no additional fitting or parameterization. This property of the model suggests two things. First, the overlap damped dispersion model is easily transferable to a range of chemical space because of the limited number of parameters. Evidence of this is shown in Figs. 10 and 11 where the S101-fitted parameters were used to predict the dispersion energy of nucleic acid base stacking interactions. Second, it hints at a physical reality behind the model. The fact that parameters generated through fitting to intermolecular electrostatic interactions, where density overlap is the determining factor in short-range interactions, work well for our dispersion model is a strong indicator that the same phenomenon is driving short-range dispersion. This physical picture of short-range dispersion is what makes the overlap damped dispersion model interpretable. There are many damping functions that can be used to correct for the behavior of the dispersion energy at short range. Several of these damping functions can likely be parameterized to yield results against symmetry adapted perturbation theory that are as accurate as those presented here. What the current model offers over the alternatives is a physical interpretation. In this model, the dispersion interaction is the result of the electrostatic interaction between the instantaneous induced dipoles of two distinct charge distributions. This characteristic of the model is valuable for two reasons. First, it gives us some intuition about the nature of intermolecular interactions. Second, the interpretability of the model makes it easier to extend the model to new areas of chemical space. We make no claim that the 18 atom classes used in this paper will accurately describe all of the variety of chemistries in organic molecules. What is clear, however, is that the interpretation of the damping parameter as a measure of an atom's charge distribution gives a clear path to determining new parameters where necessary. In this way, the overlap damped dispersion model is systematically improvable. As advanced molecular mechanics models evolve, this property will be important to their ongoing development. As models grow to explicitly take into account the short-range interactions between molecules, this dispersion model fits neatly into that framework. Accurately modeling the short-range interactions between molecules is important to making trustworthy predictions on a range of biomolecular problems. Drug binding,63 intrinsically disordered protein behavior,64 and nucleic acid structure65 are all areas where advanced force fields have been shown to be necessary for correct predictions. As models get more complex, the tendency is to accumulate additional parameters and with them empiricism. In the case of dispersion, this paper shows that a simple physical model can be employed that adds no new parameters while reducing the need to rely on the cancellation of errors. Moreover, combined with a dispersion particle mesh Ewald implementation, the evaluation of the necessary equations can be achieved as fast or faster than the standard implementation of simple non-damped models. This yields a simple physically interpretable model ready for the next generation of advanced molecular mechanics models. We are currently working on incorporating this model, along with the previously published charge penetration function, into a complete force field. SUPPLEMENTARY MATERIAL See supplementary material for plots showing the dispersion energy for ten different DNA base steps, including the energy as a function of rise, twist, shift, slide, roll, and tilt. ACKNOWLEDGMENTS J.W.R and RR. gladly acknowledge support from the U.S. National Institutes of Health under Award Nos. NIGMS R01 GM114237 and NIGMS R01 GM106137. J.A.R. wishes to acknowledge a 2018 Phase-I MolSSI Software Fellowship funded by the U.S. National Science Foundation. APPENDIX: HIGHER-ORDER DISPERSION TERMS As indicated in the text, the full dispersion interaction between two real atoms also includes higher-order components that give rise to 1/r8, 1/r10, etc. terms. These terms come from instantaneous higher-order multipole interactions between atoms. Similar to the 1/r6 term, these can be derived from a simple Drude oscillator model of atomic polarizability. As the derivation of the origin of these terms is not readily available in the literature, we present here a derivation that continues the series started in the text. 1. Dipole-quadrupole dispersion The derivation of the dipole-quadrupole dispersion energy starts from Eq. (5) in the text, where, instead of the dipole-dipole energy, the dipole-quadrupole interaction energy now enters into the Schrodinger equation, 1 d2x¥ 1 dZxV 2v M dz2 M dz2 h2 1 1. + — \E--kzf--kzi 2 J '^dipole-quadrupole^^ ~ 0- (Al) For a Drude oscillator dipole, interacting with a linear Drude oscillator quadrupole, the energy of the interaction is given by 3(fn&j + fij&t) Udipole—quadrupole ~ VV^7Uchg—chg ~ 4 ■ (A2) If we assume the magnitude of the dipole and quadrupole moments on i and j to be identical, combining Eqs. (Al) and (A2) yields 1 d2xV 1 dlx¥ 2v M dz2 M dz2 h2 ' j 1 1. 2*? " 2^ ¥ = 0. (A3) We now follow a similar transformation of variables as in the text and define A, = Zi+Zj_ V2 ' Ai = V2 (A4) 084115-16 Rackers et al. J. Chem. Phys. 149, 084115 (2018) and rewrite Eq. (8) as 1_^ + 1^2 / _1 1 M dz M dz h \ 2 2 1 1 where ki - k + ki-k- (A6) Equation (A6) is again a transformed version of the independent harmonic oscillator problem. It can be solved in the same manner giving where E(r) = ^h(ü>i +w2), Cl>2 - (A7) Applying the binomial expansion Vl +x X + ■ the total energy becomes 36ho)0 E(r) - Juüo - „ 9 + ■ (A9) (A10) We then subtract the energy of infinitely separated atoms. This gives the dipole-quadrupole term of the dispersion potential energy = + (All) It should be noted here that the spring constants, k, and frequencies, co, are not the same as those for the dipole-dipole interaction. Therefore, just as with the dipole-dipole term, parameters are introduced to give the dipole-quadrupole dispersion model energy, jjdisper C'CJ (A12) 2. Quadrupole-quadrupole dispersion At the risk of repetition, the quadrupole-quadrupole derivation follows almost exactly the formulation above. The Schrodinger equation now reads 1 d2*¥ 1 d2*¥ 2 I 1,1 + M^M£-2fa<--2fa; .2 M dz2 1 J ~Udipole-quadrupole}^ ~ 0- (A13) In this case, we have two linear Drude oscillator quadrupoles, interacting with each other. The energy of the interaction is given by Uquadrupole-quadrupole ~ ^^^^ Ucng—cng — r r Combining Eqs. (A13) and (A14) yields 1 d2x¥ 1 d2x¥ 2 1 1 2 1 2 6®j®j M~dzJ + M'dz2 + Ü2 \ 2 1 ~ 2 i ~ (A15) The transformation of variables is identical to the dipole-quadrupole case (A14) ¥ = 0. A, = V2 A2 = Zi ~ Zj V2 which gives 1 d2x¥ 1 dlx¥ 2v M dz2 where M dz2 (A16) Y = 0, (A17) k + 6Q2ZiZj k2-k- 6Q2ZiZj (A18) Equation (A 17) again gives us the independent harmonic oscillator problem with the solution E(r) - ^h(cji + cj2), where a>i - a>2 - (A19) (A20) Applying the binomial expansion as before, the total energy becomes E{r) — fujjQ ■ 36Q^z2z2nm 1 J (A21) %rwk2 Subtracting the energy of infinitely separated atoms gives the quadrupole-quadrupole term of the dispersion potential energy E{r) = - 36Q4z2z2hoj0 8k21 10 (A22) Again, the spring constants, k, and frequencies, co, are placeholders specific to this quadrupole-quadrupole interaction. To generalize, parameters are introduced to give u: dispersion (quadrupole-quadrupole) C' c1 „10 ' (A23) the quadrupole-quadrupole dispersion energy. As is apparent from the above sequence of derivations, this pattern of even power dispersion coefficients continues indefinitely for as many higher-order multipole moments as one wishes to include. (We should note that at 1/r10 terms and higher, multiple multipole interactions start to be included in terms. The 1/r10 term, for example, involves a dipole-octopole component as well as quadrupole-quadrupole.) This pattern allows us to extrapolate to the full parameterized dispersion expansion jjdispersion even Qi q] k k z k>6 (A24) 'C. Bergonzo, N. M. Henriksen, D. R. Roe, and T. E. Cheatham, "Highly sampled tetranucleotide and tetraloop motifs enable evaluation of common RNA force fields," RNA 21, 1578-1590 (2015). 2F. Fan, S. Huang, H. Yang, M. Raju, D. Datta, V. B. Shenoy, A. C. T. Van Duin, S. Zhang, and T. Zhu, "Mechanical properties of amorphous LixSi alloys: A reactive force field study," Modell. Simul. Mater. Sci. Eng. 21, 074002 (2013). 3J. G. McDaniel and J. Schmidt, "Next-generation force fields from symmetry-adapted perturbation theory," Annu. Rev. Phys. Chem. 67, 467-488 (2016). 084115-17 Rackers eř a/. J. Chem. Phys. 149, 084115 (2018) 4M. Tafipolsky and K. Ansorg, "Toward a physically motivated force field: Hydrogen bond directionality from a symmetry-adapted perturbation theory perspective," J. Chem. Theory Comput. 12, 1267-1279 (2016). 5J. Schmidt, K. Yu, and J. G. McDaniel, "Transferable next-generation force fields from simple liquids to complex materials," Acc. Chem. Res. 48, 548-556 (2015). 6 A. J. Misquitta, R. Podeszwa, B. Jeziorski, and K. Szalewicz, "Inter-molecular potentials based on symmetry-adapted perturbation theory with dispersion energies from time-dependent density-functional calculations," J. Chem. Phys. 123, 214103 (2005). 7D. E. Taylor, F. Rob, B. M. Rice, R. Podeszwa, and K. Szalewicz, "A molecular dynamics study of l,l-diamino-2,2-dinitroethylene (FOX-7) crystal using a symmetry adapted perturbation theory-based intermolecular force field," Phys. Chem. Chem. Phys. 13, 16629-16636 (2011). SM. S. Gordon, L. V. Slipchenko, H. Li, and J. H. Jensen, "The effective fragment potential: A general method for predicting intermolecular interactions," Annu. Rep. Comput. Chem. 3, 177-193 (2007). 9G. A. Cisneros, "Application of Gaussian electrostatic model (GEM) distributed multipoles in the AMOEBA force field," J. Chem. Theory Comput. 8, 5072-5080 (2012). 10L. El Khoury, S. Naseem-Khan, K. Kwapien, Z. Hobaika, R. G. Maroun, J.-P. Piquemal, and N. Gresh, "Importance of explicit smeared lone-pairs in ansiotropic polarizable molecular mechanics. Torture track angular tests for exchange-repulsion and charge transfer contributions," J. Comput. Chem. 38, 1897-1920 (2017). "j. A. Rackers, Q. Wang, C. Liu, J.-P. Piquemal, P. Ren, and J. W. Ponder, "An optimized charge penetration model for use with the AMOEBA force field," Phys. Chem. Chem. Phys. 19, 276-291 (2017). 12M. Kolář, T. Kubař, and P. Hobza, "On the role of London dispersion forces in biomolecular structure determination," J. Phys. Chem. B 115, 8038-8046 (2011). 13K. E. Riley and P. Hobza, "Investigations into the nature of halogen bonding including symmetry adapted perturbation theory analyses," J. Chem. Theory Comput. 4, 232-242 (2008). 14K. E. Riley and P. Hobza, "The relative roles of electrostatics and dispersion in the stabilization of halogen bonds," Phys. Chem. Chem. Phys. 15, 17742-17751 (2013). 15A. N. Leonard, A. C. Simmonett, F. C. Pickard, J. Huang, R. M. Venable, J. B. Klauda, B. R. Brooks, and R. W. Pastor, "Comparison of additive and polarizable models with explicit treatment of long-range Lennard-Jones interactions using alkane simulations," J. Chem. Theory Comput. 14, 948-958 (2017). 16L.-P. Wang, K. A. McKiernan, J. Gomes, K. A. Beauchamp, T. Head-Gordon, J. E. Rice, W. C. Swope, T. J. Martinez, and V. S. Pande, "Building a more predictive protein force field: A systematic and reproducible route to AMBER-FB15," J. Phys. Chem. B 121, 4023-4039 (2017). 17J. Huang, S. Rauscher, G. Nawrocki, T. Ran, M. Feig, B. L. de Groot, H. Grabmiiller, and A. D. MacKerell, "CHARMM36m: An improved force field for folded and intrinsically disordered proteins," Nat. Methods 14, 71-73 (2017). 1ST A. Halgren, "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 (1992). "P. Ren and J. W. Ponder, "Polarizable atomic multipole water model for molecular mechanics simulation," J. Phys. Chem. B 107, 5933-5947 (2003). 20F. C. Brooks, "Convergence of intermolecular force series," Phys. Rev. 86, 92-97 (1952). 21R. AhLrichs, R. Penco, and G. Scoles, "Intermolecular forces in simple systems," Chem. Phys. 19, 119-130 (1977). 22K. T. Tang and J. P. Toennies, "An improved simple model for the van der Waals potential based on universal damping functions for the dispersion coefficients," J. Chem. Phys. 80, 3726-3741 (1984). 23S. Grimme, "Density functional theory with London dispersion corrections," Wiley Interdiscip. Rev.: Comput. Mol. Sci. 1, 211-228 (2011). 24Q. Wu and W. Yang, "Empirical correction to density functional theory for van der Waals interactions," J. Chem. Phys. 116, 515-524 (2002). 25J.-D. Chai and M. Head-Gordon, "Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections," Phys. Chem. Chem. Phys. 10, 6615-6620 (2008). 26E. R. Johnson and A. D. Becke, "A post-Hartree-Fock model of intermolecular interactions," J. Chem. Phys. 123, 024101 (2005). 27L. V. Slipchenko and M. S. Gordon, "Damping functions in the effective fragment potential method," Mol. Phys. 107, 999-1016 (2009). 2SI. Adamovic and M. S. Gordon, "Dynamic polarizability, dispersion coefficient Cg and dispersion energy in the effective fragment potential method," Mol. Phys. 103, 379-387 (2005). 29P. Verma, B. Wang, L. E. Fernandez, and D. G. Trahlar, "Physical molecular mechanics method for damped dispersion," J. Phys. Chem. A 121, 2855-2862 (2017). 30G. C. Maitland, M. Rigby, E. B. Smith, and W. A. Wakeham, Intermolecular Forces (Oxford University Press, 1981). 31 A. J. Misquitta and A. J. Stone, "Dispersion energies for small organic molecules: First row atoms," Mol. Phys. 106, 1631-1643 (2008). 32N. Gresh, D. Perahia, B. de Courcy, J. Foret, C. Roux, L. El Khoury, J.-P. Piquemal, and L. Salmon, "Complexes of a Zn-metalloenzyme binding site with hydroxamate-containing ligands. A case for detailed benchmarkings of polarizable molecular mechanics/dynamics potentials when the experimental binding structure is unknown," J. Comput. Chem. 37, 2770-2782 (2016). 33L. Bytautas and K. Ruedenberg, "Correlation energy and dispersion interaction in the ab initio potential energy curve of the neon dimer," J. Chem. Phys. 128, 214308 (2008). 34A. D. Buckingham, "Theory of long-range dispersion forces," Discuss. Faraday Soc. 40, 232-238 (1965). 35P. Xu, F. Zahariev, and M. S. Gordon, "The R-7 dispersion interaction in the general effective fragment potential method," J. Chem. Theory Comput. 10, 1576-1587 (2014). 36E. B. Guidez, P. Xu, and M. S. Gordon, "Derivation and implementation of the gradient of the R-7 dispersion interaction in the effective fragment potential method," J. Phys. Chem. A 120, 639-647 (2016). 37 Q. Wang, J. A. Rackers, C. He, R. Qi, C. Narth, L. Lagardere, N. Gresh, J. W. Ponder, J.-P. Piquemal, and P. Ren, "General model for treating short-range electrostatic penetration in a molecular mechanics force field," J. Chem. Theory Comput. 11, 2609-2618 (2015). 38 C. Narth, L. Lagardere, E. Polack, N. Gresh, Q. Wang, D. R. Bell, J. A. Rackers, J. W. Ponder, P. Y. Ren, and J. P. Piquemal, "Scalable improvement of SPME multipolar electrostatics in anisotropic polarizable molecular mechanics using a general short-range penetration correction up to quadrapoles," J. Comput. Chem. 37, 494-506 (2016). 39C. Coulson, "Two-centre integrals occurring in the theory of molecular structure," Math. Proc. Cambridge 38, 210-223 (1942). 40P. Ren, SAPT Database Between Organic Molecules, Protein Side Chain Analogs, available from: http://biomol.bme.utexas.edu/~ch38988/sl01x7, 2015. 41B. Jeziorski, R. Moszynski, and K. Szalewicz, "Perturbation theory approach to intermolecular potential energy surfaces of van der Waals complexes," Chem. Rev. 94, 1887-1930 (1994). T. M. Parker, L. A. Burns, R. M. Parrish, A. G. Ryno, and C. D. Sherrill, "Levels of symmetry adapted perturbation theory (SAPT). I. Efficiency and performance for interaction energies," J. Chem. Phys. 140, 094106 (2014). 43T. H. Dunning, Jr., "Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen," J. Chem. Phys. 90, 1007-1023 (1989). 44R. A. Kendall, T. H. Dunning, Jr., and R. J. Harrison, "Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions," J. Chem. Phys. 96, 6796-6806 (1992). 45 A. Halkier, T. Helgaker, P. J0rgensen, W. Klopper, H. Koch, J. Olsen, and A. K. Wilson, "Basis-set convergence in correlated calculations on Ne, N2, and H20," Chem. Phys. Lett. 286, 243-252 (1998). 46E. G. Hohenstein and C. D. Sherrill, "Density fitting of intramonomer correlation effects in symmetry-adapted perturbation theory," J. Chem. Phys. 133, 014101 (2010). J. M. Turney, A. C. Simmonett, R. M. Parrish, E. G. Hohenstein, F. A. Evan-gelista, J. T. Fermann, B. J. Mintz, L. A. Burns, J. J. Wilke, and M. L. Abrams, "Psi4: An open-source ab initio electronic structure program," Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2, 556-565 (2012). R. M. Parrish, L. A. Burns, D. G. A. Smith, A. C. Simmonett, A. E. DePrince III, E. G. Hohenstein, U. Bozkaya, A. Y. Sokolov, R. Di Remi-gio, and R. M. Richard, "Psi4 1.1: An open-source electronic structure program emphasizing automation, advanced libraries, and interoperability," J. Chem. Theory Comput. 13, 3185-3197 (2017). T. M. Parker and C. D. Sherrill, "Assessment of empirical models versus high-accuracy ab initio methods for nucleobase stacking: Evaluating the importance of charge penetration," J. Chem. Theory Comput. 11, 4197-4204 (2015). 50J. A. Rackers, Tinker Github Repository, DPME Branch, available from: https://github.com/JoshRackers/tinker/tree/dpme, 2017. 42' 47 4s 49, 084115-18 Rackers er al. J. Chem. Phys. 149, 084115 (2018) R. Qi, Q. Wang, and P. Ren, "General van der Waals potential for common organic molecules," Bioorg. Med. Chem. 24, 4911-4919 (2016). S. Piana, A. G. Donchev, P. Robustelli, and D. E. Shaw, "Water dispersion interactions strongly influence simulated structural properties of disordered protein states," J. Phys. Chem. B 119, 5113-5123 (2015). 'W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, "Comparison of simple potential functions for simulating liquid water," J. Chem. Phys. 79, 926-935 (1983). D. H. Brookes and T. Head-Gordon, "Family of oxygen-oxygen radial dirsribution functions for water," J. Phys. Chem. Lett. 6,2938-2943 (2015). 1Y. Mao, O. Demerdash, M. Head-Gordon, and T. Head-Gordon, "Assessing ion-water interactions in the AMOEBA force field using energy decomposition analysis of electronic structure calculations," J. Chem. Theory Comput. 12, 5422-5437 (2016). 'O. Demerdash, Y. Mao, T. Liu, M. Head-Gordon, and T. Head-Gordon, "Assessing many-body contributions to intermolecular interactions of the AMOEBA force field using energy decomposition analysis of electronic structure calculations," J. Chem. Phys. 147, 161721 (2017). J. Rezac, K. E. Riley, and P. Hobza, "S66: A well-balanced database of benchmark interaction energies relevant to biomolecular structures," J. Chem. Theory Comput. 7, 2427-2438 (2011). 'W. L. Jorgensen and D. L. Severance, "Aromatic-aromatic interactions: Free energy profiles for the banzene dimer in water, chloroform, and liquid benzene," J. Am. Chem. Soc. 112, 4768-4774 (1990). U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Ped-ersen, "A smooth particle mesh Ewald method," J. Chem. Phys. 103, 8577-8593 (1995). 60M. Frigo and S. G. Johnson, "The design and implementation of FFTW3," Proc. IEEE 93, 216-231 (2005). 61C. L. Wennberg, T. Murtola, B. Hess, and E. Lindahl, "Lennard-Jones lattice summation in bilayer simulations has critical effects on surface tension and lipid properties," J. Chem. Theory Comput. 9, 3527-3537 (2013). 62C. L. Wennberg, T. Murtola, S. Pall, M. J. Abraham, B. Hess, andE. Lindahl, "Direct-space corrections enable fast and accurate Lorentz-Berthelot combination rule Lennard-Jones lattice summation," J. Chem. Theory Comput. 11, 5737-5746 (2015). 63D. R. Bell, R. Qi, Z. Jing, J. Y. Xiang, C. Mejias, M. J. Schnieders, J. W. Ponder, and P. Ren, "Calculating binding free energies of host-guest systems using the AMOEBA polarizable force field," Phys. Chem. Chem. Phys. 18, 30261-30269 (2016). 64J. Huang and A. D. MacKerell, "Force field development and simulations of intrinsically disordered proteins," Curr. Opin. Struct. Biol. 48, 40-48 (2018). 65T. M. Parker, E. G. Hohenstein, R. M. Parrish, N. V. Hud, and C. D. Sherrill, "Quantum-mechanical analysis of the energetic contributions to jt stacking in nucleic acids versus rise, twist, and slide," J. Am. Chem. Soc. 135, 1306-1316(2013).