ARTICLE pubs.acs.org/JCTC
Polarizable Atomic Multipole-Based Molecular Mechanics for Organic Molecules Pengyu Ren,*,† Chuanjie Wu,‡ and Jay W. Ponder*,‡ † ‡
Department of Biomedical Engineering, The University of Texas at Austin, Austin, Texas 78712, United States Department of Chemistry and Department of Biochemistry & Molecular Biophysics, Washington University, St. Louis, Missouri 63130, United States
bS Supporting Information ABSTRACT: An empirical potential based on permanent atomic multipoles and atomic induced dipoles is reported for alkanes, alcohols, amines, sulfides, aldehydes, carboxylic acids, amides, aromatics, and other small organic molecules. Permanent atomic multipole moments through quadrupole moments have been derived from gas phase ab initio molecular orbital calculations. The van der Waals parameters are obtained by fitting to gas phase homodimer QM energies and structures, as well as experimental densities and heats of vaporization of neat liquids. As a validation, the hydrogen bonding energies and structures of gas phase heterodimers with water are evaluated using the resulting potential. For 32 homo- and heterodimers, the association energy agrees with ab initio results to within 0.4 kcal/mol. The RMS deviation of the hydrogen bond distance from QM optimized geometry is less than 0.06 Å. In addition, liquid self-diffusion and static dielectric constants computed from a molecular dynamics simulation are consistent with experimental values. The force field is also used to compute the solvation free energy of 27 compounds not included in the parametrization process, with a RMS error of 0.69 kcal/mol. The results obtained in this study suggest that the AMOEBA force field performs well across different environments and phases. The key algorithms involved in the electrostatic model and a protocol for developing parameters are detailed to facilitate extension to additional molecular systems.
’ INTRODUCTION Organic molecules are the basic constituents of biology and material science. Modeling studies involving organic compounds are widely used in many areas such as physical chemistry, biological structure and function, and nanotechnology. Progress in quantum chemistry and the availability of fast computers has empowered the routine study of small molecules with high levels of ab initio theory and large basis sets. However, first principles statistical thermodynamics sampling techniques are still not practical for use with most high-level QM methods. Thus, molecular modeling based on empirical potentials is widely used for theoretical inquiries into microscopic and macroscopic phenomena across chemistry and biology. Atom-based force field models such as MM3,1 AMBER,2 CHARMM,3 OPLS,4 and GROMOS5 have been developed for a wide range of organic compounds and biomacromolcules. These models describe electrostatic interactions with fixed point charges on atoms and treat van der Waals interactions via Lennard-Jones potentials or other simple functions. Numerous studies have shown that many of the physical properties and structures of organic molecules can be adequately reproduced with current fixed charge force fields. Increases in computing power have enabled the simulation of larger molecular systems and more precise investigation of their properties. However, there are acknowledged shortcomings of the current generation of fixed charge potentials. They assume the atomic charges derived from training systems are approximately transferable to systems in different chemical environments. Explicit accounting of many-body effects is required for a general potential to capture the electrostatic response to different r 2011 American Chemical Society
molecular environments: homo- or heterogeneous, low or high dielectric, nonpolar or highly polarizable. Polarization effects were initially used in the description of molecular refractivity and other chemical phenomena nearly 100 years ago.6 Early in the era of modern computational chemistry, polarization was applied to the study of enzymatic reactions7 and incorporated into prototype molecular dynamics algorithms.8 Recently, there have been increasing efforts toward developing polarizable force fields for molecular simulation, based on a variety of empirical models for induction, such as classical induced dipoles,2,922 fluctuating charges,2330 and Drude oscillators.9,3135 Detailed discussions of the various polarization models can be found in recent reviews of polarizable force field development.3640 The performance of different approaches in accounting for polarization has been compared in the study of ion and small molecule interactions.41,42 The modeling of neat organic liquids, including alcohols, acids, amides, and aromatics, has also been reported using polarizable potentials.11,22,35,4350 Restriction to fixed atomic point charges constrains the flexibility of a model in representing the electrostatic potential around a molecule51,52 and thus limits the accuracy of the treatment of molecular interactions. Improvement can be achieved by adding extra charge sites, typically at bond centers or lone pair positions. For example, the TIPxP series of water models, TIP3P,53 TIP4P,53 and TIP5P,54 adopts increasing numbers of charge sites. Recently, the extra site approach was introduced into a Drude Received: May 2, 2011 Published: August 31, 2011 3143
dx.doi.org/10.1021/ct200304d | J. Chem. Theory Comput. 2011, 7, 3143–3161
Journal of Chemical Theory and Computation
ARTICLE
oscillator-based polarizable model as a way to address the anisotropy in atomic charge distribution due to lone pair electrons.50 Alternatively, one can directly incorporate higher order moments, such as dipole and quadrupole moments, at the atomic centers to improve the representation of the charge distribution. The convergence advantage of using multipoles distributed over atomic sites, as opposed to a single molecule-centered set of moments, has been discussed in the literature.55,56 Over two decades ago, Buckingham and Fowler57,58 were the first to apply distributed multipole moments to structural modeling of small molecule complexes. Their proposed intermolecular potential, consisting of hard sphere repulsion and atomic multipole-based electrostatics, was able to reproduce a number of experimental equilibrium geometries and orientational preferences. Recently, coarse-grained potentials with point multipoles have been used successfully in modeling hydrogen-bonded molecular liquids.59,60 The AMOEBA (Atomic Multipole Optimized Energetics for Biomolecular Applications) force field was initially developed for water.18,20 The current study reports the extension of the AMOEBA model to organic compounds including alkanes, alcohols, amines, sulfides, aldehydes, carboxylic acids, amides, and aromatics. A cornerstone of the AMOEBA force field is an improved electrostatic potential based on atomic multipoles and classical induced dipole moments. The atomic multipole moments are obtained from high-level ab initio calculations on gas phase monomers. An empirical atomic dipole induction model describes the many-body polarization effects important in clusters and condensed phase environments. A small, consistent set of atomic polarizability parameters is used to treat intermolecular polarization as well as intramolecular polarization between functional group fragments. van der Waals (vdW) parameters are refined via gas-phase homodimer molecular orbital calculations and molecular dynamics simulation of liquid properties. Additional gas-phase and liquid-phase computations, including hydrogen bonding in gas-phase heterodimers, the dielectric and diffusion constants of neat liquids, and hydration free energies of organic compounds have been utilized to validate the resulting force field.
’ METHODS Potential Energy Model. The interaction energy among atoms is expressed as
U ¼ Ubond þ Uangle þ Uba þ Uoop þ Utorsion perm
þ UvdW þ Uele
ind þ Uele
ð1Þ
where the first five terms describe the short-range valence interactions: bond stretching, angle bending, bondangle cross term, out-of-plane bending, and torsional rotation. The last three terms are the nonbonded interactions: van der Waals, permanent electrostatic, and induced electrostatic contributions. The individual terms for these interactions have been described in detail in a previous publication.61 Some additional methodology, introduced to treat electrostatic polarization in molecular systems beyond water, will be detailed below. Polarization effects in AMOEBA are treated via Thole’s interactive induction model that utilizes distributed atomic polarizability.62,63 According to this interactive induction scheme, induced dipoles produced at the atomic centers mutually polarize all other sites. A damping function is used at short-range to eliminate the polarization
catastrophe and results in correct anisotropy of the molecular response (i.e., diagonal components of the molecular polarizability tensor) starting from isotropic atomic polarizabilities. Thole damping is achieved by screening of pairwise atomic multipole interactions and is equivalent to replacing a point multipole moment with a smeared charge distribution.13 The damping function for charges is given by F¼
3a expð au3 Þ 4π
ð2Þ
where u = rij/(αiαj)1/6 is the effective distance as a function of interatomic distance rij and the atomic polarizabilities of atoms i (αi) and j (αj). The coefficient a is the dimensionless width of the smeared charge distribution and controls the damping strength. The corresponding damping functions for charge, dipole, and quadrupole interactions were reported previously.18 The Thole model is able to reproduce the molecular polarizability tensors of numerous small molecules with reasonable accuracy using only element-based isotropic atomic polarizabilities and a single value for the damping factor.62 In our water study, it was discovered that the dependence of molecular polarizability on the damping coefficient is weak, but the polarization energy is much more sensitive to the strength of damping. After fitting the interaction energies of a series of small water clusters, we have chosen a universal damping factor of a = 0.39, rather than the value of 0.572 suggested by Thole. We adopt the atomic polarizabilities (Å3) as originally given by Thole, i.e., 1.334 for carbon, 0.496 for hydrogen, 1.073 for nitrogen, and 0.837 for oxygen. The only exception is for aromatic carbon and hydrogen atoms, where we found that the use of larger values greatly improves the molecular polarizability tensor of benzene and polycyclic aromatics. The AMOEBA values for atomic polarizability are given in Table 1. In addition, for metal dications, we have found it necessary to use stronger damping (a < 0.39) to better represent the electric field around the ions.21,61,64 Intramolecular Polarization. For a large molecule such as a multifunctional organic or a biopolymer, polarization arises not only from the electric field of other molecules but also from distal portions of the same molecule. It is crucial to describe the intraand the intermolecular response in a consistent manner. In prior work, we investigated the effect of intramolecular polarization on the conformational dependence of the electrostatic potential surrounding a dipeptide.15 As observed by others,65 the electrostatic parameters derived for alanine dipeptide vary significantly depending upon the conformation used to derive the values. A simple average over multipole moments obtained from a set of conformers does not transfer well between conformers, i.e., gives poor electrostatic potentials on different conformers. Furthermore, when short-range polarization between bonded atoms is ignored, use of intramolecular polarization yields only marginal improvement over current nonpolarizable potentials. To overcome this problem, a group based intramolecular polarization scheme has been devised.15 The “groups” are typically functional groups with limited conformational degrees of freedom, such as an amide group or phenyl ring. In this scheme, the permanent atomic multipoles (PAM) polarize between, and not within, groups. For a small molecule consisting of a single polarization group, such as water or ammonia, permanent atomic multipoles do not polarize sites within the same molecule, while “mutual” induction occurs among all polarizable sites as described above. This design offers a clear connection between the treatment of small molecules and that of the analogous fragments inside a larger 3144
dx.doi.org/10.1021/ct200304d |J. Chem. Theory Comput. 2011, 7, 3143–3161
Journal of Chemical Theory and Computation
ARTICLE
Table 1. vdW Parameters and Atomic Polarizabilities for AMOEBA Atom Classes atom
ε (kcal/mol)
R0 (Å)
description
polarizability (Å3)
C
alkane (CH3 or CH2)
3.820
0.101
1.334
H
alkane (CH3)
2.960
0.024 (0.92)
0.496
H
alkane (CH2)
2.980
0.024 (0.94)
0.496
C
alkane (CH
