Solvent Polarization

14492

J. Phys. Chem. 1996, 100, 14492-14507

QM/MMpol: A Consistent Model for Solute/Solvent Polarization. Application to the Aqueous Solvation and Spectroscopy of Formaldehyde, Acetaldehyde, and Acetone Mark A. Thompson EnVironmental Molecular Sciences Laboratory, Pacific Northwest National Laboratory, Richland, Washington 99352 ReceiVed: March 6, 1996; In Final Form: May 10, 1996X

We present a hybrid quantum mechanical/molecular mechanical theory that consistently treats many-body polarization effects in the MM region (QM/MMpol). The method described here is based on our previous theory for single-point ground and excited states. Here we extend the formalism to include analytical gradients, making the approach useful for molecular dynamics (MD) simulations. Central to our development is a consistent treatment of the interaction of the MM polarizable dipoles with the full QM wave function in a fashion that allows for energy conservation in the MD. We present implementation details for NDDO semiempirical QM Hamiltonians of the MNDO type. We also present an analysis of the solvent and substituent effects on the spectroscopic blue shift of the n-π* electronic excited state of a series of carbonyl-containing solutes: formaldehyde, acetaldehyde, and acetone. We use the AM1 semiempirical Hamiltonian to treat the solute and the POL1 polarizable model for the solvent. For formaldehyde, we present optimized structures, binding energies, and MD results for CH2O/(H2O)n clusters (n ) 1-6) and describe the effect of specific water/chromophore interactions. The MD results show that CH2O/(H2O)3 exhibits the most favorable solvation with the solute and the three waters forming ringlike structures with near optimal hydrogen bonding interactions. The n ) 3 cluster also has the largest solvent-induced blue shift of the clusters studied, with a value of 730 cm-1 at 150 K. Aqueous bulk-phase MD simulations of formaldehyde, acetaldehyde, and acetone at 300 K show the average number of waters in the first solvation shell of the carbonyl oxygen is 1.2 (formaldehyde), 1.8 (acetaldehyde), and 2.2 (acetone). The induced dipole moments in the solutes are calculated as follows: +0.77 D (formaldehyde), +1.09 D (acetaldehyde), and +1.13 D (acetone). The simulation average solventinduced blue shift in the vertical n-π* electronic excited state of the three solutes is 1148 cm-1 (formaldehyde), 1531 cm-1 (acetaldehyde), and 1607 cm-1 (acetone). The acetone n-π* blue-shift result is in good accord with experiment and provides a benchmark for assessing the accuracy of our prediction for the blue shifts of formaldehyde and acetaldehyde. We suggest the present formalism is the logical starting point for including more accurate ab initio approaches in QM/MM methods.

I. Introduction Solvent effects play a central role in chemistry and biochemistry. Recently, microscopic solvation models that combine quantum mechanics (QM) with molecular mechanics (MM) have begun to be used with good success.1-18 These hybrid QM/ MM methods treat part of the system with QM (solutes, enzyme active sites, substrates, etc.) and the remainder of the system with MM (ie., solvent). The two regions are then coupled by a QM/MM interaction Hamiltonian. An advantage of the QM/ MM approach is that one can begin to address phenomena that are inherently quantum mechanical, such as bond-breaking, electronic excited states, charge transfer, and resonance stabilization. As well, QM/MM methods can provide a description of induced electronic polarization in the solute.6,10 Presently, nearly all uses of QM/MM methods employ simple pairwise-additive MM models to treat the solvent. Recently, MM potentials with many-body polarization effects have begun to gain widespread use.19-27 One way to include these higherorder electronic effects is the use of atom-centered polarizable dipoles (hereafter referred to as MMpol). MMpol approaches have been successfully used to describe ion/water interactions,21,26,28 polarization in nonaqueous solvents,29 and modeling clusters.30 Traditionally, discussions about the quality of molecular simulations center on two main factors, accuracy of the potential X

Abstract published in AdVance ACS Abstracts, August 1, 1996.

S0022-3654(96)00690-9 CCC: $12.00

and sufficient sampling of phase space. With the advent of QM/MM methods, a third issue has become import, the consistency of the potential. QM/MM methods allow for increasing the generality and accuracy of a part of the Hamiltonian describing the QM subsystem yet retain the simplicity of a MM approach for the remainder of the system. This can result in an imbalance in the interactions of the system being studied. This problem can be particularly severe in the case of neutral organic solutes where specific solute/solvent interactions may be important and where the energy of these interactions is comparable to the solvent/solvent interaction energy. Improving the accuracy of only part of the system without first addressing the consistent treatment of electrostatics does little to enhance our understanding of the solution chemistry of QM solutes. A more consistent treatment of many-body effects (polarization) in both the MM and QM regions is the central theme of this and our previous work.13 In an early QM/MM study, Warshel and Levitt included effects from atom-centered polarizable dipoles.1 Later Luzhkov and Warshel (LW) refined this model for both ground and vertical excited states.4,5 While computationally appealing, the LW approach did not consistently treat the interaction of the QM wave function with the MM polarizable dipoles. They convert the QM wave function to atom-centered charges for calculating the interaction of the QM region with the MM region during the dipole SCF calculation. This gives a scheme where the forces are not consistent with the energy. Previously, we © 1996 American Chemical Society

QM/MMpol Model for Solute/Solvent Polarization

J. Phys. Chem., Vol. 100, No. 34, 1996 14493

qm Zaqm el,stat H ˆ QMMM ) -∑ ∑ + ∑ ∑ i m rim a m ram

(2)

The two terms on the right hand side of eq 2 are the interaction of the MM charges with the QM electrons and nuclei, el,pol respectively. H ˆ QMMM contains the interaction between the MM atom-centered polarizable dipoles and the QM electrons and nuclei, respectively. Figure 1. Schematic representation of the fundamental components of a QM/MMpol system: Za, QM nuclei (or effective cores); Ψ, the QM wave function (electrons); qm, MM atoms with charge q; and µm, the MM atom-centered polarizable dipoles. The interactions between the components are represented as numbered arrows which correspond to terms in our electrostatic system Hamiltonian (see text).

presented an alternative QM/MM formalism that treats the QM wave function/MM dipole interaction in a consistent manner.13 We presented the theory for single-point ground and excited states, with fixed nuclear geometries in both the solute and solvent, and used it to study the excited states of the bacteriochlorophyll-b dimer in the photosynthetic reaction center of Rhodopseudomonas Viridis using an INDO/s QM Hamiltonian. In this paper, we extend the QM/MMpol formalism to include analytical energy gradients, which makes the approach useful for geometry optimizations and molecular dynamics simulations. In section II, we review the basic QM/MMpol theory and present the analytical gradient formalism. At this level, the QM/MMpol method is equally applicable to semiempirical, ab initio, and density functional (DFT) levels of QM. In section III, we describe details necessary for implementing this approach for a semiempirical-based QM/MMpol approach. In section IV, we present general calculation details. In section V we present results of a study on the aqueous solvation of formaldehyde, acetaldehyde, and acetone, and in section VI, we present conclusions. II. Theory QM/MMpol Ground-State Hamiltonian and Energy Gradient Formalism. Previously, we described the theory for the ground- and excited-state QM/MMpol Hamiltonian, to which the reader is referred for complete details.13 Here we summarize the salient features of the theory that lead to the final groundstate energy expression which we use to further develop the analytical gradients. Our model system consists of QM and MM components (Figure 1). The QM components consist of the electrons and nuclei (or effective cores). The MM components consist of charged atomic centers with atom-centered polarizable point dipoles.19,23,24 As before, we assume there are no bonds between QM and MM atoms. The treatment of polarizable MM atoms on such bonds requires additional considerations. The interactions, depicted by the arrows in Figure 1, correspond with the electrostatic terms in our system Hamiltonian given below. Thus, we write an effective Hamiltonian for this system (Figure 1) 0 el,stat el,pol vdW el,stat H ˆ sys ) H ˆ QM +H ˆ QMMM +H ˆ QMMM +H ˆ QMMM +H ˆ MM + el,pol bonded vdW +H ˆ MM +H ˆ MM (1) H ˆ MM 0 is the usual gas-phase time-independent, BornWhere H ˆ QM Oppenheimer, nonrelativistic QM Hamiltonian.31 In the equations to follow, the subscripts (i,j), (a,b), and (m,n) refer to QM el,stat is electron, QM nuclei, and MM atoms, respectively. H ˆ QMMM the interaction between the QM system and the MM atomcentered charges (the only electrostatic coupling term in simpler QM/MM theories).

el,pol ) ∑∑∇i H ˆ QMMM i

m

() 1

rim

‚µm - ∑∑∇a a

m

() Za

ram

‚µm

(3)

el,stat is the usual MM charge/charge interaction (Coulomb) H ˆ QMMM term

el,stat ) H ˆ MM

1

∑∑ 2 m n*m

qmqn (4)

rmn

el,pol contains pure MM terms involving the MM polarizable H ˆ MM dipoles.

el,pol )∑ H ˆ MM m

∑ n*m

()

µm‚∇m

qn

rmn

+

1

() 1

∑ ∑ µm‚∇m∇n r 2 m n*m

µn +

mn

µm‚µm

∑ m 2R

(5)

m

The three terms on the right hand side of eq 5 are as follows: (1) the point charge/dipole interaction, (2) the dipole/dipole interaction, and (3) the self-energy for creating the dipoles, respectively.24 The isotropic polarizability of the MM atom is given by Rm. The form of the self-energy assumes linear vdW ˆ QMMM and response in the MM polarizable dipoles.24,32 H vdW H ˆ MM describe the van der Waals interactions between QM/ MM and MM/MM nonbonded atom pairs, respectively, which we treat as a Lennard-Jones 6-12 potential described previbonded ˆ MM contains bonded terms for the MM interously.10,11,13 H actions (e.g., bonds, bond angles, torsions, etc.) and can be implemented in the forms commonly used in many MM force fields: Morse bonds, harmonic bond angles, truncated Fourier expansion, respectively.33 The total system energy is given by el vdW vdW bonded + EQMMM + EMM + EMM Esys ) Esys

(6)

el 0 el,stat el,pol el,stat + 〈Ψsys|H ˆ QM +H ˆ QMMM +H ˆ QMMM |Ψsys〉 + EMM + Esys el,pol (7) EMM el describes the energy from the electrostatic contriWhere Esys el,stat el,pol and EMM are the energies corresponding butions only. EMM to eqs 4 and 5, respectively. The last three terms on the right hand side of eq 6 are nonelectrostatic terms and do not directly affect Ψsys or {µm}. Their energy contribution and energy gradient components are solved in standard fashion and are added after the electronic terms are solved. For systems where we wish to compute energy gradients or employ MC simulations, H ˆ sys (eq 1) and Esys (eq 6) are the appropriate expressions. For simplicity in the development below, we will not explicitly include the nonelectrostatic terms in the equations, though their contributions are certainly included in the final energy and gradient expressions.

14494 J. Phys. Chem., Vol. 100, No. 34, 1996

Thompson

We solve for a ground state, with a fixed nuclear configuration, that is stationary with respect to Ψsys and {µm} by first assuming instantaneous response between Ψsys and {µm} and then employ a coupled-SCF procedure, which we described previously.13 The coupled-SCF procedure consists of a single cycle of MM-SCF (to solve for the {µm}), followed by a single cycle of QM-SCF (to solve for Ψsys). The system-SCF consists of a single cycle of MM-SCF plus QM-SCF. In solving the coupled SCF, we noted the importance of two identities that must be obeyed in the MM- and QM-SCF procedures. el,stat 〈Ψsys|H ˆ QMMM |Ψsys〉 ) ∑qmφm(QM)

