Inclusion of charge and polarizability fluxes provides needed physical accuracy in molecular mechanics force fields

Chemical Physics Letters 429 (2006) 628–632 www.elsevier.com/locate/cplett

Inclusion of charge and polarizability fluxes provides needed physical accuracy in molecular mechanics force fields Kim Palmo, Berit Mannfors, Noemi G. Mirkin, Samuel Krimm

*

Biophysics Research Division and Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA Received 20 June 2006; in final form 17 August 2006 Available online 25 August 2006

Abstract Physical accuracy in MM calculations requires maximally correct forces in addition to structures and energies. This requires inclusion of charge and polarizability fluxes in the energy function. Using our spectroscopically determined force field, which is designed to include these, we present three examples of the improved physical predictions that result from the incorporation of charge fluxes: the reproduction of the water angle opening on going from the isolated molecule to the liquid, the reproduction of the w peptide torsion potential with only a single threefold (almost zero barrier) Fourier term, and reproduction of the quantum-mechanical MD u, w map of a dipeptide analog.  2006 Elsevier B.V. All rights reserved.

1. Introduction Since the earliest efforts at producing molecular mechanics (MM) energy functions for computational simulation, the goal has been to include in these classical descriptions enough physics content to reproduce experimental observations such as molecular structures and energies. To reach this goal it has, until recently, been deemed sufficient to base these MM force fields on simple models, using diagonal valence terms, Coulomb fixed charge electrostatics, Lennard-Jones or Buckingham van der Waals (vdW) potentials, and Fourier series descriptions of torsion potentials [1]. The energy parameters were optimized primarily to observed liquid and crystal properties, thus giving rise to the concept of ‘empirical force field’. As the necessity of improved predictive accuracy became apparent, and optimization to higher-level quantum-mechanical (QM) data became more feasible, the force fields have been frequently reoptimized. Although some new features have been introduced, such as electrical atomic multipoles (or off-atom charges) and three-body electrostatics through *

Corresponding author. Fax: +1 734 764 3323. E-mail address: [email protected] (S. Krimm).

0009-2614/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2006.08.087

polarization, these efforts have mostly been made within the empirical force field strategy, and have not particularly emphasized consistency, transferability, and compatibility of the terms representing the different types of interactions. This is closely related to the significant problem that another basic physical component that must be an integral part of a successful MM function has been neglected [2]. What has been missing was pointed out almost 70 years ago by Feynman [3]: ‘Many of the problems of molecular structure are concerned essentially with forces.’, which are treated only ‘indirectly through the agency of energy and its changes with changing configuration of the molecule.’ In fact, it is possible for an MM function to give correct relative energies without obtaining correct forces, the latter depending critically on the shape of the potential surface. Since the forces are crucial for producing realistic molecular dynamics (MD) simulations, it is important that the incorporation of accurate force description be a basic pillar in developing physically accurate MM force fields. The only current methodology to adopt this philosophy, and to stress the importance of self-consistency in the energy function, is our spectroscopically determined force field (SDFF) [4]. It does so by incorporating the physical effects needed to reproduce normal mode frequencies,

K. Palmo et al. / Chemical Physics Letters 429 (2006) 628–632

eigenvectors, and infrared and Raman intensities. Vibrational information has been underused in previous MM force field development, although spectra can be accurately measured and the modes reflect the preferred motional patterns (dynamics) of molecular equilibrium structures. Spectroscopic accuracy is guaranteed by the SDFF analytical transformation from QM structures and Hessians, which preserves the QM geometry and force constants, including all of the significant valence cross-terms [5]. In this Letter, we demonstrate an important source of the gap in reproduction of complete forces in current MM functions, outline how the SDFF, a ‘physics-driven’ MM energy function, closes this gap, and present three examples of the improved results of including conformation-dependent atomic charges (a previously ignored effect) in the energy function: the reproduction of the opening of the water angle on going from the gas to the liquid phase (current MM functions predicting the opposite); the reproduction of the w peptide torsion potential with only a single threefold (almost zero barrier) Fourier term; and the reproduction of the QM MD u, w map of a dipeptide analog (in contrast to the predictions of current standard force fields). 2. Complete electrical forces require charge and polarizability fluxes It is instructive to recall the origin of electrical forces experienced by a molecule. The potential energy of a neutral molecule in a general electric field, E (due to an internal charge distribution plus an external field, if present), can be written as 1 V ðEÞ ¼ lstat  E  lind  E ð1Þ 2 P where lstat ¼ i qi ri is the static dipole moment resulting from the atomic partial charges qi at locations ri and lind = aE is the induced dipole moment resulting from the molecule’s polarizability a. The forces acting on the molecule due to the field are then (in Cartesian coordinates xk) 

