QM/MM-Based Fitting of Atomic Polarizabilities for Use in Condensed-Phase Biomolecular Simulation C. Ruben Vosmeer, † Arië n S. Rustenburg, † Julia E. Rice, ‡ Hans W. Horn, ‡ William C. Swope, ‡ and Daan P. Geerke* ,† † Leiden/Amsterdam Center for Drug Research, Division of Molecular and Computational Toxicology, Department of Chemistry and Pharmaceutical Sciences, Faculty of Sciences, VU University Amsterdam, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands ‡ IBM Almaden Research Center, 650 Harry Road, San Jose, California 95120, United States ABSTRACT: Accounting for electronic polarization effects in biomolecular simulation (by using a polarizable force field) can increase the accuracy of simulation results. However, the use of gas-phase estimates of atomic polarizabilities α i usually leads to overpolarization in condensed-phase systems. In the current work, a combined QM/MM approach has been employed to obtain condensed-phase estimates of atomic polarizabilities for water and methanol (QM) solutes in the presence of (MM) solvents of different polarity. In a next step, the validity of the linear response and isotropy assumptions were evaluated based on the observed condensed-phase distributions of α i values. The observed anisotropy and low average values for the polarizability of methanol’s carbon atom in polar solvents was explained in terms of strong solute-solvent interactions involving its adjacent hydroxyl group. Our QM/MM estimates for atomic polarizabilities were found to be close to values used in previously reported polarizable water and methanol models. Using our estimate for α O of methanol, a single set of polarizable force field parameters was obtained that is directly transferable between environments of different polarity. 1. INTRODUCTION Since their emergence in the 1950s and 1960s, classical molecular simulation methods have been successfully applied to study physical, chemical, and biological systems, ranging from pure liquids to large complexes such as proteins and cell membranes. 1,2 The increase of computational power has made possible longer simulation times (on the order of milli- seconds) 3,4 for systems and processes of high complexity (involving, e.g., protein activation or membrane transport). 5,6 However, the reliability of molecular simulation outcomes still heavily depends on the accuracy of the description of nonbonded interatomic interactions by the applied potential energy model or force field. Most force fields used in biomolecular simulation employ a pairwise additive Coulomb function to describe the electrostatic interactions, 7,8 leaving out electrostatic effects such as electronic polarization. 9 In past decades, several nonpolarizable biomolecular force fields (such as Amber, 10-12 CHARMM, 13-15 OPLS, 16-18 and GRO- MOS 19-21 ) have been under constant development and have reached high accuracy in describing thermodynamic and other relevant properties of biomolecular systems. However, the omission of explicitly accounting for electronic polarization may lead to limitations in correctly describing the balance of molecular interactions in environments of different polarity, 22 which plays a key role in determining thermodynamic equilibria typically studied in biomolecular simulation. 2 For example, Oostenbrink et al. showed the difficulty of generating a single set of force-field parameters for condensed-phase protein simulation, resulting in two distinct (nonpolarizable) versions of the GROMOS force field, GROMOS 53A5 and 53A6. 20 As another example, polarizable force fields typically bring clear improvement in correctly describing kinetic and dielectric properties of polar solvents, such as water 23 or N-methyl acetamide. 24 The inducible dipole (ID) or point polarizable dipole (PPD) model 25-27 and the Drude oscillator (DO) or charge-on-spring (COS) approach 28,29 are two of the most widely applied methods to include explicit electronic polarization in polarizable force fields. A third one, designated as the fluctuating charge (FQ) model, 30,31 includes electronic polarization effects by allowing atomic point charge distributions of molecules (or molecular building blocks) to adapt to the external electric field. In contrast, the PPD/ID and DO/COS methods assign an inducible dipole moment μ⃗ i to (heavy) atoms i, which adapts its magnitude and direction in response to the external electric field E ⃗ i , depending on its polarizability α i . Linear response is usually assumed, and the μ⃗ i ′s are determined using eq 1. 9,32 α πε μ⃗ = ⃗ E (4 ) i i i 0 (1) In the PPD method, μ⃗ i is represented by a point dipole, whereas the COS method introduces two additional point charges attached by a spring with a force constant directly depending on α i . In both methods, additional interactions are to be evaluated during simulation (i.e., those involving the introduced dipole moments or point charges, respectively), which raises computational costs when compared to the use of nonpolarizable force fields. In both the COS/DO and PPD/ID models, as well as in the FQ method, an iterative scheme 33,34 (or an extended Lagrangian formalism) 31,35-37 is usually employed to energy minimize (or follow in time) the molecular Special Issue: Wilfred F. van Gunsteren Festschrift Received: February 1, 2012 Published: April 27, 2012 Article pubs.acs.org/JCTC © 2012 American Chemical Society 3839 dx.doi.org/10.1021/ct300085z | J. Chem. Theory Comput. 2012, 8, 3839-3853