(8)

el,pol ˆ QMMM |Ψsys〉 ) -∑µm‚Em0 (QM) 〈Ψsys|H

(9)

m

m

the corresponding ground-state MM dipoles, and Em0,g is the permanent field from the QM for the final converged wave function. At system-SCF convergence, eqs 10 and 11 are equivalent. Energy gradients for the first and second terms on the righthand side of eq 11 are straightforward; gradients for el,stat 0 〈Ψg|H +H ˆ QMMM |Ψg〉 are calculated in the same fashion as ˆ QM standard QM/MM methods,3 but with the appropriate QM/ el,stat MMpol wave function, and gradients for EMM are calculated in the usual fashion for charge/charge Coulomb interactions.33 For the third and fourth terms, recall that the energy gradient of a collection of induced dipoles (obeying linear response) is given by Ahlstrom.24

Gk,ind ) -∑µmg (∇kEm0 (MM) + ∇kEm0,g(QM)) m

φm(QM) is the electrostatic potential from the QM (wave function and nuclei) at MM atomic center m, and Em0 (QM) is the corresponding permanent external electric field from the QM. Equation 8 is the interaction energy of the MM atom-centered point charges with the QM, and eq 9 is the interaction energy of the MM atom-centered dipoles with the QM. The form of eq 9 assumes instantaneous response between Ψsys and {µm} (i.e., Em0 (QM) is the permanent electric field experienced by el,pol contains the frozen MM dipoles). MM dipoles and H ˆ QMMM The left hand sides of eqs 8 and 9 are the forms of these electrostatic interactions used in the QM-SCF procedure, while the right hand sides are used in the MM-SCF. Note that eq 9 is rigorously true only at convergence of the system-SCF. We stress the importance of eqs 8 and 9 as checks on a proper solution of the coupled-QM/MMpol system. Furthermore, these identities are central to the development of the analytical gradient formalism. Utilizing the identities in eqs 8 and 9, we obtained our final working equation for the ground-state energy at system-SCF cycle n.13 el,stat el,n n 0 el,stat el,pol,n-1 n ) 〈Ψsys |H ˆ QM +H ˆ QMMM +H ˆ QMMM |Ψsys 〉 + EMM Esys 1 1 µmn‚Em0 (MM) + ∑µmn‚Em0,n(QM) (10) ∑ 2m 2m el,pol,n-1 where H ˆ QMMM indicates that we used {µmn-1} to solve for the QM wave function at cycle n. Equation 10 is the total energy expression we use to follow system SCF convergence for ground-state calculations at a given nuclear configurations. There are two equivalent approaches to formulating analytical gradients of eq 10. In both approaches, we use the identity in eq 9 to reduce the number of terms which we must differentiate in eq 10. We outline both of these approaches below. We will explicitly discuss the gradients for just the electronic portion of the total system energy (eq 7). The vdW and MM-bonded terms are calculated and added in the usual fashion.10,11,33 Analytical Gradients: Approach 1. After converging the ground-state energy, we can rewrite eq 10

el 0 el,stat el,stat ) 〈Ψg|H ˆ QM +H ˆ QMMM |Ψg〉 + EMM Esys 1 1 µmg ‚Em0 (MM) - ∑µmg ‚Em0,g(QM) (11) ∑ 2m 2m el,pol term with the We have used eq 9 to combine the H ˆ QMMM fourth term on the right-hand side of eq 10. Ψg represents the final ground-state QM wave function, reflecting the presence of both the MM atom-centered charges and dipoles, {µmg } are

∇k(µmg ‚Tmk‚µgk ) ∑ m

(12)

where Gk,ind is the gradient for MM atom k and Tmk is the standard dipole tensor from all other MM dipoles.24,32 Note the presence of the extra term, Em0,g(QM), in the permanent field that represents the external field from the QM wave function and nuclei. Em0,g(QM) is given by

Em0,g(QM) ) -∇mφm(QM)

(13)

where φm(QM) is the molecular electrostatic potential at point m, from the QM electrons and nuclei, and is given by

〈

|r |

φm(QM) ) ∑ Ψg i

-1 im

〉

Za Ψg + ∑ a ram

(14)

Calculating energy derivatives using eq 12 will require us to take derivatives of Em0,g(QM) which will involve calculating second derivatives of the one-electron matrix elements that make up the QM electrostatic potential. These derivative integrals are certainly feasible but may be computationally too demanding and contrary to keeping the approach simple and computationally straightforward. Interestingly, the LW approach circumvents this difficulty by recasting the QM wave function and nuclear charges into atom-centered charges after each QM-SCF cycle.4,5 The MM-SCF, then, sees an Em0,g(QM) arising from purely classical point charges, which makes solving the gradients of the electric field component identical to the purely classical expressions.24,33 The LW scheme is very attractive from this computational perspective. However, the LW approach does not obey the consistency relationship in eq 9 since recasting the QM into atom-centered charges gives a result in the external field experienced by the MM-SCF not being consistent with the external dipole potential experienced by the QM-SCF. To even approximately obey eq 9 would require that the QM atomcentered charges be fit to the molecular electrostatic potential, which would then add greatly to the expense of the calculation. Empirical schemes that would approximate ESP-fitted charges, and were very fast, would be attractive. For the present, we next proceed with an alternative but equivalent approach that provides a more rigorous solution to energy gradients of eq 10 and is computationally efficient. Analytical Gradients: Approach 2. At system-SCF convergence, we can use eq 9 to alternatively rewrite eq 10 as


〈 |

J. Phys. Chem., Vol. 100, No. 34, 1996 14495

| 〉

1 el,pol el 0 el,stat Esys ) Ψg H ˆ QM +H ˆ QMMM + H ˆ QMMM Ψg + 2 1 el,stat - ∑µmg ‚Em0 (MM) (15) EMM 2m where we have now combined the fourth term on the rightel,pol term in the braket hand side of eq 10 with the H ˆ QMMM 1 expression (note the factor of /2). Equations 10, 11, and 15 are completely equivalent at system-SCF convergence. Calculating energy gradients for the second and third terms on the right-hand side of eq 15 involves standard approaches, as discussed above. The only new gradient expression for which we must solve is the bracket expression in eq 15 which includes el,pol . We address this derivatives involving the operator H ˆ QMMM term below. The contribution to the fully converged ground-state energy el,pol term is given by (eq 3) (eq 10) from just the H ˆ QMMM el,pol ˆ QMMM |Ψg〉 ) Tr(PHpol) - ∑∑∇a 〈Ψg|H a

m

() Za

ram

‚µm (16)

where P is the one-particle density matrix corresponding to Ψg, el,pol ˆ QMMM operator Hpol is the one-electron matrix piece of the H evaluated over the basis set of the QM molecule, and Tr signifies the trace of the matrix product. The second term on the righthand side of eq 16 is simply the classical charge/dipole interaction between the nuclei (or effective cores) and the MM dipoles. The matrix elements of Hpol have the form

(

|∑ r |σ ) rim‚µm

λa

m

3

(17)

b

III. Implementing QM/MMpol with NDDO Semiempirical Hamiltonians of the MNDO-Type

im

where λa and σb are atom-centered basis functions on QM atom b. While energy gradients of the matrix element in eq 17 are feasible, we choose to avoid computing the explicit integrals above by approximating the point dipoles as they are included in the QM. We can always approximate a classical point dipole as a collection of six point charges, oriented along the Cartesian axes.

µmx ) lim [q+x(xm + δ) + q-x(xm - δ)] δf0

(18)

where we have shown only the x-component of the dipole (similar expressions apply for the y and z axes). {q+x, q-x} are a pair of charges having the same magnitude but opposite sign, placed along the x-axis at a distance δ equidistant from the origin. Assuming a finite value of δ, we can now reexpress the operator in eq 3 el,pol = ∑∑ H ˆ QMMM i

m

6

qp

∑ p)1, r (p∈m)

imp

- ∑∑ a

a point dipole into six point charges, and maintain high accuracy with the formal theory, than it is by recasting the QM waVe function into atom-centered charges. Our approach still obeys the formal theory, albeit with the approximation of pointcharge dipoles; we retain the full QM wave function and the MM dipoles. Furthermore, our approach allows us to simply replace the expression in eq 19 with the proper analytical matrix elements (eq 17) in the future, with no additional changes required in our formalism or the overall implementation strategy. Last, we show in section III that the approximation in eq 19 obeys eq 9 to a high degree of accuracy. We employ eq 19 in both the ground-state system SCF (eq 10) and the modified energy expression we use for calculating gradients (eq 15). Given the above expressions for the matrix elements of the el,pol H ˆ QMMM operator, calculating the gradients is now similar to el,stat in stanhow we treat the point-charge expression in H ˆ QMMM dard QM/MM approaches. The only additional concern is that el,pol will be given explicitly with energy derivatives of H ˆ QMMM respect to the six point charges in the dipole expansion. However, by implicit differentiation, we wish to have gradients with respect to the MM atom itself, not the six point charges. Since the point charges are distributed at fixed distances and in a fixed orientation about the MM atom, it is straightforward to distribute the force contribution to the MM atom from these six point charges. How we do this for semiempirical Hamiltonians is discussed in section III. The theory, as developed to this point, is equally applicable to single determinant Hartree-Fock (semiempirical and ab initio) and DFT levels of QM theory. We describe implementation details for semiempirical MNDO-type QM Hamiltonians below.

m

6

∑ p)1, (p∈m)

Zaqp rimp

(19)

This approach can always be made sufficiently accurate by adjusting the distance δ and the magnitude of the charges, provided that the induced dipoles are not too large and that a QM and MM atom pair do not get too close together. (An error analysis from using eq 19 is given below.) The approximation we have made in approach 2 is reminiscent of the LW approach. But here we haVe recast the atomcentered dipoles into a set of six point charges as opposed to recasting the QM waVe function into atom-centered charges. The appeal to our approach is that it is much easier to recast

A. Basic Expressions. We now describe details specific to implementation of the QM/MMpol equations using MNDObased semiempirical QM Hamiltonians. (Hereafter, we will refer collectively to MNDO, AM1, and PM3 Hamiltonians as MNDO methods.) We begin by assuming the separate QM and MMpol methods are adequate to describe their respective subsystems and will not be further modified. Referring to eq 0 el,pol bonded vdW vdW ,H ˆ MM ,H ˆ MM ,H ˆ MM , and H ˆ QMMM in the usual 1, we treat H ˆ QM 3,10,33 Here, we limit fashion for both energy and derivatives. discussion to the two QM/MM electrostatic coupling terms el,stat el,pol el,stat and H ˆ QMMM . While the implementation of H ˆ QMMM has H ˆ QMMM 3 been described previously, we find there are certain internal consistency requirements, imposed by the addition of el,pol el,stat , which affect how we treat H ˆ QMMM . For review, recall H ˆ QMMM el,stat from Figure 1 that H ˆ QMMM contains interactions 7,8, the interactions between the MM atom-centered charges with the el,pol consists QM nuclei and QM electrons, respectively. H ˆ QMMM of 9,10, the interactions between the MM atom-center dipoles and the QM nuclei and electrons, respectively. el,stat . The QM electron/MM charge We begin with H ˆ QMMM expression formally involves matrix elements of the form