oV ðEÞ olstat 1 olind 1 oE ¼ Eþ  E þ ðlstat þ lind Þ  oxk 2 oxk 2 oxk oxk

ð2Þ

From the above definitions, the derivatives of the dipoles lstat and lind are, respectively, olstat X ori X oqi ¼ qi þ ri ð3Þ oxk oxk oxk i i and olind oE oa ¼a þ E oxk oxk oxk

ð4Þ

In a vibrating molecule, the intensity of an IR absorption band is proportional to the squares of the o(lstat + lind)/ oxk that contribute to the normal mode, and the intensity of Raman scattering is proportional to the squares of the oa/oxk of the mode. In standard MM functions, the first term on the righthand side of (3) is calculated but the second term would

629

not be, because no oqi/oxk are included in the force field. In non-polarizable force fields there is, of course, no olind/oxk, but even in polarizable force fields only the first term in (4) would be calculated. The SDFF incorporates oqi/oxk terms, i.e., charge fluxes [4], and makes provision (presently being implemented) for the inclusion of polarizability fluxes, oa/oxk. Thus, the SDFF starts out with the goal of providing correct forces (of course as well as correct structures and energies). It achieves this goal by its attention to reproducing vibrational frequencies and intensities to spectroscopic standards [4]. In other words, the accurate prediction of vibrational frequencies and intensities is a defining test of whether a force field has achieved substantively correct force reproduction. It is important to note that the charge fluxes represent the existence of conformation-dependent charges, i.e., the well-known intrinsic electrical anharmonicity associated with deformation of internal coordinates [6–9], and should not be confused with the ‘fluctuating charge’ polarization model [10]. Polarization may indeed be important in the description of electrostatic interactions, and many models have been proposed for introducing this effect [11]. We have formulated a non-iterative group polarizability model for this purpose [12], which is accurate, stable, and easy to parallelize. However, there is evidence that the effect of charge fluxes can be comparable to or even significantly larger than that of polarization [8,13]. We have found that the best way to model charge fluxes from ab initio dipole derivatives is to require a balanced charge treatment of intramolecular Coulomb interactions [14] and to possibly incorporate dipole flux representations [15]. With regard to polarizability flux, the parameters can be determined as follows. For each model molecule in a set, we start by calculating the ab initio molecular polarizability derivative tensor oa/oxk in Cartesian coordinates. We then use our MM polarization model (based on the same ab initio level) to compute the contribution to the Cartesian polarizability derivatives resulting from the static polarizability ao only, i.e., oao/oxk (note that this is not zero if the molecule contains anisotropic groups). By subtracting that contribution from the ab initio (total) polarizability derivative tensor, we obtain the matrix (oa 0 /oxk = oa/ oxk  oao/oxk) that needs to be reproduced by the MM polarizability flux parameters. Transforming to internal coordinates, si, we get a set of 3 · 3 matrices Ai, defined by oa0 X oa0 oxk ¼ : ð5Þ Ai ¼ osi oxk osi k Each Ai matrix is still associated with the same (arbitrary) Cartesian orientation as the original a, and its elements are, therefore, not suitable to be stored as transferable MM parameters. It needs to be further transformed into a local orthonormal coordinate system with base vectors e1, e2, e3 that rotate with si. (For example, an angle coordinate is chosen to have e1 along the bisector, e2 perpendicular to the bisector but in the angle plane, and e3 = e1 · e2.) The transformation that takes Ai into its local basis is given by

