522
J. Phys. Chem. B 2011, 115, 522–531
Extending Classical Molecular Theory with Polarization Tom Keyes* and Raeanne L. Napoleon Department of Chemistry, Boston UniVersity, Boston, Massachusetts 02215, United States ReceiVed: June 17, 2010; ReVised Manuscript ReceiVed: NoVember 15, 2010
A classical, polarizable, electrostatic theory of short-ranged atom-atom interactions, incorporating the smeared nature of atomic partial charges, is presented. Detailed models are constructed for CO monomer and for CO interacting with an iron atom, as a first step toward heme proteins. A good representation is obtained of the bond-length-dependent dipole of CO monomer from fitting at the equilibrium distance only. Essential features of the binding of CO to myoglobin (Mb) and model heme compounds, including the binding energy, the position of the minimum in the Fe-C potential, the Fe-C frequency, the bending energy, the linear geometry of FeCO, and the increase of the Stark tuning rate and IR intensity, are obtained, suggesting that a substantial part of the Fe-CO interaction consists of a classical, noncovalent, “electrostatic bond ”. The binding energy is primarily polarization energy, and the polarization energy of an OH pair in water is shown to be comparable to the experimental hydrogen bond energy. I. Introduction The subject of this paper is the application of classical electrostatics to the calculation of molecular electric multipoles, vibrational spectroscopy, molecular geometry, and some aspects of chemical bonding: a theory valid on the angstrom length scale is required. Of course, one can simply perform a quantal calculation. However, our aim is to explore how far classical electrostatics may be extended into what is normally considered the realm of quantum mechanics. The general chemistry procedure of assigning fixed partial point charges to atoms is the basis of ionic bonding, in which the Coulomb energy at the equilibrium geometry approximates the binding energy (BE). It also allows the calculation of molecular electric multipoles. A glaring defect of ionic bonding is that it cannot bind atoms with like partial charges. Also, the Coulomb potential cannot exhibit a minimum corresponding to an equilibrium molecular geometry. The geometry arises from the combination of the Coulomb attraction and repulsive forces. A hint of what is missing in the partial charge approach may be found in CO, where the dipole points toward O, requiring an unphysical positive charge on oxygen. From the perspective of distributed multipole analysis,1 molecular multipoles are built up from atomic multipoles derived from charge distributions assigned to the atoms. Using partial charges only ignores the asymmetry of the distributions. Including atom dipoles, the larger on C pointing toward O and the smaller on O pointing toward C, yields2 the correct molecular dipole, even with a negative charge on O. More complete electrostatics provides a significant improvement over point charge models. However, we will argue that, for short-ranged interactions, a zero-field theory is inadequate. Molecular electrostatic properties respond to the environment through polarization by the local electric field. We believe that the associated polarization energy, varying as the square of the field, is important for chemical bonding. Fields are large at bonding pair separations. While polarization is broadly recognized as, e.g., causing the increase of the mean molecular dipole of water upon condensation from 1.85 to 2.6-3.1 D,3 very-short-ranged contributions
to the polarization energy, notably those between bonded atoms, are often ignored. Two reasons are that a bond is invariably already described by a bond potential, and the presence of the “polarization catastrophe”, a singularity in the induced dipole at a critical atom-atom distance. However, (1) a bond potential is either a zero-field energy from a quantal calculation or an empirical effective potential for a specific environment, while bonded atoms in any interesting situation are always in fluctuating local fields, which change their interaction energy, and (2) recognition of the nonpoint or “smeared” nature of atomic charge distributions damps the interactions, as originally described by Thole,4 averts the catastrophe, and yields realistic energies at distances where charge clouds overlap and the multipole expansion breaks down. The Thole approach may be justified in two ways. First, electron distributions are indeed “smeared”, and the properties of charge distributions are quite different from those of point charges or point dipoles at short range. The result is a proper expression of classical electrostatics for the case of overlapping charge distributions. Second, if charge distributions overlap, exchange, which is purely quantal, is also important. Within a classical approach, the smeared charge functions provide an empirical fit to the exchange energy. Thus, we are extending the usual picture both by providing a version of classical electrostatics that is valid at short range, and by empirically including related quantal effects. The classical formalism remains unchanged, so we continue to refer to our approach as such. Short-ranged polarization has probably been most carefully studied in water. The most successful water potentials use Thole damping, or variants.5-12 Among other things, they predict cluster energies which are in excellent agreement with quantal calculations, expressing the many-body energy via polarization. All of these potentials include intramolecular effects, but, except for our POLIR,12 they lack the contribution to the polarization energy from bonded atoms. While the induced dipole on an atom gets a contribution from the field from the induced dipole on a bonded neighbor, there is no term in the polarization energy for its interaction with the field from the partial charge on the bonded neighbor. This approach prevents
10.1021/jp105595q 2011 American Chemical Society Published on Web 12/23/2010
Classical Molecular Theory with Polarization overcounting of interactions already present in the bond potential, but in doing so discards the very real field dependence of the bond potential. Our thesis is that short-ranged classical polarization, including all interactions of bonded atoms, is responsible for several phenomena generally considered to be quantum mechanical. Along with everything else, quantum mechanics provides a complete description of polarization. We believe that sometimes the theory of polarization is the most important contribution of the quantal formalism, and it may be implemented classically. For specific examples, consider that in our first two papers on this topic2,13 we obtained a classical description of the vibrational Stark effect (VSE, response of a vibrational frequency to an electric field) of CO monomer and the vibrational IR spectrum of CO liquid. Turning to ligand-heme protein interactions, we did not attempt to construct a complete potential. The idea was that including polarization on the iron and on the ligand would capture some effects of ligation, with the geometry taken from experiment, as in ionic bonding. The resulting theory13 yielded more nominally quantal results: noncovalent bonding of CO, NO, and O2 to iron, the reduction of the binding energy of carbonmonoxy myoglobin, MbCO, due to interaction with the rest of the protein, and the ≈3.5× increase of the VSE14-16 of CO upon binding to the iron in Mb. We proposed that classical, polarization-energy-dominated “electrostatic bonds” could describe interactions intermediate in character between ionic and covalent. The polarization energy can bind atoms with like partial charges. Both C and Fe have positive partial charges. They nevertheless form an electrostatic bond because the polarization energy produces an attraction that is stronger than the Coulomb repulsion for the appropriate distances. The energy of FeOC was found to be ≈10 kcal/mol higher than that of FeCO, because O is less polarizable than C. An Fe-C attraction can also be found with complete nonpolarizable electrostatics, not limited to point charges. The question then is simply the magnitude of the different contributions. For our primary results in the following, the polarization energy is dominant. In addition to any model-dependent findings, the action of polarization is indicated by spectroscopic changes. The dipole correlation determines the IR spectrum, polarization governs the response of the dipole to the environment, and a large increase in absorption intensity is a clear signal of the presence of polarization. Ahlborn et al.17 showed that the ≈18× increase of the OH stretching intensity of ambient water compared to the monomer arose from induced dipoles, i.e., from polarization. Note also that one mechanism for the enormous intensity enhancements observed in SERS is18 the mutual polarization of the molecule being studied and a metal surface. In quantal calculations of the spectrum, the environment is included through“non-Condon effects”,19 in which the transition dipole depends upon the coordinates of the neighbor molecules. The transition dipole is proportional to the dipole derivative and has been shown19 to correlate very well with the local electric field. With a dipole derivative responding to a local field, it would seem that non-Condon effects are capturing the effect of polarization, although the precise relationship remains to be formulated. When CO binds to Fe, in addition to the enhancement of the VSE mentioned above, the IR intensity of the CO stretch increases by ≈34×.20 From our unified perspective on spectroscopy and energetics, this indicates that the focus should be
J. Phys. Chem. B, Vol. 115, No. 3, 2011 523 on the polarization energy. We emphasize that a theory should predict spectra as well as the more usual proberties. The spectrum is not a peripheral matter, but an essential test of whether polarization is properly described. Conversely, optimizing an empirical potential for the condensed-phase spectrum is an excellent procedure for fine tuning the effects of polarization. An open question21 is whether classical mechanics is capable of describing condensed-phase vibrational frequency shifts. A hybrid alternative is to perform a classical simulation, but to use the VSE21-25 to estimate frequencies from the local electric field. Given that we found the VSE to be dominated by polarization, and given the findings17 for the intensity in liquid water, we considered that shifts might be dominated by polarization forces and constructed the POLIR water potential12 accordingly, including the polarization energy of bonded neighbors. Straightforward classical simulation with POLIR, and with other recent polarizable potentials,6,7 yields excellent vibrational frequencies and absolute intensities of clusters, the liquid, and ice, answering the question in the affirmative. In this paper, our ideas about short-ranged polarization are further developed. A detailed electrostatic description of CO is obtained, including the bond-length-dependent dipole, normally computed with quantum mechanics.26 With bonded atoms polarizing each other, a change in the bond length is a change in the environment, and the atom dipoles respond accordingly. Using an FeCO model to represent the most significant polarizable atoms on MbCO or (model heme)-CO compounds, we obtain binding and bending energies, and the linear geometry, in agreement with quantal calculations,27,28 and the FeC frequency15 and the increases of the VSE14-16 and the IR intensity20 upon binding, in agreement with experiment. In sharp contrast to the Coulomb energy, the electrostatic potential possesses a minimum reflecting the molecular geometry, as a function of both θ, the deviation from linearity, and RFeC. Polarization energy is the largest part of the total energy near the minimum, i.e., of the binding energy. The minimum vs RFeC arises because, as the charge distributions begin to overlap, damping causes the energy to rise, in an apparent repulsion. The electrostatic minimum should lie inside the true minimum and be pushed out by the van der Waals repulsions, which we do not include. Nevertheless, if we match the expected location and depth of the true well with electrostatics, the calculated Fe-C stretching frequency of 505 cm-1 is in excellent agreement with experiment.15 With an empirical model it is essential to avoid mere parametrization. Note the following in the work to be presented below: (1) the physical mechanism of polarization by charge distributions determines the form of all the functions to be fit; (2) the parameters obtained have physically reasonable values; (3) in FeCO, the bending energy, the linear geometry, the IR intensity, and the Fe-C stretching frequency are not fit and represent very successful predictions; (4) in CO, the direction of the VSE and the r dependence of the dipole are predictions. We will emphasize the “value added” by the theory throughout. We submit that reproducing and predicting the many physical properties which we will demonstrate below makes a strong case for the importance of polarization in ligand-protein binding, even though our calculations are for FeCO only. Similar analysis of the POLIR potential argues that hydrogen bonds are electrostatic bonds. Thus, electrostatic bonding offers the prospect of a classical, computationally undemanding description of a significant part of an expanded class of systems, including molecular recognition, hydrogen bonding, and spectroscopy.
524
J. Phys. Chem. B, Vol. 115, No. 3, 2011
Keyes and Napoleon
II. Molecular Electrostatic Model 1. Induced Dipoles and the Electrostatic Energy. The atoms are assigned partial charges, qi, which may depend on the configuration, permanent dipoles µpi , and isotropic polarizabilities, Ri. The total dipole is the sum of the induced dipoles, permanent dipoles, and the “charge dipole”, calculated as usual
µ˜ ) µ˜
ind
+ µ˜ + µ˜ p
q
(1)
where the tilde indicates a quantity in a 3N-dimensional space, indexed by atom label i and Cartesian direction; µq(ix) ) xiqi. Summing the 3-dimensional pieces of µ˜ corresponding to each atom gives the 3-dimensional system dipole, µ. The induced dipoles obey
µ˜ ind ) (I˜ - r˜ · T˜)-1 · r˜ · (E˜int + E˜ext)
(2)
The atomic polarizability matrix, r˜ , is diagonal, with Ri in all three positions corresponding to atom i. The dipole tensor matrix, T˜, is the 3 × 3 dipole tensor for atoms i and j in the corresponding 3 × 3 block, and zero in the diagonal blocks. The external field is denoted E˜ext, the internal field, from all partial charges and permanent dipoles, but not induced dipoles, is E˜int, and E˜(ix) ) Ex(ri), where E(r) is the usual field. The quantity, r˜ sy ≡ (I˜ - r˜ · T˜)-1 · r˜ , is the 3N -dimensional system polarizability, with R˜ sy(ix,jy) giving the induced dipole on atom i in the x direction from a unit field applied to atom j in the y direction, and the sum over all its 3 × 3 blocks is the 3-dimensional system polarizability tensor, rsy. The electrostatic energy is equal to the energy in the absence of polarization plus the polarization energy, given by -1/2 × the sum of all induced dipoles dotted into the local fields, excluding fields from induced dipoles
1 Uel ) Ucp - (µ˜ p + µ˜ c) · E˜ext - (µ˜ ind,i + µ˜ ind,x) · (E˜int + 2 E˜ext) (3) where Ucp is the usual interaction energy energy of a collection of charges and permanent dipoles, and µ˜ ind,i and µ˜ ind,x denote the dipoles induced by E˜int and E˜ext, respectively, according to eq 2. The zero-external-field electrostatic energy is the sum of Ucp and the zero-field polarization energy
1 Upol,0 ) - µ˜ ind,i · E˜int 2
with the arrow indicating the usual 3-dimensional expression for a slowly varying field. The linear term is more interesting. What happens to the textbook -µ · Eext in the presence of internal polarization? We have
1 ˜ ext + µ˜ ind,x · E˜int) Uel,1 ) -(µ˜ p + µ˜ c) · E˜ext - (µ˜ ind,i · E 2 (7) The polarization contribution is rewritten using the molecular polarizability
˜ ext + µ˜ ind,x · E ˜ int ) E˜int · r˜ sy · E ˜ ext + E˜ext · r˜ sy · E˜int ) µ˜ ind,i · E 2E˜ext · µ˜ ind,0 (8) where we used the symmetry of rsy, and the zero-field induced dipole is
˜ int µ˜ ind,0 ) r˜ sy · E
(9)
The upshot is
˜ ext f -µ0 · Eext Uel,1 ) -(µ˜ p + µ˜ c + µ˜ ind,0) · E˜ext ) -µ˜ 0 · E (10) with the arrow indicating constant electric field. Thus, the usual result is recovered, even when the zero-field total dipole has contributions from internal polarization. However, the decomposition of the energy into polarization and nonpolarization contributions will differ, a point which is relevant to the discussion of the Fe-C-O bending energy below. 2. Smeared Charges. Following Burnham,7 short-ranged damping due to smeared charges is incorporated by using the following expressions in evaluating the above quantities: Coulomb potential from a unit positive charge
φ(r) ) [f2((r/γl);b) + (r/γl)Γ((b - 1)/b)Γ((b 1)/b, (r/γl)b)]
1 r
(11)
Coulomb field from a unit positive charge
(4) E(r) ) f2((r/γl);b)
The energy that is quadratic in the field is strictly due to polarization
rˆ r2
(12)
and the dipole tensor el,2
U
)U
pol,2
1 ) - µ˜ ind,x · E˜ext 2
(5) T(r) ) [3rˆrˆf1((r/γl);b) - If2((r/γl);b)]
and expressing the induced dipoles with the polarizability
1 1 Uel,2 ) - E˜ext · r˜ sy · E˜ext f - Eext · rsy · Eext 2 2
(6)
1 r3
(13)
Reduced units of angstrom and proton charge, with factors of 4πε0 set equal to unity, are employed in the calculations, and converted to experimental units as necessary.
Classical Molecular Theory with Polarization
J. Phys. Chem. B, Vol. 115, No. 3, 2011 525 el,0 el UMXY ) B(r) - UXY (r) + UMXY (RMX, r)
For an interaction between the ij atom pair
l ) (RiRj)1/6
(14)
The damping functions are
f1(x, b) ) 1 - (1 + xb)e-x
b
(15)
and
f2(x, b) ) 1 - e-x
b
(16)
and Γ with one and two arguments represents the gamma function of the first and second kind, respectively. The damping parameters are γ and b; γ scales the natural distance, l, for the onset of charge cloud overlap, and b governs the sharpness of its falloff with increasing r. The parameters are different for each atom pair. Burnham has considered models with b ) 3 and b ) 4; we regard both γ and b as adjustable. According to eqs 12 and 16, the electric field varies as rb-2 for small r. The field from a charge distribution ∼rp varies as rp+1 inside the charge-containing region, and b ) p + 3. The case of b ) 3 thus corresponds to p ) 0, a uniform charge distribution. Larger values of b, as utilized in refs 7 and 12 and in the following, indicate a shell-like distribution which vanishes at the origin. 3. Adding a Bond Potential. General Considerations and Diatomic Case. To include inherently quantal effects, our approach is to combine classical electrostatics with an accurate, zero-field bond potential, B(r˜), for all the covalently bonded atoms. The zero-field bond electrostatic energy is contained in the bond potential, so the total energy obeys
U(r˜) ) B(r˜) - S0(r˜) + Uel(r˜)
(17)
where S0 corrects for double counting. If V is a bare potential containing only the effects absent from our model, then U ) V + Uel, and V ) B - S0. For a single diatomic, the bond potential provides a complete description at zero field, and S0 ) Uel,0
U(r) ) B(r) - Uel,0(r) + Uel(r)
(18)
Then, our contribution in this simplest case is the field dependence, encompassing both an external field and the local fields from the neighbors in condensed phases. MXY Complex. A hint of the richness of polarizable molecular electrostatics is obtained by adding just one more atom. Consider a complex of an atom M, e.g., a metal in a protein, with a diatomic ligand, XY, having bond potential B(r), where r ≡ RXY. No bond potential is used for the M-X interaction although at least a van der Waals repulsion is required for a complete description. For large RMX, the energy equals B(r), and the bare potential el,0 (r), as for the diatomic. If V does not change equals B(r) - UXY upon binding
(19)
el allows for a nonlinear complex. where writing r in UMXY In contrast to the diatomic, the bond potential does not contain all the zero-field electrostatic energy. Some of it goes into producing an effective r-dependent bond potential, given by e ) for the bound complex, with a minimum position, U(r;rˆ,RMX ec r , shifted from the value in the monomer, re, and a shifted classical frequency. The remainder of the zero-field electrostatic energy becomes the BE, the difference in energy of the unbound and bound complexes
el,0 ec el,0 e BE ) (B(re) - B(rec)) + (UXY (r ) - UMXY (RMX , rec)) (20)
If binding changes the bare potential, we will simply introduce a bound bond potential, with a modified minimum position, denoted reb, which yields an optimal bound effective XY potential when used in eq 19. 4. Vibrational Stark Effect. Now we outline the classical VSE. The results below were already in the literature;29,30 our contribution is implementation with a detailed classical dipole and polarizability model, dominated by polarization. For simplicity, we consider a diatomic, but the results are applicable to any normal coordinate, r, of a larger system. For one molecule in an external field, E, applied parallel to the dipole, and hence the bond
1 U(r) ) B(r) - µ(r)E - R|(r)E2 2
(21)
The bond potential has a minimum at re, and the squared classical vibrational frequency is B′′(re)/m, where m is the reduced mass. The applied field changes the position of the minimum to re(E), such that U′(re(E)) ) 0, and changes the force constant. Denoting the change in the minimum position, re(E) - re as δr(E), and expanding everything around re, the first- and secondorder coefficients in E may be obtained
δr(1) )
δr(2) )
µ′ B''
(22)
µ′ 1 1 µ′′µ′ 1 R′ + - B′′′ B'' 2 | B'' 2 B''
(
2
( ))
(23)
where the derivatives with no argument are to be evaluated at re. The energy at the minimum is
U(re(E)) ) B(re) - µ(re)E -
(
)
1 (µ2 2 R|(re) + E 2 B''
(24)
with the rightmost term giving the effect of the change in position of the minimum. Evaluating U′′ by expanding around re and using the above, we obtain the force constant in the presence of the field
526
J. Phys. Chem. B, Vol. 115, No. 3, 2011
U′′(re(E)) ) B'' +
Keyes and Napoleon
- µ )E + O(E ) ( µ′B′′′ B''
2
(25)
which yields the classical ω(E). A quantum anharmonicity correction31 must be applied to ω(E), because classical mechanics sees the harmonic frequencies at the minima of the Born-Oppenheimer surface, which are not identical to the quantal transition frequencies. To first order in E the zero-field correction still applies, and the classical frequency shift equals the quantal shift. The Stark tuning rate (STR), ∆µ, is defined as the linear term in the redshift, and we obtain
ω(0) µ′B′′′ ∆µ ) µ′′ 2B'' B''
(
)
(26)
The notation “∆µ” came about because the difference in transition dipoles between the ground and excited states gives one contribution to the quantal STR, but ∆µ is now commonly used to represent the entire experimental quantity. The potential of a diatomic softens as the bond length is stretched beyond the minimum, corresponding to negative B′′′. If the dipole were represented with point charges only, µ′ ) q and µ′′ ) 0, where q is the partial charge, and the result is positive STR, a red shift. Polarization both causes µ′′ to be nonzero and influences the value of µ′. It appears21-25 that the VSE, responding to the local field, describes a substantial part of the IR frequency shifts from the monomer observed in condensed phases. Given the role of polarization in µ′ and µ′′, and given that the quadratic contribution to the classical frequency shift
δω(2) 1 )[2R′′|(B′′)3 - 2R′|(B′′)2(B′′′) + ω(0) 8(B′′)4 3(µ′B′′′)2 - 2(µ′)2B''B′′′′ - 6µ′µ′′B''B′′′ + (µ′′B′′)2 + 4µ′µ′′′(B′′)2]E2 (27) contains the polarizability derivatives, it is clear that a careful treatment of polarization is a necessary condition for the classical prediction of the STR and of classical vibrational frequencies. In particular, since polarization determines the response to the environment, it is essential for any large changes upon condensation or binding. The applicability of classical methods cannot be evaluated without regard to this condition. III. Carbon Monoxide We fit our original CO model2,13 to a selection of experimental data, plus total atom dipoles calculated with Gaussian 03,32 and used an early version of Thole damping. Here we construct an improved model, seeking to reproduce more molecular properties, and incorporating Burnham’s latest damping functions.7 We proceed by minimizing the summed squared deviation of the calculated and desired values. The parameters to be varied are the following: the carbon partial charge, qe (qO ) -qC), at the equilibrium bond length; two atom polarizabilities, RC and RO; two atom permanent dipoles, µCp and µOp; two damping parameters, γ and b; γq, which governs the bond-length dependence of the charge, q(r) ) qe exp[-γq(r - re)]. Regarding γq, the motivation was to see if classical electrostatics could predict the full µ(r), not just fit µ(re), and fixed charges will simply lead to a dipole that grows with increasing
TABLE 1: Parameters of the Molecular Electrostatic Model for CO Monomera qe RC RO µCp µOp γ b γq
0.262 1.79 0.0995 0.143 -0.110 1.32 5.60 1.21
a Units are Å3 for polarizabilities and e Å for dipoles; damping constants are dimensionless.
r, failing to describe dissociation to neutral atoms. The atom dipoles depend upon r through intramolecular polarization. The initial goal was to reproduce, dipole and first two derivatives, µ(re), µ′(re), and µ′′re); molecular parallel and perpendicular polarizabilties, R| and R⊥; polarizability derivatives, R′|(re and R′⊥(re)); and principal quadrupole, Qzz(re). Fitted properties were experimental except for the polarizability derivatives, computed with Gaussian 03,32 and were chosen as follows. We seek a CO model which will respond correctly to the local environment. An accurate dipole surface is crucial33 for the response of a water molecule to its environment, and there is no reason to believe that the situation is any different for CO. Having the dipole and its first two derivatives correct at re can be considered as first steps to building up the complete zerofield, nonlinear dipole surface, while the polarizability and its derivative describe the field dependence. The need for the derivatives in obtaining frequency shifts is explicitly seen in the previous section on the VSE. We estimate the “experimental” µ′′ from the experimental STR15 and eq 26; that is, the STR is built in, so the vibrational frequency will respond properly to local fields.14-16,21-25 The VSE is the only monomer property requiring the bond potential. We employ the34 Huffaker potential, as in our prior2 work. The electric quadrupole is important for the structure of the liquid,35 and for determining the preferred configurations of CO in the heme pocket.36-38 Additionally, µ′ and R′ are minimal requirements for a spectroscopic theory, as they determine the IR and Raman intensities. Again, we argue that electrostatics, energetics, and spectroscopy should be treated in a unified fashion. It was impossible to reproduce R′⊥(re) while keeping good values of the other properties. Since R′⊥(re) is the least important quantity, entering neither the STR nor the MbCO calculation for the prevailing linear geometry, we abandoned it and attained a representation of the remaining properties, detailed in Tables 1 and 2, which is correct to within experimental and computational error. Over and above the conceptual framework and the perfect representation of the selected properties, what is gained from the model? Note that all our fitting has been at the equilibrium bond length. Figure 1 shows the complete r-dependent dipole, along with representative points estimated from Figure 1 of Harrison’s26 paper (length unit changed to Bohr a0 for the comparison). Note that the quantal calculations are not exact, overestimating µ(re) by ≈2×. We obtain qualitatively reasonable agreement with the quantal calculations for all r from fitting at re only. The characteristic behavior of the dipole near re follows from intramolecular classical polarization. For r slightly greater than re the dipole points toward C, the expectation from partial charges. However, carbon is much more polarizable than
Classical Molecular Theory with Polarization TABLE 2: Calculated CO Monomer Properties at the Equilibrium Bond Length, 1.128 Åa R| R⊥ R′| µ µ′ µ′′ ∆µ Qzz
2.3239 1.7939 2.90 0.022840 -0.67041 -2.75 -0.69015 -0.51242
a Carbon is at the origin of coordinates with oxygen on the +z axis. Units are constructed from the basic Å length unit and proton charge unit, except ∆µ is cm-1/(MV/cm). The achieved representation of the selected properties is essentially perfect, so only one number is shown.
Figure 1. Bond length (in bohr) dependent dipole (+ is toward O) of CO monomer from polarizable electrostatic model. The vertical line denotes equilibrium separation, and the stars are from the quantal calculation in ref 26.
oxygen, and the field from the negative charge on O, pointing toward O, induces an atom dipole on carbon which exceeds the charge dipole, as r f re and the field grows. Harrison26 also discusses competing charge and atom dipoles, but with the latter obtained from quantal calculations; we have shown that there exists a classical alternative. Furthermore, while the magnitude of the Stark tuning rate is fit to Boxer’s15 result, the negative sign is a prediction, as the model will not give a positive ∆µ and maintain the desired properties. So, a field pointing from O to C gives a red shift, as was argued14 by Park et al. Finally, we consider that the primary “value added” of the CO model is its ability to predict the response to the local field, e.g., a nearby iron, the topic of the next section. Many nonpolarizable models have been published43-46 with a focus on explaining experiments47,48 on CO in the unbound docking sites of the heme pocket. Several have yielded the experimental dipole, but not higher multipole moments.43-45 Straub and Karplus accurately reproduced the quadrupole moment of CO and various ab initio interaction energies with Mb.36 By fitting the r-dependence of the atomic charges, Meuwly and co-workers37,38 obtain a good representation of µ(r) and Q(r) for 0.8-1.5 Å and have had success in predicting the IR spectrum and preferred docking sites of unbound MbCO. Perhaps such fitting incorporates polarization empirically, since, as we demonstrate next, the response of polarizable CO to local fields is crucial to the Fe-CO interaction; notably, the binding energy is mostly polarization energy. IV. Heme-CO Binding We aim to establish that classical electrostatic bonding, dominated by polarization, plays a significant role in ligation. Note that we are referring to the protein-ligand bond itself,
J. Phys. Chem. B, Vol. 115, No. 3, 2011 527 not the influence of nearby residues. The heme binds CO, O2, NO, and CN-, and the relative affinities are of immense physiological consequence. In future work we will consider cobalt-containing proteins,49 which scavenge the impurities CNand H2S in crude oil. Thus, we will bring CO, obeying the model just constructed, into close proximity with an iron atom, with no Fe-C bond potential. The geometry and iron properties are chosen to be representative of MbCO and (model heme)-CO compounds. The result is a classical description of many aspects of CO binding and the associated spectroscopy. Perhaps most strikingly, the Fe-CO bond itself is obtained. It seems that the signature of polarization is strongly written on the heme-CO interaction. Since we are not treating any other atoms or other terms in the potential, and since we seek to broadly represent heme-CO compounds, a qualitative description is anticipated. The idea is to see how much can be obtained from polarizable electrostatics alone, and that Fe, C, and O are the most important polarizable atoms for phenomena related to ligand binding. Nevertheless, some of the resulting predictions are in remarkable agreement with experiments on MbCO and quantal calculations on heme-CO, suggesting that, if the observables are indeed dominated by polarization, the details are not so significant and our simple approach may yield quantitative results. Of course, a large body of work exists in this area. While is is routinely said that recognition of ligands by proteins is noncovalent, much theory of ligation is based upon quantum mechanics.50-52 Steric hindrance from specific residues is a primary contributor to selective recognition.50 More recently, it has been argued that the residues enable recognition of diatomics by generating specific electrostatic environments.51 Biophysical examples of ionic bonding have been presented, e.g., the case of zinc ions attracting partially negative atoms in a peptide.53 However, there is no ionic bond between partially positive iron and partially positive carbon. The development of detailed, all-atom polarizable potentials for biophysical applications is an active area,54-56 and includes metal complexes and ligation. In light of such efforts, it is important to highlight the contributions of our model, namely, that it (1) obtains nontrivial results with extreme simplicity; (2) is focused upon the broad goal of extending classical theory beyond the ionic bond, using as a “demonstration project ” a bond generally considered to be covalent; (3) includes the polarization energy of the induced dipoles on C and O in the fields from the permanent charges and dipoles on their bonded neighbors; (4) treats spectroscopy on an equal footing with energetics; and (5) seeks a qualitative physical understanding of the action of short-ranged polarization. 1. Properties of Heme-CO Compounds. The C-O bond length increases upon binding to Mb, with estimates of the new value as high as 1.20 Å in the literature.27,57 The range for the Fe-C bond is 1.73-1.92 Å.27 Similar results are found for model heme compounds.27 BE for model compounds are27 25-35 kcal/mol, which are the relevant values for this work, since experimental BE for MbCO also include the effects of bringing the ligand into the protein. The minimum energy configuration of Fe-C-O is linear, and the bending energy has been estimated as ≈1.0 kcal/mol for the first 10 degrees.27,28 Vibrations of the ligand are a useful probe of the binding interaction. Bound MbCO has been has been investigated with conventional and nonlinear IR spectroscopy,22,48,58 and vibrational Stark spectroscopy.14-16
528
J. Phys. Chem. B, Vol. 115, No. 3, 2011
The bound “A” state of MbCO has three absorption bands, designated A0, A1, and A3, at 1965, 1944, and 1930 cm-1, respectively,22,58 constituting a red shift of ≈200 cm-1 from the monomer. While much effort22,58 has been devoted to explaining the splittings of the A state in terms of interactions with nearby residues of the protein, such considerations are obviously beyond the scope of this paper, since we do not include those residues. The STR and the IR intensity are enhanced14-16 by ≈3.5× and ≈34×,20 respectively. 2. Polarizable FeCO Model: Parameters. For meaningful results, it is essential to minimize the number of adjustable parameters. Thus, we leave CO monomer unchanged, except for allowing binding to shift the minimum of the bond potential, re f reb. The new parameters to be determined are reb, qFe, RFe, p µFe , and, for the Fe-C damping, bFeC and γFeC. In the quest for simplicity, bFeC is held at the value found for CO, bFeC ) 5.6, p and Fe is not given a permanent dipole, µFe ) 0. These simplifications are consistent with the qualitative theory which we are constructing, in contrast to the CO monomer case. Some guidelines exist for qFe and RFe. The heme iron is believed to remain in the +2 state after binding CO, in contrast to the O2 case. However, the charge is spread over the porphyrin ring. Some recent estimates for the remaining Fe partial charge are 1.21 for ferrous and 1.33 for ferric,59 and 0.95 when binding O2.60 The CHARMM61,62 charge for bound CO is 0.24, which is an effective charge in the context of the rest of the CHARMM potential. We will be content with charges representative of this range. Studies on iron clusters63 indicate that RFe is in the vicinity of 9 Å3. However, we are unaware of information on the heme iron. In our previous paper,13 we found RFe ) 4.39 Å3. 3. Polarizable FeCO Model: Calculations. The proper procedure with an empirical model is to select some properties to fit, to determine any unknown parameters, and to measure success by the prediction of additional properties. The following calculations all build in the CHARMM MbCO bond lengths of RFeC ) 1.92 Å and RCO ) 1.17 Å, a representative BE of 0.1 reduced unit (1 reduced unit ) 332 kcal/mol),27 and the experimental MbCO STR of 3.5× that of CO monomer. In a further effort to be parameter-lean, we adopt our previous13 value of RFe ) 4.39 Å3, eliminating another possible adjustment of the theory. First, we follow the ionic bonding procedure of rigidly fixing RFe ) 1.92 Å plus the known Fe-C-O angle of 180°. However, the desired RCO ) 1.17 Å is taken from the minimum of the bond plus electrostatic potential. The parameters qFe, reb, and γFeC are then varied to obtain the indicated BE, RCO, and STR, leading to qFe ) 0.45, γFeC ) 1.041, and reb ) 1.195 Å. BE and the VSE. Again, we do not claim credit for quantities that are built in. It should be noted, however, that essentially perfect fits are easily obtained with reasonable parameter values. Existing calculations obtain the BE and the enhanced VSE from quantum mechanics,27,64 but it is eminently plausible that they are dominated by classical polarization, allowing greatly simplified computation in some instances. Infrared Intensity. Our first attempt at a prediction is the infrared intensity. Estimating the intensity as proportional to the square of the dipole derivative of the FeCO complex with respect to the CO normal coordinate, we obtain a 28× increase over the monomer. A recent quantal calculation65 obtained ≈9×. It is important to include the derivative of the induced dipole on Fe; considering the dipole of the CO fragment alone gives 6×.
Keyes and Napoleon
Figure 2. Energy cost for bending Fe-C-O (upper, heavy black). The intermediate lighter black lines, top to bottom, are the polarization, Coulomb, and charge (permanent atom dipole) contributions, and the stars are from the quantal calculation in ref 27. The red line represents the interaction of the Fe charge with the total atom dipoles.
We were not surprised to obtain a large enhancement, given our perspective on the role of polarization in spectroscopy, but the close agreement with experiment (34×20) is gratifying. The increase arises from induced dipoles, as in the OH stretch upon condensation of water,12,17 in the enhancement of the VSE, and in one mechanism18 of SERS. In bound CO, the spectrum is dominated by induced dipoles; conversely, the polarization energy must be a significant part of the binding energy. FeCO Angle and the Bending Energy. We will now show that our model exhibits a potential well, and does not require an imposed geometry. Thus, we next attempt to predict the dependence of the energy upon the angle, θ, measuring the deviation of FeCO from linear. Figure 2 (heavy black line) demonstrates good agreement with ab initio calculations, giving ≈0.8 kcal/mol for the first 10 degrees, as opposed to ≈1.0 kcal/mol.27,28 The minimum at θ ) 0 correctly27,28 predicts that the equilibrium geometry is linear. The energy minimum at θ ) 0 could not be found with the point charges alone. Then, bending the complex would decrease the energy, bringing the partially negative oxygen closer to the partially positive iron. By contrast, both polarization and the atom dipoles favor linearity, since the larger R| gives a more favorable polarization energy, and the much larger dipole on C pointing toward O aligns in the field from the positive Fe. It is then of interest to know the relative importance of the different effects. In Figure 2, the lighter black lines, top to bottom, show the polarization energy, the charge-atom dipole energy, and the Coulomb energy. The bending energy is almost entirely polarization energy. However, recall (section II.1) that, with internal polarization contributing to the total atom dipoles, the definition of polarization energy is modified. The red curve is the interaction of the Fe charge with the total atom dipoles. The charge-dipole energy is now substantially increased, but is still less than half of the total energy. We conclude that classical electrostatics can describe the bending potential quite well and that polarization energy is its most important component. Fe-C Bond and Stretching Frequency. An energy minimum is also found as a function of RFeC. Although the complete Fe-C interaction has a van der Waals term which we are not considering, as an illustration we now construct the Fe-C bond from damped electrostatics only, bearing in mind that exchange is included empirically in this fashion. Only a small change from the parameters of the previous section, to RFe ) 3.5, qFe ) 0.74, and γFeC ) 1.147, is needed to position the electrostatic minimum at RFeC ) 1.92 Å. The resulting potential energy surface, U(RFeC,RCO;θ)0), is shown in Figure 3. To two decimal places, the contributions to the BE of the polarization energy, the Coulomb energy, and the energy of interaction of the charge on Fe with the permanent dipoles on
Classical Molecular Theory with Polarization
Figure 3. Potential energy surface (reduced units, 1 unit ) 332 kcal/ mol) as a function of RCO and RFeC (Å).
Figure 4. Fe-C pair polarizability trace and anisotropy, from the electrostatic model, in units of (RFeRC)1/2.
CO are 0.12, -0.04, and 0.02 atomic units, respectively. In addition to the Coulomb repulsion, there is some attraction from going beyond point charges on CO, but on the whole the potential well is formed by polarization energy: Fe binds CO with an “electrostatic bond”. With RFeC obtained from a potential, we may make another prediction, the Fe-C stretching frequency. Calculating the harmonic normal modes of the linear complex, with Fe held fixed, we obtain ωFeC ) 505 cm-1, squarely in the center of the range of ≈475-525 cm-1 for 10 MbCO variants reported by Park and Boxer.15 Source of the Potential Well for RFeC. Why does the electrostatic energy exhibit a minimum as a function of RFeC? Examination of its separate components reveals that the minimum arises in the polarization energy, and upon further reflection, is perhaps not so surprising. For many years,66 calculations have been performed on the polarizability of pairs of atoms, primarily to describe light scattering from atomic fluids. The trace (with isolated atom contributions subtracted) and anisotropy of the pair polarizability are easily calculated with the usual dipole tensor. However, since the results become unphysically large at short, but accessible, pair separations, due to the proximity of the polarization catastrophe, quantum mechanics has been applied. With decreasing r, the tensor elements are seen to rise to a maximum, and then turn over.66 Figure 4 shows the FeC pair polarizability calculated with our model. Damped polarizable electrostatics reproduces the features found in quantal calculations at short distances, including the sign reversal of the trace. The minimum in the polarization energy is a simple consequence of the maximum in the pair polarizability, occurring approximately at the Fe-C bond length of 1.92 Å. Note that, despite the qualitatively correct form, we are not claiming to have accurately reproduced the true Fe-C pair polarizability. We are arguing that the well-known maximum in the pair polarizability makes it easy to understand why there should be a minimum in the polarization energy. Fitting to quantal calculations of the pair polarizability is a possible future direction for determining the parameters in the theory. For now,
J. Phys. Chem. B, Vol. 115, No. 3, 2011 529 the relevant observation is that the calculated maximum and minimum occur at the same RFeC. Summary. In this section we have sought to illuminate the role of polarization in ligation and the associated spectroscopy, while avoiding parametrization and carefully distinguishing between fit quantities and predictions. Beginning with the CHARMM geometry built in, and the BE and STR in accord with ab initio calculations and experiment, respectively, the IR intensity is our first prediction, and the results are very satisfactory. Allowing the FeCO angle to vary, it is next seen that we obtain the linear geometry, i.e., an energy minimum at θ ) 0, and a bending energy in good agreement with quantal calculations.27,28 Upon relaxing the constraint on RFeC and slightly adjusting the parameters so that the electrostatic energy minimum occurs at RFeC ) 1.92 Å with BE ) 0.1 reduced unit, the predicted Fe-C frequency is in perfect agreement with experiment. The bending energy and FeC frequency indicate that we obtain the correct shape of the potential well, not just the location of the minimum. We have just listed four predictions, which were not fit in any way, establishing the crucial point that we are getting more out of the theory than what we put in. The potential minimum vs RFeC arises because the polarization attraction of the like-charged atoms becomes a repulsion, due to damping, as the charge distributions begin to overlap. Concomitantly, the pair polarizability decreases66 at short range. Compared to the Coulomb potential, the richness of the damped polarization energy, and its possibility of describing chemical phenomena in a novel, computationally convenient fashion, are apparent. The red shift and bond lengthening of bound CO are not obtained. These phenomena are caused by backbonding, an inherently quantal mechanism. Backbonding changes the bond potential, and is incorporated in the theory by having reb > re, as discussed in section II.3. We emphasize that we are not claiming that polarization dominates all the significant changes upon ligation. Back-bonding leads to a net negative charge on the CO fragment. It might appear, then, that we should use a different electrostatic model for bound CO. However, we are trying very hard to minimize parametrization and to maximize simplicity, and so we have chosen not to do so. This is possible because we do not calculate the Fe charge and polarizability, and the Thole γFeC, from first principles, but use them as fit parameters. Thus, they take on effective values which mimic the effect of back-bonding on the Fe-CO interaction, and we are able to utilize a single CO model. Even without the red shift, the model is usable as is to obtain CO frequency fluctuations about a bound-state mean value, as is required22 for the calculation of nonlinear IR spectra, from local field fluctuations. We obtain the experimental STR at the experimental geometry, but the STR is in general a sensitive function of the geometry. It might be said that we have a fluctuating STR, while several existing theories21-25 use a fixed STR. It will be interesting to explore the consequences for the frequency correlation function and to compare our results with calculations based upon simpler frequency/field relations. V. Implications for Hydrogen Bonding in Water The POLIR water potential12 contains a comprehensive treatment of the polarization energy, along the lines described above. POLIR has performed very well overall, with excellent
530
J. Phys. Chem. B, Vol. 115, No. 3, 2011
Figure 5. O-H pair polarization energy (reduced units) vs ROH (Å) for POLIR potential.
results for spectra, cluster energies and dipoles, the dielectric constant, the second virial coefficient, and the diffusion constant. The pair distribution function is reasonable but is not as close to experiment as some other potentials. Note, however, that empirical water potentials are routinely optimized for g(r), and POLIR is not. The only intermolecular information built into POLIR is the minimum dimer energy vs ROO, and g(r) looks much better when seen as a true prediction. Also, other contemporary potentials5-7,9 have evolved through several versions, while POLIR is in “version 0”. Let us now ask if the previous analysis can provide some insights into POLIR’s success. Specifically, we ask if a hydrogen bond might be an electrostatic bond. Figure 5 shows the POLIR O-H pair polarization energy, exhibiting a minimum at ROH ∼ 1.6 Å. Liquid-state simulations with the full POLIR potential yield a most probable hydrogen bond distance of ≈2.0 Å, which is in agreement with other estimates.67 This distance lies outside the electrostatic minimum because of the van der Waals replusion. The possibly significant observation is that, at 2.0 Å, the pair polarization energy is very close to the actual hydrogen bond energy of 24 kJ/mol. The numbers could change in the next version of POLIR, but these are rough estimates and we do not expect the qualitative conclusion to change. Note also that POLIR’s strength so far is predicting polarization-related quantities and that the argument at hand does not depend upon g(r). Thus, we suggest that hydrogen bonds in water are electrostatic bonds, dominated by polarization energy. VI. Discussion In this paper we have implemented our ideas about shortranged classical electrostatics and polarization in CO and FeCO. Previously, we assumed that RFeC was large enough so that Fe-C damping could be ignored. Here we incorporate Burnham’s damping functions for both the Fe-C and C-O interactions. Our approach to the parameters in the theory is empirical, and we have emphasized that we then have the burden of predicting results that were not built in. While the CO model is parametrized at the equilibrium bond length, reasonable agreement with quantal calculations26 of the dipole is obtained over a broad range of atom separations. The unusual behavior of the sign of the dipole results from the competition of opposing charge and atom dipoles. The direction of the STR, ∆µ, is predicted correctly. In FeCO, the predicted IR intensity enhancement, FeCO bending energy, linear geometry, and FeC vibrational frequency are in good agreement with experiments and quantal calculations on MbCO and heme-CO compounds. The significance of polarization for condensed-phase vibrational spectroscopy cannot be overemphasized. The 34×
Keyes and Napoleon increase in IR intensity for bound CO and the 18× increase for liquid water are both due to induced dipoles interacting with the radiation field. While the model is fit to the 3.5× enhancement of the STR of bound CO, so it is not a prediction, the fact remains that the STR increase also stems from polarization. Ionic bonding via the electrostatic Coulomb energy is textbook material. However, the Coulomb interaction can neither bring together like-charged atoms nor, by itself, form a potential well. In ref 13 we suggested that one could describe cases intermediate between ionic and covalent via polarization, and named the results “electrostatic bonds”. In FeCO, the likecharged atoms actually attract more strongly than the unlikecharged atoms, making the point that classical electrostatic bonding is a natural extension of classical ionic bonding. Atomic multipoles higher than charges can also cause like-charged atoms to attract, but in our model this effect is much smaller than polarization, and it cannot cause a large increase in absorption intensity. Here it was shown that, in addition to generating an attraction, damped polarization can cause like-charged atoms to form a potential minimum, which attains a reasonable depth, location, and curvature, with plausible values of the parameters. With more development, protein-ligand binding might be treated classically, with no need to switch the potential between bound and unbound functions. Turning to water, we found that the OH pair polarization energy, evaluated with the POLIR parameters, has a minimum at 1.6 Å, and a BE of 24 kJ/mol at 2.0 Å. With the damped polarization energy for typical OH distances comparable to the experimental hydrogen-bond BE, it is plausible that hydrogen bonds are electrostatic bonds. Acknowledgment. Acknowledgment is made to the donors of the American Chemical Society Petroleum Research Fund and to the National Science Foundation (grant CHE 0848427) for support of this research. We acknowledge valuable discussions with Dr. Revati Kumar. References and Notes (1) Stone, A. J. The Theory of Intermolecular Forces; Oxford University Press: Oxford, UK, 1997. (2) Mankoo, P.; Keyes, T. J. Chem. Phys. 2006, 124, 204503. (3) Guillot, B. J. Mol. Liq. 2002, 101, 219. (4) Thole, B. T. Chem. Phys. 1981, 59, 341. (5) Fanourgakis, G.; Xantheas, S. J. Chem. Phys. 2006, 124, 174504. (6) Fanourgakis, G.; Xantheas, S. J. Chem. Phys. 2008, 128, 074506. (7) Burnham, C. J.; Anick, D.; Mankoo, P. K.; Reiter, G. F. J. Chem. Phys. 2008, 128, 154519. (8) Xantheas, S.; Burnham, C. J.; Harrison, R. J. J. Chem. Phys. 2002, 116, 1493. (9) Ren, P.; Ponder, J. J. Comput. Chem. 2002, 23, 1497. (10) Kumar, R.; Wang, F.; Jenness, G.; Jordan, K. D. J. Chem. Phys. 2010, 132, 014309. (11) Bernardo, D.; Ding, Y.; Krogh-Espersen, K.; Levy, R. J. Phys. Chem. 1994, 98, 4180. (12) Mankoo, P.; Keyes, T. J. Chem. Phys. 2008, 129, 034504. (13) Mankoo, P.; Keyes, T. J. Phys. Chem. B 2006, 110, 25074. (14) Park, E.; Andrews, S.; Hu, R.; Boxer, S. J. Phys. Chem. B 1999, 103, 9813. (15) Park, E.; Boxer, S. J. Phys. Chem. B 2002, 106, 5800. (16) Suydam, I.; Snow, C.; Pande, V.; Boxer, S. G. Science 2006, 313, 200. (17) Ahlborn, H.; Ji, X.; Space, B.; Moore, P. J. Chem. Phys. 1999, 111, 10622. (18) Efrima, S.; Metiu, H. J. Chern. Phys. 1979, 70, 1602. (19) (a) Schmidt, J.; Corcelli, S.; Skinner, J. L. J. Chem. Phys. 2005, 123, 044513. Corcelli, S.; Skinner, J. L. J. Phys. Chem. A 2005, 109, 6154. (b) Auer, B.; Kumar, R.; Schmidt, J.; Skinner, J. L. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 14215. (20) Lim, M.; Jackson, T.; Anfinrud, P. A. J. Chem. Phys. 1995, 102, 4355.
Classical Molecular Theory with Polarization (21) Li, S.; Schmidt, J.; Corcelli, S.; Lawrence, C.; Skinner, J. L. J. Chem. Phys. 2006, 124, 204110. (22) Merchant, K.; Xu, Q.; Thompson, D.; Fayer, M. D. J. Phys. Chem. A 2002, 106, 8839. Merchant, K.; Noid, W.; Thompson, D.; Akiyama, R.; Loring, R. F.; Fayer, M. D. J. Phys. Chem. B 2003, 107, 4. Finkelstein, I.; Goj, A.; McClain, B.; Massari, A.; Merchant, K.; Loring, R. F.; Fayer, M. D. J. Phys. Chem. B 2005, 109, 16959. Finkelstein, I.; Ishikawa, H.; Kim, S.; Massari, A.; Fayer, M. D. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 2637. (23) Kwac, K.; Cho, K. J. Chem. Phys. 2003, 119, 2247. (24) DeCamp, M.; DeFlores, L.; McCracken, J.; Tokmakoff, A.; Kwac, K.; Cho, M. J. Phys. Chem. B 2005, 109, 11016. (25) Hayashi, T.; Jansen, T.; Zhuang, W.; Mukamel, S. J. Phys. Chem. A 2005, 109, 64. (26) Harrison, J. J. Phys. Chem. A 2006, 110, 10848. (27) Rovira, C. J. Phys. Condens. Matter 2003, 15, S1809. (28) Rovira, C.; Kunc, K.; Hutter, J.; Ballone, P.; Parrinello, M. J. Phys. Chem. A 1997, 101, 8914. (29) Lambert, D. Phys. ReV. Lett. 1983, 50, 2106. (30) Marti, J.; Bishop, D. J. Chem. Phys. 1993, 99, 3860. (31) Berens, P.; Wilson, K. R. J. Chem. Phys. 1981, 74, 4872. (32) Frisch, M.; Trucks, G.; Schlegel, H.; Scuseria, G.; Rommet, M.; et al. Gaussian 03, ReVision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (33) Burnham, C. J.; Xantheas, S. J. Chem. Phys. 2002, 116, 1500. (34) Huffaker, J. J. Chem. Phys. 1976, 64, 3175. ; J. Chem. Phys. 1976, 64, 4564. (35) Geiger, L.; Ladanyi, B. M. J. Chem. Phys. 1988, 89, 6588. (36) Straub, J. E.; Karplus, M. Chem. Phys. 1991, 158, 221. (37) Nutt, D.; Meuwly, M. Biophys. J. 2003, 85, 3612. ; Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 5998. (38) Plattner, N.; Meuwly, M. Biophys. J. 2008, 107, 120519. (39) Bridge, N.; Buckingham, A. D. Proc. R. Soc. London A 1966, 295, 334. Voisin, C.; Cartier, A. J. Mol. Struct.: THEOCHEM 1993, 286, 35. (40) Muenter, J. J. Mol. Spectrosc. 1975, 55, 490. (41) Guadagnini, P.; Oliveira, A.; Bruns, R.; Neto, B. J. Am. Chem. Soc. 1997, 119, 4224. (42) Imrie, D.; Raab, R. Mol. Phys. 1001, 74, 833. Buckingham, A. D.; Longuet-Higgins, H. C. Mol. Phys. 1968, 14, 63. (43) Elber, R.; Karplus, M. J. Am. Chem. Soc. 1990, 112, 916. (44) Kottalam, J.; Case, D. J. Am. Chem. Soc. 1988, 110, 7690. (45) Henry, E.; Levitt, M.; Eaton, W. Proc. Natl. Acad. Sci. U.S.A. 1985, 82, 2034.
J. Phys. Chem. B, Vol. 115, No. 3, 2011 531 (46) Fracassi, P.; Della Valle, R. Chem. Phys. Lett. 1984, 104, 435. (47) Lim, M.; Jackson, T.; Anfinrud, P. A. J. Am. Chem. Soc. 2004, 126, 7946. (48) Helbing, J.; Nienhaus, K.; Nienhaus, G.; Hamm, P. J. Chem. Phys. 2005, 122, 124505. (49) Mosseri, S.; Neta, P.; Harriman, A.; Hambright, P. J. Inorg. Biochem 1993, 39, 93. Godinez-Salomon, F.; Hallen-Lopez, J.; Hopfl, H.; MorelesPacheco, A.; Beltran, H.; Zamuido-Rivera, L. Green Chem. 2005, 7, 716. (50) Collman, J. Acc. Chem. Res. 1977, 10, 265. (51) Olson, J.; Phillips, G. J. Bio. Inorg. Chem. 1997, 2, 544. (52) Schulze, B.; Evanseck, J. J. Am. Chem. Soc. 1999, 121, 6444. (53) Stote, R.; Karplus, M. Proteins: Struct., Funct., Genet. 1995, 23, 12. (54) Gresh, N.; Cisneros, G.; Darden, T.; Piquemal, J.-P. J. Chem. Theory Comput. 2007, 3, 1960. de Courcy, B.; Piquemal, J.-P.; Gresh, N. J. Chem. Theory Comput. 2008, 4, 1659. (55) Gordon, M.; Freitag, M.; Bandyopadhyay, P.; Jensen, J.; Kairys, V.; Stevens, W. J. Phys. Chem. A 2001, 105, 293. (56) Yu, H.; van Gunsteren, W. Comput. Phys. Commun. 2005, 172, 69. (57) Vojtechovsky, J.; Chu, K.; Berendzen, J.; Sweet, R.; Schlichting, L. Biophys. J. 1999, 77, 2153. Li, G.; van Wynsberghe, A.; Demerdash, O.; Cui, Q. in Normal Mode Analysis: Theory and Applications to Biological and Chemical Systems; Cui, Q., Bahar, I., Eds.; Chapman and Hall/CRC Taylor and Francis: Boca Raton, FL, 2006. (58) Phillips, G.; Teodoro, M.; Li, T.; Smith, B.; Olson, J. J. Phys. Chem. B 1999, 77, 2153. (59) Autenrieth, F.; Tajkhorshid, E.; Baudry, J.; Luthey-Schulten, Z. J. Comput. Chem. 2004, 25, 1613. (60) Menyhard, D. Chem. Phys. Lett. 2004, 392, 439. (61) Kuriyan, J.; Wilz, S.; Karplus, M.; Petsko, G. J. Mol. Biol. 1986, 192, 133. (62) Kuczera, K.; Kuriyan, J.; Karplus, M. J. Mol. Biol. 1990, 213, 351. (63) Calaminici, P. Chem. Phys. Lett. 2004, 387, 253. (64) Franzen, S. J. Am. Chem. Soc. 2002, 124, 13271. (65) Polack, T.; Ogilvie, J.; Franzen, S.; Vos, M.; Joffre, M.; Martin, J.-L.; Alexandrou, A. Phys. ReV. Lett. 2004, 93, 018102. (66) Joslin, C.; Goddard, J.; Goldman, S. Mol. Phys. 1996, 89, 791. (67) Madan, B.; Sharp, K. J. Phys. Chem. B 1997, 101, 4343.
JP105595Q