(

| ∑a r

λa

|σ )

qm

a

(20)

im

For MNDO methods, we treat these integrals in the same fashion as Field et al.3 Thus, we use a two-electron expression that scales the MM charge as an s-orbital charge distribution on the MM atom.

14496 J. Phys. Chem., Vol. 100, No. 34, 1996

(||) λa

qm σ rim a

MNDO

Thompson

( || ) 1 s s rim m m

w qm λaσa

(21)

In the MNDO formalism, the two-center two-electron integrals are computed as the interaction of a finite multipole expansion of charges about the relevant atoms.34 For example (in atomic units)

( | | ) sasa

1 1 s s ) 2 QM 2 0.5 rim m m [Ram + (F0 + FMM 0 ) ]

(22)

where, for simplicity, we have shown the expression for an s-orbital distribution on both the QM and MM atoms, where the point charge expansion is simply an atom-centered charge on each center. The term in the denominator is the OhnoKloppman scaling of the charge/charge interaction which uses MM 34 In general, these empirical parameters, FQM 0 , and F0 . parameters are specific to the orbital combination on the atom. QM QM corresponding Thus, there are parameters, FQM 0 , F1 , and F2 to SS, SP, and PP shell combinations which are treated as monopole, dipole, and quadrupole expansions.34 Specific values are given by Field et al.3 for FMM 0 We also treat the QM core/MM charge term in a similar fashion to that of Field et al.3 However, here we further neglect the Gaussian expansion terms for the QM atoms that are used in AM1 and PM3 levels of theory (which are retained by Field et al.). There is no compelling reason to include the Gaussian terms in the QM/MM interaction as much of their physics is embodied in the QM vdW parameters. Thus, our QM core/ MM charge interaction is given by

( | | )

qmφm(QM,core) ) Zaqm sasa

1 s s [1 + ram m m

core exp(-Rcore a ram) + exp(-Rm ram)] (23)

(to avoid confusion, we use Rmcore (the MNDO exponential scale factor) to distinguish it from the MM atom’s isotropic polarizability, Rm). Here, φm(QM,core) is the electrostatic potential from the MNDO atomic core at MM atom m. el,pol . Recall, we use Next, we describe out treatment of H ˆ QMMM the expression given in eq 19 where we expand the MM dipoles as a finite collection of six point charges. The charge distribution we use for the MM dipoles is exactly analogous to the charge distribution used in the MNDO methods for an SPorbital expansion.34 We exploit this similarity in solving for el,pol . For the the matrix elements and derivatives of H ˆ QMMM el,pol we can electron/dipole interaction matrix elements of H ˆ QMMM write

(

λa

rim‚µm

|r | 3

t)1 t∈m

im

6

qt ∑ t)1 t∈m

) ( ( | | ) ∑( 6

σa = ∑ λaσa

m

λaσa

1

ritm

|r ) qtm

MNDO

w

itm

stmstm )

λaσa

ξ)x,y,z

|r1 |s p

)

m mξ

im

′′ (24)

The first term is the exact expression. The first approximate equality is from eqs 17-19 above. In the third term we see that, in the MNDO formalism, these one-electron integrals are treated with s-orbital distributions similar to the approach used el,stat . In the fourth term, we exploit the MNDO twowith H ˆ QMMM electron integral charge distributions to write a more compact and efficient way to calculate these point charge terms over the

three cartesian components of an SP-orbital distribution on the MM atom (SPx, SPy, and SPz) rather than over six point charges. To make this last term an equality with the third requires that ) FMM in the we always set the empirical parameters FMM 1 0 34 Ohno-Kloppman scaling term. We denote the integral in the fourth term with a (′′) to signify that the different charge magnitude, sign, and distance from the atomic center for the MM dipole charges must be accounted for in the value of the integral (these values are different from the standard values used in the derivation of the original two-electron integral expressions of Dewar and Thiel).34 In the MNDO integral formalism, the distance of the point charges from the atomic center is given by the distance parameter D1 for the SP-orbital combination.34 To account for our MM dipole charges having a different distance from the atomic centers, we merely replace D1 with the distance we have chosen for our MM dipole charge distances. On the basis of numerical tests discussed below, we employ 0.1 Å as the separation of the MM dipole charges from their corresponding atomic center. To account for the different charge magnitude and sign, recall that, in the MNDO approach, only the symmetry unique integrals are calculated in the local diatomic frame.34 These integrals are then transformed back into the global coordinate frame. We exploit this transformation to scale the “native” MNDO integrals to the QM/MMpol electron/dipole integrals in the global frame. For example, for a QM atom with an S-orbital distribution interacting with another atom having an SP-orbital expansion, we have for the standard MNDO integrals.

[ ] [ (sasa|sbpxb) (sasa|sbpyb) (sasa|sbpzb)

(sasa|sbpσb) ) T‚ 0 0 global

]

(25) local

where the 3 × 3 matrix T is the local f global transformation matrix composed of standard Eulerian angles.35 We have shown this example for simplicity. For higher moment expansions on both atoms, there are both left- and right-hand transforms, but this does not change our approach. To account for the specific values of the MM dipole charges, we modify the local f global transformation

[

(sasa|smpxm)′′ (sasa|smpym)′′ (sasa|smpzm)′′

]

[

) C‚T‚ QMMMpol,global

(sasa|smpσm) 0 0

]

(26) local

where we have an added 3 × 3 matrix C in the transformation. The only nonzero elements of C are the diagonal elements Cxx, Cyy, and Czz. The values of these elements are simply the ratio of the MM dipole charge (including sign) with the specific charge used in the standard MNDO method for the native integral (i.e., we perform a simple linear scaling). For an MNDO sp-orbital distribution, charges of +0.5 (+x axis) and -0.5 (-x axis) are employed in the standard expressions.34 Thus, if our MM dipole had an x-component (in the global frame) composed of a -0.2 au charge along the +x axis (and a +0.2 charge on the -x axis), we would have CXX ) -0.2/ 0.5. Note, we take the ratio of corresponding +x axis charges for both MNDO native integral and the MM dipoles in our example. One could obviously use the corresponding negative axis values as well. Thus, the matrix elements and their derivatives are calculated using the standard MNDO approach, where the only modifications are the change in the MM charge distance, D1, and the scaling factors above. el,pol , we use For the QM nuclear core/MM dipole part of H ˆ QMMM the expression


(

J. Phys. Chem., Vol. 100, No. 34, 1996 14497

)

6 Zaqt MNDO ram‚µm m 1 ′′′ Za =∑ w ∑ Za sasa smpmζ (27) 3 r r t)1 ξ)x,y,z ram atm am t∈m

| |

Again, the first approximation is due to eqs 17-19. In the third part of eq 27, we express the interaction with the dipole charges as scaled SP expansions centered on the MM atom, as we did above. Note, our use of (′′′) is meant to signify that we perform the same integral scaling as done in eqs 24-26, but here we further modify the interaction of the nuclear core with each point charge in the MM SP distribution by the term [1 + exp(-Rcore a ratm)], which is the same as that used in eq 23. This latter scaling term is important since, in the MNDO formalism, is a property of the QM the exponential parameter Rcore a nuclear core and thus the exponential term must be included in all expressions that involve the electrostatic potential and electric field from the QM nuclear cores. Note, we do not have a corresponding exp(-Rmcoreratm) term in eq 27 since Rmcore is a property of the MM atomic center, not the MM dipoles. It is very appealing to use the SP-orbital expansion for representing the six dipole point charges on the MM atoms. From a theoretical perspectiVe, we see how the addition of MM dipoles treats electrostatics in a fashion that is completely consistent with the way that QM electrostatics are treated in the MNDO methods. MNDO (and all ZDO-based semiempirical methods) is a Hartree-Fock approach that provides an ansatz for writing the elements of the Fock matrix directly in an orthogonal basis framework.36 However, the two-center matrix elements of the MNDO Fock matrix are very “molecular mechanics-like” in their formulation. One is tempted to view the MNDO methods as an elegant MM method; one-center atomic components are QM, and two-center electrostatic terms are classical. From a computational perspective, our approach to treating el,pol is very efficient since we compute fewer integrals in H ˆ QMMM the local frame, and we also get the energy gradient contribution from the six charges automatically distributed to the MM atomic center since this implicit differentiation is already part of calculating derivative two-electron integrals in our MNDO twoelectron integrals code. Thus, we have implemented the QM/ MMpol method with virtually no new integrals code needing to be written beyond what we already have in our standard QM/ MM method.10-13,37 All we have to do is slightly modify the existing integrals code to account for the addition of the MM dipoles and assemble these terms separately into the corresponding one-electron matrices. Implementing QM/MMpol is quite straightforward, given an existing QM/MM implementael,pol is tion. The additional computational cost of the term H ˆ QMMM negligible compared to the added cost of doing the MM solvent with a MMpol method. We are able to treat the correct el,pol (within a finite dipole approximation) in formalism for H ˆ QMMM an efficient manner. B. Parametrization Issues. The QM/MMpol parameters MM core specific to the MM atoms include FMM 0 , F1 , and Rm . Here, MM we have the additional parameter, F1 , due to our use of the SP-orbital distribution to treat the MM dipoles. After Field et ) 0 and Rmcore ) 5.0 Å-1 in al., we choose to initially set FMM 0 all expressions. As mentioned above, we are required to always ) FMM to obey the last equality in eq 24. The QM/ set FMM 1 0 MMpol parameters specific to the QM atoms are the LJ parameters, and σ (the same as for simple QM/MM methods). Our approach to developing values of the QM/MMpol parameters for specific cases is to employ small clusters that contain the essential interaction in which we are interested and

Figure 2. Error in eq 9 vs the water oxygen-Cs+ distance. The error is defined as the QM form for the electron/dipole interaction (with dipoles expressed as a finite charge distribution) minus the MM form for the electron/dipole interaction (using the exact point-dipole expression) (see text). Shown are the errors for dipole finite charges at 0.1 and 0.2 Å from the MM atomic center.

to fit QM/MMpol structures and binding energies using accurate ab initio results (or an experimental result). We then use these parameters to treat systems both at the cluster level and fully condensed phase. As we show in Results and Discussion, the QM/MMpol Hamiltonian is remarkably robust in obtaining good agreement with accurate ab initio calculations. C. Error Analysis for H ˆ el,pol QMMM. Because we expand the point dipoles from the MM-SCF calculation as six point charges, we will incur some error relative to the exact expression. As δ gets larger (eq 18), the point-charge distribution becomes a poorer approximation to the perfect point-dipole field. As δ gets too small, numerical instabilities would occur if the resulting point charges become too large. In the MM-SCF we use the exact point-dipole operator in all expressions.13,24,32 At the end of each MM-SCF, we recast each atom-centered MM dipole into a set of six point charges at a fixed distance from the MM atomic center. We can systematically test the error in this approximation by calculating both sides of eq 9; the left side using the point-charge approximation and the right side using the exact point-dipole terms. We choose as a “worst case” example, H2O/Cs+, where we treat H2O with AM1 and Cs+ with the MMpol parameters of Smith and Dang.30 Note, Cs+ has a large value for its polarizability (2.571 Å3), which should help to magnify any errors in eq 9. We tested two values of δ ) 0.1 and 0.2 Å (these values are comparable to the distance parameters used in the MNDO two-electron integrals).34 In Figure 2, we show the error in eq 9 as we vary the Owater-Cs+ distance along the C2 axis. At each point, we calculate both sides of eq 9 from a fully converged QM/MMpol system SCF. The error is defined as the QM side minus the MM side of eq 9. Not surprisingly, at long distances the error is very small and nearly identical for both values of δ (a point-charge expansion is a good approximation to the dipole field at large distances).32 For both values of δ, the error is essentially zero at distances typical of nonbonded interactions (>3.5 Å). For distances closer than ∼2.5 Å, the difference in error for the two values of δ becomes larger as expected, with δ ) 0.2 Å having the larger error. However, even at these close distances, the error is always quite