630

K. Palmo et al. / Chemical Physics Letters 429 (2006) 628–632

A0i ¼ Pti Ai Pi

ð6Þ

where Pi is a matrix whose columns are e1 e2 e3, and the superscript t denotes transpose. The elements of A0i can be directly used as MM parameters (a maximum of six for each internal coordinate). These parameters describe how the polarizability changes locally when the internal coordinate si is deformed. A0i contains off-diagonal elements unless e1, e2, e3 happen to be its eigenvectors. To use the parameters and calculate the MM oa/oxk, the transformations described above are carried out in reverse order. 3. Opening of the water angle The HOH water angle is experimentally known to be slightly larger (1–2 on the average) in the liquid phase than in the gas phase [16]. This behavior is in agreement with QM based calculations [17,18] but in disagreement with previous force field based MD simulations. Such force fields instead predict a closing of the HOH angle by 4 even with the inclusion of polarization [19]. However, our SDFF water model [20] correctly reproduces the opening of the HOH angle when charge fluxes are included. The SDFF water model also correctly predicts the small closing of the HOH angle with increasing temperature. In fact, we have found that the erroneous gas-to-liquid closing trend of the HOH angle can be reversed (and corrected) for virtually any water model by using appropriate charge fluxes. Three different charge flux parameters are involved (all of them influencing the flow of charge along the OH bonds as a function of geometry): (1) the flux along an OH bond when that bond is deformed, (2) the flux along an OH bond when the other OH bond is deformed, and (3) the flux along the OH bonds when the HOH angle is deformed. Table 1 Gas-to-liquid deformation of the HOH water angle at 300 K

SDFF TIP3P-F QMa a

w/charge flux

w/o charge flux

1.93 1.95 1.93

3.43 4.74

Ref. [17].

Not surprisingly, the last mentioned type (angle charge flux) has by far the greatest impact on the gas-to-liquid HOH angle change. We have investigated these properties by doing liquid water MD simulations with our (non-polarizable) SDFF water model and comparing these results with those of a QM calculation on liquid water [17] and of a flexible version of TIP3P water (named TIP3P-F). The TIP3P [21] water molecule was made flexible by applying the same stretching, bending, and cross-term force constants as in the SDFF model, while preserving the original TIP3P structure as the intrinsic equilibrium geometry. The TIP3P electrostatic and van der Waals potentials were also left unchanged. In 100 ps MD simulations (0.5 fs time step) on a droplet of 120 equilibrated water molecules (subject to soft spherical boundary conditions [22]) we first determined charge fluxes, for both the SDFF and TIP3P-F models, that approximately yield the QM derived opening of the HOH angle at 300 K. (The required optimization from starting QM dipole derivative values probably reflects the present lack of inclusion at this stage of polarizability fluxes, as well as the assumed interaction centers on the atoms for the charge fluxes.) This is shown in Table 1. We then determined the temperature dependence of the angle deformation in additional 100 ps simulations at 275–375 K. The results are shown in Fig. 1. As can be seen, the small closing trend of the HOH angle with increasing temperature, as given by QM, is well captured by the SDFF and TIP3P-F models when charge fluxes are used. Without charge fluxes, the temperature behavior of the MM models is incompatible with the QM results. Test calculations show that flexible versions of TIP4P and TIP5P (and probably any water model) behave in the same way in this regard. These simulations confirm that the charge flux is a nonnegligible physical property that must be included in MM force fields if a realistic representation of molecular forces is desired. We expect that these parameters will need to be re-optimized when polarizability fluxes are also included. The importance of structure dependent charges has also been pointed out from related modeling of the water angle change [23].

-3.0

2.2

HOH deform. (deg)

HOH deform. (deg)

2.4

2.0 1.8 1.6 1.4 1.2 260

280

300

320 340 Temp (K)