14498 J. Phys. Chem., Vol. 100, No. 34, 1996 small relative to the total energy of the cluster. We also show in Figure 2, the Owater-Cs+ distance for the energy minimized structure. (We have not fit any QM/MMpol parameters for the H2O/Cs+ cluster, nor do we propose treating water with AM1 in general. Rather, we merely use this cluster as a test case for the present error analysis.) At this distance, we see that δ ) 0.1 Å shows a very small error, having a magnitude of ∼0.001 kcal/mol. At distances typical of QM atom/MM atom interactions for neutral organics in polar solvents, we find the error is even smaller, being ∼0.0001 kcal/mol (using δ ) 0.1 Å) (data not shown). Thus, we choose δ ) 0.1 Å in all of our subsequent calculations. For all practical purposes, the relationship in eq 9 is satisfied in our present approach. Further tests of eq 9 involved running microcanonical MD simulations of the H2O/ Cs+ cluster. These tests gave excellent energy conservation with no discernible energy drift (data not shown). IV. Methods All QM/MM and QM/MMpol calculations used the Argus computer program version 3.0.37 The current version of Argus supports QM/MM and QM/MMpol hybrid methods with the following semiempirical Hamiltonians: MNDO,34,38 AM1,39 PM3,40-42 and INDO/S (fixed geometry spectroscopic calculations only).43,44 Both MD simulation and later analysis of the trajectories for spectroscopic averages employed the AM1 Hamiltonian for the solute and the POL1 method for the polarizable water.25 The electronic excited-state calculation consisted of a configuration interaction (CI) of the 12 lowest single-excited configurations (SCI) from the single-determinant ground state of the QM solute. All SCF calculations were of the closed shell restricted Hartree-Fock (RHF) type. Calculated QM-state dipole moments retained all one-center terms but omitted two-center contributions, consistent with the ZDO approximation.36 Transition moments were calculated with the use of the dipole length operator in a similar fashion. In general, the AM1 Hamiltonian underestimates the absolute energy of the n-π* excited state of carbonyl compounds but has recently been shown to reproduce the solvent-induced blue shift in acetone in good agreement with experiment.14 Since we are interested in reporting relative shifts in the excited states (caused by both solvent and substituent effects), the AM1 Hamiltonian will suffice for calculations of n-π* excited state. However, note that while we use an AM1:POL1 approach to generate the trajectory for the ground-state system, we are not constrained to perform the subsequent spectroscopic analysis using the same methods. Indeed we could utilize an ab initio-based QM/ MMpol approach for calculating the excited states. We converge the MM dipoles to 0.001 D in the MM-SCF part of the system-SCF. The total system energy (eq 10) is converged to at least 10-10 au. Energy minimizations were carried out until the maximum component of the gradient for any atom was < 0.0001 au. SHAKE constraints45 were used for rigid POL1 H2O and AM1 CH2O, CH3CHO, and (CH3)2CO C-H bonds. For the QM solute, bond lengths were constrained to the AM1 optimized gas-phase values. Velocities were scaled by the method of Berendsen to simulate contact with an infinite thermal bath.46 Nonbonded interactions employed neutral group cutoffs. The protocol for molecular dynamics simulations of CH2O/ (H2O)n clusters employed a time step of 1 fs. The cluster simulations were performed as follows: (1) the geometry of the system is optimized by using steepest descents for 10 steps; (2) the system is annealed to 10 K using a temperature coupling constant τ ) 0.01 ps; (3) the system is slowly warmed to the target temperature in three equal phases of 2 ps each; (4) the

Thompson data collection phase of the simulation is run for 200 ps (τ ) 0.4 ps) with coordinates and velocities saved to disk every 100 fs for later analysis. The target temperature for the CH2O/ (H2O)n simulations (n ) 1, 2, 3, 4, 5, 6) was 150 K. At higher temperatures, the clusters routinely lost solvent molecules (which was also observed by Fukinaga and Morokuma47). In all cluster simulations, a half-harmonic restraining potential was employed to gently return any waters that had left the cluster. The fully condensed-phase simulations were carried out at constant NVT in a periodic system using the minimum image convention.48 In all cases, the combination of periodic box size and nonbonded cutoffs were such that the QM solute did not interact directly with its image in neighboring boxes. Periodic box size was 20 × 20 × 20 Å. The protocol for preparing these systems and collecting simulation data is similar to our previous QM/MM study of 18-crown-6.10 Both MM/MM and QM/MM nonbonded lists were constructed with 10.0 Å cutoffs, while interactions calculated from these lists were truncated at 9 Å. The time step was 2.0 fs with both NB lists being reconstructed every 15 steps. The total simulation time for the data collection phase was 100 ps. Coordinates and velocities were saved to disk every 100 fs for later analysis. A temperature coupling constant τ ) 0.2 ps was employed. The target temperature was 300 K. To calculate the electronic excitation spectra of the solvated chromophores, the saved trajectories were reprocessed to calculate excited-state energies in the field of the POL1 water molecules. Thus, we performed a QM/MMpol single-point ground- and excited-state calculation on each frame analyzed, where the QM solute was treated with the AM1 method and the solvent was treated with the POL1 method. V. Results and Discussion Our approach to using the QM/MMpol method for clusterand condensed-phase studies is to first fit the free QM/MM parameters to reproduce structures and binding energies of model solute/solvent clusters, obtained at high-level ab initio calculations. If the system Hamiltonian has sufficient flexibility, then the same parameters should be applicable to larger clusters and fully condensed-phase simulations of both the original QM solute, and solutes with related functional groups. Also, by fitting to ab initio values, we can begin to develop a more systematic description of the accuracy of the QM/MM approach. Here, we are interested in describing solvent effects on a simple chromophore, containing a single carbonyl group, from clusters to the fully solvated system. Below we examine the solvent structure and spectroscopy of the series, CH2O (formaldehyde), CH3CHO (acetaldehyde), and (CH3)2CO (acetone). The results we present will be compared to those previously reported by Fukunaga and Morokuma (FM),47 Blair et al.,49,50 Gao,14 and DeBolt and Kollman.51 These four studies reported results on CH2O/H2O clusters (FM using MM), CH2O in the condensed phase (Blair using MM), (CH3)2CO in the fully condensed phase (Gao using QM/MM), and both CH2O and (CH3)2CO in the condensed phase (DeBolt and Kollman using MM), respectively. A. CH2O/(H2O)1: QM/MMpol Parameter Development. We begin by developing a set of QM/MM parameters to describe the structure and binding energies of CH2O interacting with one water molecule. The structures we chose were those recently reported by Dimitrova and Peyerimhoff (DP), shown in Figure 3.52 Structure 1, with the water having one strong H-bond to the carbonyl carbon, is the global minimum. We use the AM1 QM semiempirical Hamiltonian to describe CH2O and the POL1 model25 to describe H2O. As in previous QM/


J. Phys. Chem., Vol. 100, No. 34, 1996 14499

Figure 3. Comparison of QM/MMpol and ab initio binding energies and structural elements for CH2O/(H2O)1 cluster geometries used to develop the QM/MMpol parameters. Ab initio results are from Dimitrova and Peyerimhoff.52

TABLE 1: QM/MMpol Parameters Used in This Study atom Rcore, Å-1 F0, Å F1, Å

σ, Å

, kcal/mol

qb

C O H

Formaldehyde, Acetaldehyde, and Acetone (QM) AM1a AM1 AM1 3.4 0.080 nad AM1 AM1 AM1 3.0 0.190 na AM1 AM1 AM1 2.0 0.010 na

O H

5.0e 5.0

0.0 0.0

Water (MMpol) 0.0 3.205 0.160 0.0 0.0 0.0

R,c Å3

Figure 4. Comparison of QM/MMpol and MP254 (H2O)n, n ) 2, 3, cluster structural elements and the QM/MMpol CH2O/(H2O)n, n ) 2, 3, clusters. (The (H2O)n results are not QM/MM as they are described entirely with POL1.)

na na na

-0.730 0.528 +0.365 0.170

a Standard AM1 values.39 Note, Rcore is a parameter in the corecore repulsion term for AM1, while R is the isotropic atom polarizability we use for the MMpol atoms. b Atom centered charges.25 c Isotropic atomic polarizabilities.25 d Not applicable. e Field et al.3

MM studies, we do not alter either the QM or MM method, but limit our attention to fitting the parameters in the QM/MM coupling terms of the Hamiltonian. The parameters we are free to vary are the Lennard-Jones (LJ) and σ for the QM atomic MM core AM1 parameters for the cores and the FMM 0 , F1 , and Rm 34,39 We choose to use FMM ) 0 Å and Rmcore ) 5 MM atoms. 0 Å-1 after Field et al.3 Also, by the consistency arguments given ) FMM above, we are constrained to set FMM 0 1 . Thus, we varied only the QM LJ parameters in our fit. The structures and binding energies to which we fit are given by DP and are calculated as follows: structures HF/6-311++G(2d,2p) and binding energies MP2/6-31++G(2d,2p) (at the HF geometry) with correction for BSSE.52 All fits of the LJ terms were carried out manually and began with the oxygen, carbon, and hydrogen parameters previously used by us in a QM/MM study of 18crown-6.10 We varied the LJ terms only slightly for the oxygen and carbon atoms of CH2O from our previous values. In Figure 3, we summarize the ab initio results of DP with our QM/ MMpol results. In Table 1, we show the final QM/MMpol parameters employed for this study. The agreement for both structure and binding energies is impressive, given our hand-fit of the parameters and the inherent limitations of the AM1 and POL1 methods. Of particular interest is the relative binding energy for structures 1 and 3. These two interactions may be important to correctly describing the short-range and long-range solvent effect on the n f π* blue shift for CH2O in H2O. The water dimerization energy for POL1 H2O is 5.6 kcal/mol.25 The binding energy of structures 1 and 2 (Figure 3) are comparable to the water dimerization energy, while that for structure 3 is 3 (Vide infra). The largest error is in the H-bonding distance in structure 1, where the QM/MMpol value is ∼0.26 Å shorter than the ab initio result. This is typical for describing hydrogen

Figure 5. Comparison of QM/MMpol and MP254 (H2O)n, n ) 4, 5, cluster structural elements and the QM/MMpol CH2O/(H2O)n, n ) 4, 5, clusters. (The (H2O)n results are not QM/MM as they are described entirely with POL1.)

bonding with QM/MM methods employing semiempirical QM Hamiltonians.7 B. Formaldehyde/(H2O)n (n ) 2-6): Energy Minimized Structures. Using the above QM/MM parameters, we examine the structure and energetics of larger water clusters of CH2O. We obtained energy minimized structures for (H2O)n and CH2O/ (H2O)n clusters (n ) 2-6) by an MD simulated annealing calculation. We first equilibrate the cluster to 150 K and then continuously lower the target temperature to 0.1 K while at the same time continuously scaling the temperature coupling constant from 0.2 to 0.001 ps. The velocities were rethermalized to the target temperature at 100 equally spaced time intervals during the annealing. The total simulation time for the annealing of each cluster was 50 ps. We employ an annealing approach in an attempt to locate the global minimum for these clusters. We checked our protocol by repeating the CH2O/(H2O)3 cluster annealing starting at 300 K; the final result was the same as the 150 K starting temperature. The structural results are summarized in Figures 4-6. We compare the (H2O)n structures with ab initio MP2 and DFT results of Xantheas53,54 and ab initio MP2 results of Tsai and Jordan.55 Note, in Figures 4-6, the (H2O)n results are not QM/MM since the (H2O)n clusters are described purely with MMpol. We present the (H2O)n data alongside the corresponding QM/MMpol CH2O/(H2O)n clusters in order to compare changes caused by adding the QM solute. For the corresponding CH2O/(H2O)n clusters, we do not have

14500 J. Phys. Chem., Vol. 100, No. 34, 1996

Figure 6. Comparison of QM/MMpol and MP255 (H2O)6 cluster structural elements and the QM/MMpol CH2O/(H2O)6 cluster. Also shown is a stereoview of the CH2O/(H2O)6 cluster. (The (H2O)n results are not QM/MM as they are described entirely with POL1.)

ab initio or experimental results with which to compare. It would be interesting to see comparisons of our CH2O/(H2O)n cluster results with ab initio calculations. CH2O/(H2O)2. Figure 4 shows the results for the n ) 2, 3 clusters. The Ow-Ow (Ow ) O water) distance is calculated as 2.800 Å (POL1) and 2.911 Å (MP2). The angle between Ow-Hw of the donor water and Hw‚‚‚Ow of the acceptor water is 176.0° (POL1) and 174.3° (MP2). We calculate the water dimer binding energy as -5.6 kcal/mol (POL1) relative to -5.3 kcal/mol (MP2). The CH2O/(H2O)2 cluster results shows an interesting ring structure where each molecule in the cluster acts as both an H-bond donor and acceptor. The acceptor water (from the water dimer) has rotated along the Ow-Hw‚‚‚Ow axis (with the donor water) and one of its hydrogens now H-bonds with Os (the carbonyl oxygen; “s” ) solute). Our CH2O/(H2O)2 cluster structure is similar to that described by the MD results of FM.47 CH2O/(H2O)3. The (H2O)3 cluster, also shown in Figure 4, exhibits a ring structure in which each water accepts and donates a single H-bond. The POL1 structure is planer, while the MP2 results of Xantheas and Dunning shows each water being slightly rotated with one of its hydrogens out of the plane defined by the three oxygens.53,54 Possible reasons for this difference between POL1 and MP2 could include the constrained internal geometry for POL1 and a more complete treatment of manybody effects in the MP2 calculations. The three Ow-Ow distances are calculated as 2.779, 2.778, and 2.780 Å by POL1 which compare to the MP2 values of 2.798, 2.799, and 2.800 Å.54 The agreement in these values is quite good. The internal Ow-Hw‚‚‚Ow angles are 160, 161, and 161° compared with MP2 values of 148, 151, and 151°. In the CH2O/(H2O)3 cluster we also see a ring structure in which CH2O and each of the waters acts both as a H-bond donor and acceptor. The major structural change is that the (H2O)3 cluster has opened up, allowing CH2O to insert into the water cluster. Note, how the Ow-Hw‚‚‚Ow angles are now closer to linear H-bonds, relative to (H2O)3. The exception is the Cs-Hs‚‚‚Ow H-bond between CH2O and one of the water oxygens, which has a value of 144°. This is expected given the much weaker strength of the H-bond at this end of CH2O (Figure 3). Note the Os‚‚‚Hw-Ow angle is calculated to be 131°, which is near the optimal angle of

Thompson 127° calculated for the QM/MMpol global minimum CH2O/ (H2O)1 cluster shown in Figure 1. The cluster MD results of FM were also consistent with the CH2O/(H2O)3, forming a ringlike structure.47 CH2O/(H2O)4. In Figure 5 we show results for the n ) 4, 5 clusters. The n ) 4 water cluster forms a square structure that is in good accord with the MP2 results.54 The major difference is the MP2 result shows one hydrogen from each water alternately pointing up and down, while the POL1 result gives a planer cluster. The Ow-Ow distance for POL1 is calculated as 2.707 Å relative to an MP2 value of 2.743 Å.54 The CH2O/ (H2O)4 cluster shows CH2O binding to the “outside” of the cluster, relative to the square (H2O)4 ring. The Ow-Ow distance between the two waters that H-bond to CH2O has opened up ∼0.04 Å relative to the (H2O)4 cluster. CH2O binding to the outside of the cluster is in contrast to the n ) 3 result above where we saw CH2O insert into the water cluster. The reason for this is likely the more favorable energetics for inserting CH2O into the (H2O)3 cluster, which makes the OwHw‚‚‚Ow H-bonds significantly more linear as well as giving a Hw-Ow angle much closer to the value for the optimized CH2O/(H2O)1 cluster (Figure 3). Since the Ow-Hw‚‚‚Ow angles are nearly optimal in the (H2O)4 cluster, opening the ring to insert CH2O does not incur as much of an energetic advantage. CH2O/(H2O)5. The (H2O)5 cluster forms a pentagonal structure (Figure 5) qualitatively similar to the MP2 structure obtained by Xantheas.54 Again, POL1 gives a planer structure, while the MP2 structure shows non-H-bonding hydrogens pointing slightly up and down from the plane of the oxygens. The unique Ow-Ow distances are calculated to be 2.687 and 2.689 Å for POL1, while the MP2 values are 2.726, 2.729, and 2.731 Å.54 As with the n ) 3, 4 clusters, POL1 gives the OwOw distances systematically shorter by ∼0.02-0.04 Å relative to MP2. As in the n ) 4 cluster, we see that CH2O/(H2O)5 shows CH2O binding to the outside of the water cluster (Figure 5), where, here, CH2O binds to the outside and slightly above the plane of the water cluster. CH2O/(H2O)6. Last, in Figure 6, we show the results for n ) 6. The (H2O)6 cluster is hexagonal, in qualitative agreement with the MP2/aug-cc-pVDZ optimized structure of Tsai and Jordan.55 In the ab initio cluster, each water has a hydrogen that points alternately up and down around the ring. The POL1 structure shows two hydrogens pointing down, two pointing up, and two in the plane of the water oxygens. The two unique Ow-Ow distances are calculated as 2.683 and 2.684 Å for POL1, while the MP2 result gives a single distance of 2.725 Å. We see an interesting structural change in the CH2O/(H2O)6 cluster. Here, CH2O binds to one face of the plane of the water cluster. Two waters bind directly with CH2O and have been “pulled” in, away from the hexagonal geometry. The resulting structure of the six water oxygens is now two fused water tetramers, each very similar to the structure reported above for (H2O)4. The water that H-bonds with the carbonyl oxygen has rotated one of its hydrogens up to interact with CH2O. By forming two tetramers, the waters have formed an extra H-bond between the bridging waters. There are now seven H2O-H2O H-bonds in addition to the CH2O/H2O interactions, as opposed to six H-bonds in the (H2O)6 cluster. All H2O-H2O H-bonds are nearly linear in both the (H2O)6 and CH2O/(H2O)6 clusters. We next compare structural trends in the clusters. In Figure 7 we plot the average Ow-Ow distances in the (H2O)n clusters to compare the POL1 and ab initio results described above. In all cases the POL1 values are somewhat smaller than the ab initio results. The worst agreement is for the water dimer, which


J. Phys. Chem., Vol. 100, No. 34, 1996 14501 TABLE 2: Energetics for CH2O Binding with (H2O)n Clusters na

∆Edist,b (kcal/mol)

∆Estab,c (kcal/mol)

∆Ebind,d (kcal/mol)

2 3 4 5 6

+0.4 +3.2 +0.5 +0.4 +1.0

-8.3 -10.0 -5.4 -5.6 -7.4

-7.9 -6.8 -4.9 -5.2 -6.4

a Number of waters in cluster. b The energy penalty to distort the optimized water cluster to the structure it has in the CH2O/(H2O)n cluster, but with CH2O still at infinite separation. c The stabilization energy for CH2O binding with the distorted (H2O)n cluster. d ∆Ebind ) ∆Edist + ∆Estab (see text).

Figure 7. Interatomic distances plotted for (H2O)n and CH2O/(H2O)n, n ) 2-6, clusters. Distances shown are as follows: average water oxygen-oxygen (Ow-Ow) from (H2O)n clusters, average Ow-Ow from CH2O/(H2O)n clusters, and CH2O carbon (Os)-Ow distances in CH2O/(H2O)n clusters. QM/MMpol luster geometries are obtained by simulated annealing. Ab initio values are from Xantheas54 and Tsai.55

is expected. The POL1 gas-phase monomer dipole moment is 2.02 D, which is slightly larger than the accurate ab initio value of 1.88 D54 or experimental value of 1.854 D.56 The slightly larger POL1 dipole contributes to the shorter Ow-Ow distances for all the water clusters, and most severely for the water dimer. However, the trends in Ow-Ow distances are in qualitative agreement between the POL1 and ab initio results. In all cases, the Ow-Ow distance gets smaller as the cluster size goes from n ) 2 to n ) 6. Also shown in Figure 7 are the Ow-Ow distances from the CH2O/(H2O)n clusters. Here we plot the average Ow-Ow distance between the waters that still form H2O-H2O H-bonds (neglecting the two waters that interact directly with CH2O). Also shown is the Os-Ow distance between the carbonyl oxygen and the oxygen of the water to which it H-bonds. We observe that Ow-Ow distances are consistently smaller in the CH2O/(H2O)n clusters relative to the (H2O)n clusters. This is indicative of the favorable interaction between the solute (CH2O) and the water cluster, allowing for an enhanced polarization of both the waters and the solute, making the cluster geometry more compact. The Os-Ow distance shows an interesting trend. This distance gets smaller from n ) 2 to its smallest value for the n ) 3 cluster. Then, it increases ∼0.1 Å for the n ) 4 cluster. For n ) 5 and 6, the Os-Ow distance again gets progressively smaller. This latter trend may be indicative of increasing bulk behavior in the OsOw distance. Indeed the first peak in the Os-Ow radial distribution function for aqueous CH2O simulation occurs at 2.7 Å (Vide infra). It is apparent that the n ) 3 cluster enjoys a special stability with respect to the Os-Ow distance. In Table 2, we examine the energetics of the above cluster calculations. We tabulate the energy differences for the following transformations ∆Edist

∆Estab