360

380

-3.5 -4.0 -4.5 -5.0 260

280

300

320

340

360

380

Temp (K)

Fig. 1. Temperature dependence of the HOH angle deformation in liquid water: (a) w/charge flux and (b) w/o charge flux. Notation: () QM (Ref. [17]); (h) SDFF; (d) TIP3P-F.

K. Palmo et al. / Chemical Physics Letters 429 (2006) 628–632

4. w Torsion potential Proper parameterization of the u (CNCaC) and w (NCaCN) torsion potentials is critical to representing the correct conformational flexibility in the protein chain, yet thus far only ad hoc modifications of standard Fourier representations have been implemented [11,24–26]. Such ‘corrections’, which usually are accompanied by non-transferability of parameters, may actually hide inaccuracies elsewhere in the potential. We also encountered this problem in trying to obtain expected transferable intrinsic u and w torsion parameters. This led us to a more detailed ab initio study of Nmethylacetamide (NMA) and alkyl chain analogs, such as CH3CH2CONHCH3, CH3CONHCH2CH3, and CH3CH(CH3)CONHCH3. [27]. The initial indications were twofold: (1) the alkyl chain atomic charges depend significantly on conformation, i.e., charge fluxes cannot be neglected, and (2) using our charge flux model, u, w torsion potential energy barriers consistent with ab initio can be obtained with the same intrinsic single 3-fold terms as in the Cand N-methyl torsions of NMA (and only of a few tenths of a kcal/mol). An example for the w potential is shown in Fig. 2. As can be seen, the torsion energy variations generally follow the charge variations (Figs. 2a, b), thus indicating that the latter contribute significantly to the torsion energy. Our SDFF calculations confirm this (Fig. 2c). Similar results apply to the u potential. As supported by the results presented in the next section, there is reason to believe that the charge fluxes also account, in

631

a natural and physical way, for the so-called u/w crossterm effects in peptides [28]. 5. Peptide conformational dynamics In the absence of experimental data, the ultimate validation of an MM force field is close agreement with QM results (although possible only for small systems), and the most demanding application is molecular dynamics. Encouraged by the experiences of the alkyl substituted NMA chains, we have optimized charge flux parameters to the entire u, w energy surface of an alanine dipeptide analog (end methyl groups replaced by hydrogen atoms). We chose this molecule because QM molecular dynamics (at the density functional BP/DZVP level of theory, 85 ps, 1.21 fs time step) was compared with classical MD as given by a few standard force fields [29]. The resulting u, w population patterns are shown in Fig. 3 for QM (left) and CHARMM (middle) as well as our corresponding results for the SDFF (right). As can be seen, the QM and SDFF trajectories cover almost the same regions of space and contain barrier crossings between the two minima (C7eq and C5), whereas the CHARMM trajectory remains close to the C7eq minimum at all times (reportedly for nanoseconds, as did the trajectories of the other standard force fields). This example shows that the conformational dynamics produced by the SDFF is superior to that of current force fields, and is far closer to the physical reality as a result of better force prediction and a more accurate description of the motional modes.

Fig. 2. Ab initio (MP2/6-31++G(d,p)) and SDFF results for CH3CaH(CH3)CONHCH2CH3 [27]: (a) Variation of the CHELPG charge on the Ca atom as a function of the w (NCCaC) torsion angle. (b) Ab initio w torsion energy. (c) SDFF w torsion energy.

Fig. 3. Population maps of QM (left), CHARMM (middle), and SDFF (right) from MD trajectories of an alanine dipeptide analog (QM and CHARMM maps reprinted with permission from Ref. [29]).

632

K. Palmo et al. / Chemical Physics Letters 429 (2006) 628–632