CH2O + (H2O)n 98 CH2O + (H2O)n* 98 CH2O/(H2O)n* Here ∆Edist is the energy for distorting the optimized isolated water cluster to the structure it will have in the final CH2O/ (H2O)n cluster (but with CH2O still at infinite separation), and ∆Estab is the subsequent stabilization caused by the interaction of CH2O with the distorted water cluster (there is a very small distortion cost for CH2O in ∆Edist of only ∼0.02 kcal/mol for

all the clusters). The total binding energy of CH2O with the cluster is given as ∆Ebind ) ∆Edist + ∆Estab. The trends in the energies (Table 2) reflect the structural results above (Figures 4-6). Thus, we see that the energy penalty to distort the water cluster (∆Edist) is relatively small for the n ) 2, 4, 5, 6 clusters, being between 0.4 and 1.0 kcal/mol. For the n ) 3 cluster, there is a substantially higher energy penalty, where ∆Edist ) 3.2 kcal/mol. The larger ∆Edist reflects the cost of breaking one of the water-water H-bonds (but partially strengthening the other two H-bonds by making them linear). The values of ∆Estab reflect the trends in the Os-Ow distances shown in Figure 7. Thus, we see that n ) 3 has the largest stabilization energy, being -10.0 kcal/mol. The n ) 2 case is next largest with a value of -8.3 kcal/mol. The n ) 4 stabilization energy is only about half that of the n ) 3 cluster, being -5.4 kcal/mol. The n ) 4, 5, 6 stabilization energies show a trend of increasingly negative values from -5.4 to -7.4 kcal/mol, reflecting with the shortening Os-Ow distances in Figure 7. Interestingly, the total CH2O binding energy for the clusters shows that the n ) 2 cluster is slightly more stable than the n ) 3 case (∆Ebind ) -7.9 for n ) 2 vs ∆Ebind ) -6.8 for n ) 3). This is due to the much larger value of ∆Edist for the n ) 3 cluster. The overall agreement between the POL1 and ab initio (H2O)n structures is good. For the corresponding CH2O/(H2O)n clusters we observe, for n > 3, that CH2O preferentially binds to the outside of the water cluster. For n ) 2 there is no “inside” or “outside” of the cluster. For n ) 3, there appear to be unique structural advantages to opening up the water cluster and allowing CH2O to insert between two waters. That the larger clusters show CH2O binding to the outside is reasonable given the relative weakness of H-bonding to the CH2 end of CH2O compared to the water-water interaction. It has been proposed that hydrophobic solvation in liquid water is dominated by pentagonal and hexagonal ring structures.57 It is intriguing to speculate on the possible significance of ringlike structures being important to the solvation of CH2O liquid water. The data in Figure 7 and the energetic results in Table 2 suggest that CH2O/ (H2O)n (n ) 2, 3) ring structures may be present, though at 300 K these clusters may not dominate due to thermal fluctuations. C. Formaldehyde/(H2O)n Molecular Dynamics: Solvent Structure and Spectroscopy. We next performed a series of MD simulations to examine the solvent structure of CH2O/ (H2O)n clusters at a higher temperature. All cluster simulations were run at 150 K for 200 ps (see Methods). These results allow us to compare to the cluster simulations of FM47 and to examine specific solvation effects at the cluster level that may be important to the condensed-phase results. For brevity, we limit discussion of specific solvent structure to the n ) 2-4 clusters. In Figures 8 and 9 we show the distribution of the water around formaldehyde for the CH2O/(H2O)2 cluster simulation.

14502 J. Phys. Chem., Vol. 100, No. 34, 1996

Figure 8. (a) Os-Hw, (b) Hs-Ow, and (c) Os-Ow simulation average distributions from CH2O/(H2O)2 cluster MD simulation at 150 K.

Figure 9. (a) cos(Os-Hw-Ow) and (b) cos(Cs-Os-Ow) simulation average distributions from CH2O/(H2O)2 cluster MD simulation at 150 K.

We show distributions for the oxygen of formaldehyde (Os) (s ) solute) with hydrogen (Hw) and oxygen (Ow) of the water molecules, as well as the hydrogen of formaldehyde (Hs) with Ow. The Os-Hw distribution (Figure 8a) shows a well-defined first peak centered at 1.9 Å that integrates to one Os-Hw interaction. Next is a prominent peak, centered near 3.5 Å, with shoulders at 4.5 and 5.6 Å. The latter three peaks overlap strongly, but each essentially integrates to ∼1 Os-Hw interaction. The 1.9 and 3.5 Å peaks are consistent with one of the waters maintaining a single good hydrogen bond with the carbonyl oxygen of formaldehyde (similar to Figure 3 structure 1). Figure 8b, the Hs-Ow distribution, has a weak shoulder at 2.5 Å and a prominent peak centered at 4.5 Å. There is also a broad shoulder between 5.5 and 8.0 Å. The shoulder at 2.5 Å is indicative of the optimized Hs-Ow interaction (Figure 3 structure 3). The total integrated coordination number for this shoulder is ∼0.25-0.75 Hs-Ow interactions. In Figure8c, the Os-Ow distribution shows a prominent peak centered at 2.9 Å that is also consistent with one of the waters maintaining a single good hydrogen bond to Os. The other water oxygen is broadly distributed at distances up to 6.5 Å away from Os, though there is a shoulder at ∼3.5 Å. In Figure 9 we examine the orientation of the H-bond for the water molecule that interacts with Os. Parts a and b of Figure 9 show the simulation average distributions of the cosine of the Os-Hw-Ow and Cs-Os-Ow angles, respectively. The cos(Os-Hw-Ow) distribution clearly shows that one Hw has a near-linear H-bond with Os (peak at -1), while the other Hw maintains an angle of ∼50°. This is again consistent with the

Thompson


H-bonding interaction shown in the global minimum structure in Figure 3 (structure 1). Figure 9b shows a much broader distribution for the cos(Cs-Os-Ow) function, which suggests the H-bonding water moves reasonably freely about a region defining a cap off the carbonyl oxygen. The distribution is noisy, but there may be a slight preference for angles of ∼120140°, which compares to the optimized angle of 127° for structure 1 in Figure 3. Our result in Figure 9 is similar to that reported by Blair et al. for the same distribution applied to just the first solvation waters around formaldehyde in their aqueous bulk-phase simulations at 300 K.50 Previously, FM reported an Hs-Ow distribution having a more prominent peak at 2.5 Å in their MM study of the CH2O/ (H2O)2 cluster.47 They concluded that formaldehyde and the two waters form a ringlike structure. Visualization of our simulation trajectory shows two general motifs contribute to the structural results for CH2O/(H2O)2: the cyclic structure (Figure 4) and an extended linear configuration, where CH2O H-bonds to one end of the water dimer (not shown). We did not observe any configurations where both waters simultaneously contributed a single H-bond to Os. There may be an increased tendency for our method to have the linear configurations since we do slightly underestimate the binding energy of structure 3 in Figure 3, where a H2O H-bonds to Hs. In Figure 10 we show the Os-Hw, Hs-Ow, and Os-Ow distributions for the CH2O/(H2O)3 simulation. In Figure 10a, the peak at 1.9 Å is more prominent than the same peak from the CH2O/(H2O)2 cluster (Figure 8a), and the minimum between the first and second peaks is deeper. This is consistent with the water that H-bonds directly to Os, being more strongly bound in this cluster. The Hs-Ow distribution shows the 2.5 Å peak is more prominent than the same feature in CH2O/(H2O)2. In CH2O/(H2O)3, the 2.5 Å peak integrates to one Hs-Ow interaction. In Figure 10c, the Os-Ow distribution shows a prominent peak at 2.9 Å which integrates to one water molecule. The second broader peak centered near 4.0 Å shows the remaining two oxygens maintaining similar distances from the carbonyl oxygen. This latter peak is very different from the same region in CH2O/(H2O)2 (Figure 8c). Visualization of the simulation results shows the CH2O/(H2O)3 cluster maintains compact ringlike configurations very similar to the energy



minimized structure shown in Figure 4. Our result for the n ) 3 cluster simulation is similar to that reported by FM.47 In Figure 11, we show the Os-Hw, Hs-Ow, and Os-Ow distributions from the CH2O/(H2O)4 cluster simulation. The distributions are all consistent with the formaldehyde now lying essentially outside the H2O cluster, as we saw in the energy minimized structure (Figure 5). Thus, in Figure 11b, we see there is essentially no structure at 2.5 Å in the Hs-Ow distribution, indicating that the hydrogen end of formaldehyde no longer has significant interactions with the solvent at this temperature. The first peaks in Figure 11a,c indicate that one water maintains a good H-bond to the carbonyl oxygen. However, the remainder of the waters are distributed away from the formaldehyde molecule. Visualization of the simulation results shows configurations similar to the energy minimized structure (Figure 5) predominate. While FM did not study the CH2O/(H2O)4 cluster, they commented that formaldehyde should not be able to effectively compete with the water/water interactions in clusters larger than three waters and would occupy surface states. This statement is consistent with both our simulation results and the 0 K optimized cluster results above. The pattern of formaldehyde solvation with cluster size is further reflected in the solvent-induced blue shift of the n-π* excited state presented below. We next examine the solvent-induced shift in the n-π* vertical excited state of formaldehyde from the cluster simulations. For each configuration saved in the MD simulation, we calculated the ground and vertical excited state of formaldehyde in the field of the cluster of polarizable water molecules. We model the excited state by a configuration interaction (CI) calculation consisting of all single excitations (SCI) from the three highest occupied molecular orbitals (MO) into the four lowest unoccupied virtual MOs of the solute. We used the AM1 method for the spectroscopic calculation and included the polarizable solvent for the ground state in the same manner as that used in the MD simulations. The contributions of the polarizable solvent to the excited states was included in the same fashion as that described in our previous QM/MMpol study.13 Thus, for the excited state, we calculate the QM-wave function/ MM-dipole relaxation only to first order. However, for the n-π* transition, the ground-excited-state dipole moment change is only ∼1.5 D, and thus, we do not expect there to be large changes in the excited-state energies if we were to do a full

J. Phys. Chem., Vol. 100, No. 34, 1996 14503

Figure 12. Calculated blue shift in the vertical n-π* excited state of formaldehyde in CH2O/(H2O)n, n ) 1-6, clusters. Shown are the values for QM/MMpol energy minimized clusters (0 K) and simulation average values from cluster simulations at 150 K. Blue shift is defined in eqs 28 and 29.

QM/MMpol SCF on the excited state. We emphasize, again, that, for ground-state MD energy and gradient calculations, we include the QM-wave function/MM-dipole relaxation to infinite order. We define the solvent-induced shift for a vertical electronic excited state as

〈∆ν〉 )

1

N

∑∆νi N i)1

CI CI ) - E*(Ψgas ) ∆νi ) E*(Ψsol

(28) (29)

where 〈∆ν〉 is the simulation average solvent shift, N is the total number of configurations analyzed, and ∆νi is the energy difference for configuration i, consisting of the energy difference CI )) and between the excited state in the presence (E*(Ψsol CI absence (E*(Ψgas)) of the solvent. In Figure 12 we show the solvent-induced blue shift for the n-π* excited state of formaldehyde for the 150 K cluster simulations, for cluster sizes n ) 1-6. Also shown is the blue shift for the 0 K minimized structures (Figure 3, structure 1, and Figures 4-6). The 150 K results display an interesting trend in the solvent effect on the n-π* excited-state energies. There is a near-linear increase in the calculated blue shift as the cluster increases from one to three waters, where it peaks at 730 cm-1. The blue shift then decreases to ∼470 cm-1 for the n ) 4-6 cluster sizes. The stability of the n ) 3 cluster appears to be reflected in the larger blue shift. The 0 K blueshift values are all greater than the 150 K values. This negative temperature effect is expected as the clusters will explore less stable solute/solvent configurations at higher temperatures. Interestingly, the n ) 3 cluster has a small temperature effect, relative to the n ) 1, 2 clusters, again emphasizing the special stability of the n ) 3 cluster afforded by the cooperative H-bonding with near-optimal H-bonding angles. The simulation results showed the n ) 3 cluster maintained its ringlike structure at 150 K. The larger temperature effect for the n ) 2 cluster is certainly caused by this cluster visiting linear configurations at 150 K. The greater stability of the n ) 2 cluster at 0 K is reduced at higher temperature, relative to n ) 3, possibly by the nonoptimial H-bonding angles in its cyclic structure. The n ) 4, 5 clusters also have an apparent small temperature effect,

14504 J. Phys. Chem., Vol. 100, No. 34, 1996 suggesting they also may have some special stability. This is not the case, however, since formaldehyde occupied a surface configuration in the n ) 4, 5 clusters at 0 K and already showed a relatively weak blue shift. For the n ) 6 cluster, we again see a large temperature effect of ∼300 cm-1. We did not analyze the n ) 6 MD results in great detail, due to the complexity of the cluster geometries visited at 150 K. However, visualization of the simulation results showed CH2O spends most of its time outside the cluster with just the carbonyl oxygen interacting with the cluster. The stability of the fused tetrameric structure, observed at 0 K, is apparently lost at higher temperatures. Simulations for clusters of up to 14 waters gave qualitatively similar results to those shown here; CH2O remains largely outside the water cluster, and its n-π* blue shift is similar to that seen for the n ) 4-6 values (data not shown). Our structural results are in overall agreement with the cluster results of FM. However, FM predict blue shifts considerably larger than our values. They calculate a blue shift of 1696 cm-1 for the n ) 3 cluster at 150 K.47 For the fully condensed phase, FM calculate a blue shift of 2557 cm-1 at 323 K, considerably higher than the experimental value for even acetone of 1560 cm-1.58 We calculate the n-π* blue shift for the gas-phase optimized CH2O/(H2O)1 structure (Figure 3, structure 1) as 767 cm-1. This value compares well with the same value reported by Blair et al. (with an ab initio CISD level of theory) of 600 cm-1 for a similar structure.49 Dimitrova and Peyerimhoff calculate a blue shift of 1140 cm-1 for this same CH2O/(H2O)1 structure using an ab initio CISD calculation.52 The reasons that the DP value is somewhat higher than our and Blair’s values are not clear to us. DP also calculate a very large blue shift of 3800 cm-1 for a CH2O/(H2O)3 optimized cluster. However, their cluster has two waters H-bonding directly to the carbonyl oxygen of formaldehyde and the third water optimally Hbonding to the H2C end of formaldehyde. There do not appear to be any direct H2O/H2O interactions in the DP CH2O/(H2O)3 cluster, in contrast to our results. In our annealed structures and MD simulations, we did not observe multiple waters simultaneously H-bonding to the carbonyl oxygen of CH2O. Furthermore, we see significant water/water interactions in all of the clusters we studied. The DP CH2O/(H2O)3 cluster is perhaps a local higher-energy minimum. There are no experimental results for CH2O/(H2O)n clusters of which we are aware. Yet, based on our overall agreement with FM’s cluster simulation structural results and Blair’s ab initio calculated blue shifts, the present approach likely describes the major structural and spectroscopic characteristics of these clusters correctly. Comparison to experiments would be very interesting. D. Formaldehyde, Acetaldehyde, and Acetone: Molecular Dynamics Simulations of the Aqueous Bulk Phase. We present results of MD simulations of the aqueous condensed phase for formaldehyde, acetaldehyde, and acetone. Recall, we have fit the QM/MMpol parameters to the ab initio CH2O/ (H2O)1 cluster results (Figure 3). We assert that these parameters should be transferable to the series R1,R2CO (where R1 and R2 can be H or CH3), where we can study the substituent effects on the carbonyl group. Furthermore, while there are no experimental values for the n-π* blue shift for monomeric formaldehyde or acetaldehyde in aqueous solutions, acetone does have a well-described n-π* blue shift of 1560 cm-1.58 By comparing our results to the acetone value, we may get some estimate of the accuracy of our formaldehyde result. Also, we should observe a trend in the calculated results for the series formaldehyde, acetaldehyde, and acetone, since successive methyl substitution on the carbonyl group should increase the

Thompson

Figure 13. (a) Os-Hw and (b) Os-Ow radial distribution functions for QM/MMpol aqueous condensed-phase CH2O MD simulations at 300 K.

ground-state dipole moment, and consequently the chromophore’s solvent interaction. Our calculated (and experimental) gas-phase AM1 groundstate dipole moments for the three solutes (at the AM1 optimized geometry) are as follows: formaldehyde, 2.33 D (2.33 D59); acetaldehyde, 2.61 D (2.69 D59); and acetone, 2.87 D (2.88 D59). The AM1 values are in excellent agreement with experiment. For each of the three solutes, we performed a 100 ps. MD simulation in a 20 × 20 × 20 Å box of water at 300 K and constant NVT (Methods). The saved configurations were analyzed to give the simulation average results reported below. First we examine the solvation structure and then the effects of the solvent on n-π* excited-state energies. SolVent Structure. In Figure 13, we show the radial distribution functions, g(R), for Os-Hw (Figure 13a) and Os-Ow (Figure 13b) for formaldehyde. The Os-Hw distribution shows a well-defined peak centered at 1.8 Å. This peak integrates to 1.4 Os-Hw interactions. The second peak in the Os-Hw distribution is much broader and centered near 4.2 Å. This peak integrates to ∼25 Os-Hw interactions. For the Os-Ow distribution in Figure 13b, there is a small first peak centered at ∼2.7 Å and a broader second peak, centered at about 4.0 Å. The 2.7 Å peak integrates to 1.4 waters. The second peak integrates to ∼25 waters. The g(R) distribution for Hs-Ow was unremarkable, indicating no appreciable solvent ordering at the hydrogen end of formaldehyde (data not shown). To better discern the specific solute/solvent interactions in the first solvation shell of formaldehyde, we recalculated the Os-Hw and Os-Ow pair distributions, but included only waters with oxygens within 3.25 Å from Os. These results are shown in Figure 14. In Figure14a, we see a typical Os-Hw profile for a single water molecule maintaining a single H-bond with the carbonyl oxygen. The total curve integrates to 2.3 OsHw interactions (essentially one water molecule). Figure 14b shows the Os-Ow distribution. We see a single well-defined peak, centered at 2.9 Å that integrates to 1.2 water molecules. In Figure 14c, we show the simulation average distribution for the cosine of the angle composed of Cs-Os-Ow for just the first solvation shell. We observe a broad distribution that is centered between 125 and 135°. This is similar to the value of 127° for the CH2O/H2O optimized structure 1 in Figure 3. The Os first solvation shell water molecule maintains a single good H-bond to the carbonyl oxygen, with a structure similar to that observed in the CH2O/(H2O)2 cluster (Figure 8). The result of Figure14c is similar to that reported by Blair.50 Otherwise, our results contrast somewhat with those reported by both FM47 and Blair.50 FM report 6.2 waters in the first solvation shell of the carbonyl oxygen (however, they integrate their Os-Ow distribution out to 4.6 Å). Blair reports 2.6 waters in the first solvation shell of Os (integrated out to 3.5 Å). A more direct comparison between the FM and Blair results is not possible here. Nonetheless FM conclude that their model


J. Phys. Chem., Vol. 100, No. 34, 1996 14505 TABLE 3: Simulation Results for Formaldehyde, Acetaldehyde, and Acetone

solute

no. of H2O in first solvation shella

〈µg〉,b D

induced 〈µg〉,c D

n-π* blue shift,d cm-1

formaldehyde acetaldehyde acetone

1.2 1.8 2.2

3.1 3.7 4.0

+0.77 +1.09 +1.13

1148 1530 1607

a Number of water molecules in first solvation shell of the carbonyl oxygen. b Simulation average AM1 ground-state dipole moment in the presence of the solvent. c Induced ground-state dipole moment, relative to the AM1 gas-phase optimized geometry in the absence of solvent. d The calculated blue shift in the n-π* vertical excited state relative to the gas-phase solute.

Figure 14. (a) Os-Hw, (b) Os-Ow, and (c) cos(Cs-Os-Ow) simulation average distribution functions for QM/MMpol aqueous condensed-phase CH2O MD simulations at 300 K. The distributions examine only solvent molecules whose oxygens are within 3.25 Å of the oxygen of CH2O.

Figure 15. (a) Os-Hw and (b) Os-Ow radial distribution functions for QM/MMpol aqueous condensed-phase CH3CHO MD simulations at 300 K.

Figure 16. (a) Os-Hw and (b) Os-Ow radial distribution functions for QM/MMpol aqueous condensed-phase (CH3)2CO MD simulations at 300 K.

shows formaldehyde to be more strongly solvated than the study of Blair. Blair shows a well-defined strong first peak in the Os-Ow distribution50 in contrast to the much smaller first peak we observe in Figure 13b. They do show an Os-Ow second solvation shell peak containing 19 waters similar to our value of 25. Overall, our results show formaldehyde to be less strongly solvated than the Blair results by a factor of ∼2 and apparently considerably less strongly solvated than the FM results. This difference is reflected in the calculated blue shifts (Vide infra). In Figures 15 and 16 we show the structural analysis of the acetaldehyde and acetone simulations, respectively. The OsHw distributions (Figures 15a and 16a) both show similar structure to that of formaldehyde (Figure 13a); a well-defined

first peak centered at 1.8 Å. However, for both acetaldehyde and acetone, there is a less well-defined second solvation peak in contrast to formaldehyde. The size and placement of the second peak also reflects the increasing size of the solutes (a more detailed analysis of this region of the g(R) curves was not performed). The integration of the first Os-Hw peak gives 1.8 and 2.0 for acetaldehyde and acetone, respectively, in comparison to 1.4 for formaldehyde. The stronger solvation across the series is further reflected in the relative heights of the peaks at 1.8 Å and the value of the minima at 2.5 Å in the three Os-Hw distributions (Figures 13a, 15a, and 16a). In the Os-Ow distributions, we see significant difference relative to the formaldehyde result. For both acetaldehyde and acetone, we observe a more well-defined first peak, centered at 2.8 Å (Figures 15b and 16b). Also, the 2.8 Å peak is slightly higher and the minimum at 3.0-3.25 Å is lower in the acetone curve, relative to the acetaldehyde curve indicating the stronger solvation for acetone. These 2.8 Å peaks integrate to 1.8 and 2.2 water molecules for acetaldehyde and acetone, respectively. Thus, we calculate the number of water molecules in the first solvation shell of the three solutes as follows: formaldehyde, 1.2; acetaldehyde, 1.8; and acetone, 2.2. This trend is reflected by the relative dipolar character of the three solutes. We calculate the simulation average ground-state dipole moment for the three solutes as follows: formaldehyde, 3.1 D; acetaldehyde, 3.7 D; and acetone, 4.0 D. This represents an induced dipole moment (relative to the gas-phase values above) as follows: formaldehyde, +0.77 D; acetaldehyde, +1.09 D; and acetone, +1.13 D. Both the total dipole moments and the solvent-induced dipole moments increase from formaldehyde to acetaldehyde and acetone. The trend in the substituent effect is reasonable since the methyl groups are better electron donors than hydrogen. Spectroscopy. The experimental gas-phase vertical n-π* excitation energies for the three solutes are as follows: formaldehyde, 32 828 cm-1; acetaldehyde, 34 522 cm-1; and acetone, 35 732 cm-1.60 This represents a gas-phase shift, relative to formaldehyde, of the following: formaldehyde, 0; acetaldehyde, +1694 cm-1; and acetone, +2904 cm-1. We optimized the gas-phase structure of the three solutes, using AM1, and calculated the n-π* excited state using AM1 SCI with an active space consisting of all single excitations from the three HOMOs into the four LUMOs (Methods). The resulting AM1 calculated n-π* gas-phase shift, relative to formaldehyde, is as follows: formaldehyde, 0; acetaldehyde, +1518 cm-1, and acetone, +3772 cm-1. The agreement with acetaldehyde is excellent, while for acetone we overestimate the relative shift by ∼800 cm-1. Yet, given the approximate nature of the AM1 Hamiltonian, we reproduce the overall experimental gas-phase trend well. In Table 3, we summarize the results of the spectroscopic calculations of the three solutes from the MD simulations along with some of the structural results given above. All three solutes