6. Conclusions As described here, allowing the molecular charge distribution (and ultimately the polarizability) to vary with conformation in a way consistent with QM is a necessary step toward the goal of obtaining more realistic forces from MM energy functions. This is true whether the charge variations are small, as in the case of water [20,23], or relatively large, as in the case of pyramidalization of the peptide group nitrogen atom [30]. The improved predictions of physical properties resulting from this heretofore neglected inclusion can no longer be considered ‘small effects,’ and, with the demonstrated successful balanced charge treatment of charge fluxes in the SDFF [4,14], MM force fields must now incorporate such physical effects in the electrostatic interactions. Acknowledgements This research was supported by NSF Grant Nos. MCB0212232 and DMR0239417. Additional funding was provided by the Michigan Universities Commercialization Initiative. References [1] J.W. Ponder, D.A. Case, in: F.M. Richards, D.S. Eisenberg, J. Kuriyan (Eds.), Advances in Protein Chemistry, Elsevier, Amsterdam, 2003. [2] M.E. Colvin et al., J. Mol. Struct. (THEOCHEM) 764 (2006) 1. [3] R.P. Feynman, Phys. Rev. 56 (1939) 340. [4] K. Palmo, B. Mannfors, N.G. Mirkin, S. Krimm, Biopolymers 68 (2003) 383.

[5] K. Palmo, L.-O. Pietila, S. Krimm, J. Comput. Chem. 12 (1991) 385. [6] G. Herzberg, Infrared and Raman Spectra, D. Van Nostrand, New York, 1945, p. 241. [7] D.E. Williams, Biopolymers 29 (1990) 1367. [8] U. Dinur, J. Phys. Chem. 94 (1990) 5669. [9] C.A. Reynolds, J.W. Essex, W.G. Richards, J. Am. Chem. Soc. 114 (1992) 9075. [10] S.W. Rick, S.J. Stuart, B.J. Berne, J. Chem. Phys. 101 (1994) 6141. [11] A.D. MacKerell Jr., J. Comput. Chem. 25 (2004) 1584. [12] K. Palmo, S. Krimm, Chem. Phys. Lett. 395 (2004) 133. [13] U. Dinur, J. Phys. Chem. 95 (1991) 6201. [14] K. Palmo, B. Mannfors, S. Krimm, Chem. Phys. Lett. 369 (2003) 367. [15] K. Palmo, S. Krimm, J. Comput. Chem. 19 (1998) 754. [16] K. Ichikawa, Y. Kameda, T. Yamaguchi, H. Wakita, M. Misawa, Mol. Phys. 73 (1991) 79. ˚ strand, J. Phys. Chem. A 101 (1997) 10039. [17] T.M. Nymand, P.-O. A [18] P.L. Silvestrelli, M. Parrinello, J. Chem. Phys. 111 (1999) 3572. [19] P. Ren, J.W. Ponder, J. Phys. Chem. 107 (2003) 5933. [20] K. Palmo, B. Mannfors, S. Krimm, in press. [21] M.W. Mahoney, W.L. Jorgensen, J. Chem. Phys. 112 (2000) 8910. [22] R. Sankararamakrishnan, K. Konvicka, E.L. Mehler, H. Weinstein, Int. J. Quant. Chem. 77 (2000) 174. [23] G.S. Fanourgakis, S.S. Xantheas, J. Chem. Phys. 124 (2006) 174504. [24] S. Dasgupta, W.B. Hammond, W.A. Goddard III, J. Am. Chem. Soc. 118 (1996) 12291. [25] Z.-X. Wang, W. Zhang, C. Wu, H. Lei, P. Cieplak, Y. Duan, J. Comput. Chem. 27 (2006) 781. [26] A.E. Garcia, K.Y. Sanbonmatsu, Proc. Nat. Acad. Sci. USA 99 (2002) 2782. [27] K. Palmo, N.G. Mirkin, S. Krimm, in press. [28] A.D. MacKerell Jr., M. Feig, C.L. Brooks III, J. Am. Chem. Soc. 126 (2004) 698. [29] D. Wei, H. Guo, D.R. Salahub, Phys. Rev. E 64 (2001) 011907. [30] B.E. Mannfors, N.G. Mirkin, K. Palmo, S. Krimm, J. Phys. Chem. A 107 (2003) 1825.