14506 J. Phys. Chem., Vol. 100, No. 34, 1996 show the expected blue shift upon solvation. The magnitude of the blue shift increases across the series as follows: formaldehyde, 1148 cm-1; acetaldehyde, 1530 cm-1; and acetone, 1607 cm-1. This trend is reasonable given the relative number of water molecules in the first solvation shell of the carbonyl oxygen for the three solutes (Figures 13-15). The acetone result is in reasonable accord with the experimental value of 1560 cm-1.58 Also, our acetone n-π* blue shift result compares well with a value of 1694 cm-1 reported by Gao14 and 1678 cm-1 reported by Debolt and Kollman.51 Our result for the formaldehyde n-π* blue shift is less than the value of 1900 cm-1 reported by Blair50 and considerably less than the value of 2557 cm-1 reported by FM.47 This trend reflects the relative number of waters reported in the first solvation shell of the formaldehyde carbonyl oxygen by us (1.2), Blair (2.6), and FM (6.2). If we assume our description of the blue shift trend for formaldehyde/acetaldehyde/acetone is correctly described by our QM/MMpol simulations, then our results would predict an experimental value of 1148 cm-1 for formaldehyde; within the predicted 600-1900 cm-1 given by Blair et al.50 Also, we predict the blue shift for acetaldehyde having a value of 1530 cm-1. There has been some discussion in the literature on the relative contributions of short- and long-range solvent effects on the n-π* blue shift of carbonyl chromophores.50 Thus, we recalculated the blue shift from the formaldehyde simulation but included only waters within 3.0 Å from any atom of formaldehyde. The value of the blue shift from this analysis is 560 cm-1. We predict that ∼50% of the blue shift is from shortrange and 50% from longer-range solvation for formaldehyde. Our results are in qualitative agreement with both Blair and FM, which indicate that cluster solvation is not sufficient to describe formaldehyde in solution, but rather longer-range solvent effects are important. VI. Summary and Conclusions We have presented the analytical gradient formalism for the QM/MMpol method, making this approach useful for MD simulations. Unlike previous approaches, we use the same Hamiltonian for both energy and gradients. Thus, our forces are consistent with the energy so that energy is conserved in MD simulations. Our approach to implementing the interaction of the QM electrons with the MM polarizable dipoles exploits the similarity of the treatment of two-center two-electron electrostatics in the MNDO approach with the classical interaction of polarizable dipoles. Indeed, we stress that one reason the semiempirical MNDO-type methods are appealing for QM/ MMpol is that they are very “molecular mechanics-like” in this regard. By simply addressing the consistent treatment of polarization in both the QM and MM parts, we are able to describe the interactions of CH2O/H2O well in comparison to accurate ab initio cluster calculations.52 QM/MMpol calculations show intriguing results for the CH2O/(H2O)3 cluster, which indicates it has a special stability in which CH2O inserts itself into the (H2O)3 cluster, allowing the remaining water-water H-bonds to obtain a more optimal linear structure. Indeed, the MD results at 150 K, show that the n ) 3 cluster has the largest blue shift in the vertical n-π* excited state of CH2O of all the cluster results. For clusters larger than n ) 3, we observe that CH2O binds to the outside of the cluster. The structural trends in the (H2O)n clusters that we calculate mirror more accurate MP2/ aug-cc-pVDZ results quite well. While a similar structural comparison for the optimized structures of the CH2O/(H2O)n

Thompson clusters is not possible, the structural and energetic results are reasonable and will hopefully serve as a basis for future comparison with more accurate QM studies. Using the QM/MM parameters developed for the cluster calculations, we presented results of aqueous condensed-phase simulations of formaldehyde, acetaldehyde, and acetone. The results showed a systematic trend in the solvation and spectroscopic properties of the three solutes. The first solvation shell of the carbonyl oxygen contains 1.2 (formaldehyde), 1.8 (acetaldehyde), and 2.2 (acetone) water molecules. We predict the vertical n-π* blue shift as 1148 cm-1 (formaldehyde), 1530 cm-1 (acetaldehyde), and 1607 cm-1 (acetone). The acetone result is in good accord with the experimental value of 1560 cm-1.58 Both the solvation and spectroscopy of these solutes give the expected trend with substituent on the carbonyl group. Acknowledgment. This work was performed under the auspices of the Division of Chemical Sciences, Office of Basic Energy Research, U.S. Department of Energy under Contract DE-AC06-76RLO 1830 with Battelle Memorial Institute, which operates the Pacific Northwest Laboratory, a multiprogram national laboratory. We gratefully acknowledge Gregory Schenter for a critical reading of this manuscript and useful discussions. References and Notes (1) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227. (2) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1986, 7, 718. (3) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700. (4) Luzhkov, V.; Warshel, A. J. Am. Chem. Soc. 1991, 113, 4491. (5) Luzhkov, V.; Warshel, A. J. Comput. Chem. 1992, 13, 199. (6) Gao, J.; Xia, X. Science 1992, 258, 631. (7) Gao, J. J. Phys. Chem. 1992, 96, 6432. (8) Gao, J.; Chou, L. W.; Auerbach, A. Biophys. J. 1993, 65, 43. (9) Stanton, R. V.; Hartsough, D. S.; Merz, K. M., Jr. J. Phys. Chem. 1993, 97, 11868. (10) Thompson, M. A.; Glendening, E. D.; Feller, D. F. J. Phys. Chem. 1994, 98, 10465. (11) Thompson, M. A. J. Phys. Chem. 1995, 99, 4794. (12) Thompson, M. A. J. Am. Chem. Soc. 1995, 117, 11341. (13) Thompson, M. A.; Schenter, G. K. J. Phys. Chem. 1995, 99, 6374. (14) Gao, J. J. Am. Chem. Soc. 1995, 117, 8600. (15) Muller, R. P.; Warshel, A. J. Phys. Chem. 1995, 99, 17516. (16) Hwang, J. K.; Pan, J. J. Chem. Phys. Lett. 1995, 243, 171. (17) Maseras, F.; Morokuma, K. J. Comput. Chem. 1995, 16, 1170. (18) Tunon, I.; Martins-Costa, T. C.; Millot, C.; Ruiz-Lopez, M. F.; Rivail, J. L. J. Comput. Chem. 1996, 17, 19. (19) Dang, L. X.; Rice, J. E.; Caldwell, J.; Kollman, P. A. J. Am. Chem. Soc. 1991, 113, 2481. (20) Dang, L. X. J. Chem. Phys. 1992, 97, 2659. (21) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757. (22) Bernardo, D. N.; Ding, Y.; Krogh-Jespersen, K.; Levy, R. J. Phys. Chem. 1994, 98, 4180. (23) Vesely, F. J. J. Comput. Phys. 1977, 24, 361. (24) Ahlstrom, P.; WallQvist, A.; Engstrom, S.; Jonsson, B. Mol. Phys. 1989, 68, 563. (25) Caldwell, J.; Dang, L. X.; Kollman, P. A. J. Am. Chem. Soc. 1990, 112, 9144. (26) Caldwell, J. W.; Kollman, P. A. J. Phys. Chem. 1995, 99, 6208. (27) Gao, J. J. Phys. Chem. 1995, 99, 16460. (28) Smith, D. E.; Dang, L. X. Chem. Phys. Lett. 1994, 230, 209. (29) Chang, T.; Peterson, K. A.; Dang, L. X. J. Chem. Phys. 1995, 103, 7502. (30) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 101, 7873. (31) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to AdVanced Electronic Structure Theory; MacMillian: New York, 1982. (32) Bottcher, C. J. F. The Theory of Electric Polarization; Elsevier: Amsterdam, 1973. (33) McCammon, J. A.; Harvey, S. C. Dynamics of proteins and nucleic acids; Cambridge University Press: Cambridge, U.K., 1987. (34) Dewar, M. J. S.; Thiel, W. Theor. Chim. Acta 1977, 46, 89. (35) Goldstein, H. Classical Mechanics; Addison-Wesley: Reading, MA, 1980. (36) Pople, J. A.; Beveridge, D. L. Approximate Molecular Orbital Theory; McGraw-Hill: New York, 1970.

QM/MMpol Model for Solute/Solvent Polarization (37) Thompson, M. A. Argus 3.0; Pacific Northwest Laboratory: Richland, WA, 1994 (99352, [email protected]). (38) Dewar, M. J. S.; Thiel, W. J. Am. Chem. Soc. 1977, 99, 4899. (39) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902. (40) Stewart, J. J. P. J. Comput. Chem. 1989, 10, 209. (41) Stewart, J. J. P. J. Comput. Chem. 1989, 10, 221. (42) Stewart, J. J. P. J. Comput. Chem. 1991, 12, 320. (43) Zerner, M. C.; Loew, G.; Kirchner, R.; Mueller-Westerhoff, U. J. Am. Chem. Soc. 1980, 102, 589. (44) Ridley, J.; Zerner, M. Theor. Chim. Acta 1973, 32, 111. (45) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (46) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (47) Fukunaga, H.; Morokuma, K. J. Phys. Chem. 1993, 97, 59. (48) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University Press: New York, 1987.

J. Phys. Chem., Vol. 100, No. 34, 1996 14507 (49) Blair, J. T.; Westbrook, J. D.; Levy, R. M.; Krogh-Jespersen, K. Chem. Phys. Lett. 1989, 154, 531. (50) Blair, J. T.; Krogh-Jespersen, K.; Levy, R. M. J. Am. Chem. Soc. 1989, 111, 6948. (51) Debolt, S. E.; Kollman, P. A. J. Am. Chem. Soc. 1990, 112, 7515. (52) Dimitrova, Y.; Peyerimhoff, S. D. J. Phys. Chem. 1993, 97, 12731. (53) Xantheas, S. S. J. Chem. Phys. 1994, 100, 7523. (54) Xantheas, S. S. J. Chem. Phys. 1995, 102, 4505. (55) Tsai, C. J.; Jordon, K. D. Chem. Phys. Lett. 1993, 213, 181. (56) Lovas, F. J. J. Phys. Chem. Ref. Data 1978, 7, 1445. (57) Liu, K.; Cruzan, J. D.; Saykally, R. J. Science 1996, 271, 929. (58) Bayliss, N. S.; McRae, E. G. J. Phys. Chem. 1954, 58, 1006. (59) Nelson, R. D.; Lide, D. R.; Maryott, A. A. National Standard Reference Data Set; National Bureau of Standards: Washington, DC, 1967. (60) Robin, M. B. Higher Excited States of Polyatomic Molecules; Academic Press: New York, 1985; Vol. III.

JP960690M

Solvent Polarization

Recommend Documents