Nuclear-Electronic Orbital Method within the Fragment Molecular

Oct 26, 2009 - Benjamin Auer, Michael V. Pak, and Sharon Hammes-Schiffer*. Department of Chemistry, 104 Chemistry Building, PennsylVania State ...
0 downloads 0 Views 176KB Size
Download PDF
5582

J. Phys. Chem. C 2010, 114, 5582–5588

Nuclear-Electronic Orbital Method within the Fragment Molecular Orbital Approach† Benjamin Auer, Michael V. Pak, and Sharon Hammes-Schiffer* Department of Chemistry, 104 Chemistry Building, PennsylVania State UniVersity, UniVersity Park, PennsylVania 16802 ReceiVed: July 28, 2009; ReVised Manuscript ReceiVed: September 3, 2009

The nuclear-electronic orbital (NEO) method, which treats selected nuclei quantum mechanically on the same level as the electrons, is combined with the fragment molecular orbital (FMO) method, which is an approximate scheme for electronic structure calculations of large systems. This FMO-NEO approach is implemented in conjunction with Hartree-Fock theory, density functional theory, and second-order perturbation theory using a variety of fragmentation schemes, including those that fraction covalent bonds. A multilayer FMO method, which limits the NEO calculation to a specified portion of the total system, is also implemented. These approaches are applied to four model systems: a water hexamer, a methyl-capped glycine dimer, a cluster of 32 water molecules, and a phenol molecule solvated by 16 water molecules. The FMO-NEO results are in excellent agreement with full NEO results for the calculation of properties associated with the nuclear quantum effects, such as zero-point energies, isotope effects, and vibrational excitation energies. The multilayer FMONEO method also provides reasonable results, and systematic improvements can be achieved by increasing the number of fragments in the NEO layer. I. Introduction Nuclear quantum effects impact a wide range of molecular properties, including energetics, vibrationally averaged geometries, and dynamics, as well as geometric and kinetic isotope effects. Conventional electronic structure calculations invoke the Born-Oppenheimer approximation, and the nuclei are treated as classical point charges that move on the electronically adiabatic ground-state potential energy surface. A variety of methods have been developed to include nuclear quantum effects without invoking the Born-Oppenheimer approximation.1-6 In the nuclear-electronic orbital (NEO) method,5 specified nuclei, typically hydrogen nuclei involved in hydrogen bonding or hydrogen transfer, are treated quantum mechanically on the same level as the electrons. The NEO method has been implemented in conjunction with Hartree-Fock theory (NEO-HF),5 multiconfigurational approaches,5,7,8 second-order perturbation theory [NEO-MP2(ee,ep)],9 which includes corrections for electronproton as well as electron-electron correlation, and density functional theory [NEO-DFT(ee)],10 which includes electronelectron correlation using standard electronic density functionals. To account for the significant electron-proton correlation, NEO methods using explicitly correlated wave functions11-13 and multicomponent DFT with electron-proton functionals derived from these explicitly correlated wave functions14 are under development. The extension of these NEO methods to larger systems requires fragment-based or hybrid quantum mechanical/molecular mechanical (QM/MM)15,16 methods. In fragment-based methods, the total system is divided into fragments, and the properties of the system are obtained from calculations on the constituent fragments. One successful fragment-based approach is the fragment molecular orbital (FMO) method.17,18 In the FMO method, electronic structure calculations are performed on the †

Part of the “Barbara J. Garrison Festschrift”. * To whom correspondence should be [email protected].

addressed.

E-mail:

fragments, which are immersed in the electrostatic field of the other fragments, according to a self-consistent-field procedure. The exchange interactions and correlation effects between fragments can be included by performing calculations on dimers, trimers, and higher-order fragment clusters, which are also immersed in the external electrostatic field of the other fragments. The FMO method has been shown to reproduce full electronic structure calculations to chemical accuracy when this series is expanded up to trimers.19,20 The FMO method has been implemented in conjunction with restricted Hartree-Fock (RHF),17,21 DFT,20,22,23 and correlated methods such as coupled cluster24 and MP225-27 theories. A multilayer FMO method has also been developed to allow the description of specified fragments at different levels of theory.28,29 The FMO framework is highly suitable for extending the NEO method to larger systems. In general, the nuclear wave functions and densities are much more localized than the electronic densities for molecular systems. As a result, the treatment of selected nuclei quantum mechanically does not contribute to the difficulties associated with fragmenting molecules at covalent bonds provided that these bonds do not involve a quantum nucleus. Moreover, exchange interactions between electrons and nuclei are not present, and correlation effects between electrons and nuclei are relatively short ranged compared to electronic correlation. In addition, nuclear-nuclear exchange and correlation effects are negligible in molecular systems due to the spatial localization of the nuclear densities. Thus, the quantum mechanical treatment of nuclei is not expected to significantly increase the errors introduced by the neglect of exchange and correlation effects between fragments. In many NEO applications, only a few nuclei are treated quantum mechanically, and a multilayer FMO approach may be used to limit the NEO calculations to a small subsystem. To our knowledge, only one study has combined the FMO method with nuclear wave functions for hydrogen nuclei, but this study was limited to a RHF formulation for the electrons and a single s-type Gaussian basis function for the quantum nuclei.30

10.1021/jp907193g  2010 American Chemical Society Published on Web 10/26/2009

NEO Method within the FMO Approach

J. Phys. Chem. C, Vol. 114, No. 12, 2010 5583

In this paper, we present a combined FMO-NEO approach at the NEO-HF, NEO-DFT(ee), and NEO-MP2(ee,ep) levels of theory. Nuclear basis sets containing s-, p-, and d-type Gaussian basis functions are utilized. A multilayer FMO-NEO approach, which enables the efficient treatment of a relatively small number of nuclei quantum mechanically, is also implemented. We apply these methods to four model systems: a water hexamer, a methyl-capped glycine dimer, a cluster of 32 water molecules, and a phenol molecule solvated by 16 water molecules. We compare the FMO-NEO and full NEO approaches for calculating the zero-point energies associated with the quantum nuclei, the differences between hydrogen and deuterium vibrational ground-state energies, and the vibrational excitation energies. An outline of the paper is as follows. Section II summarizes the theory underlying the FMO-NEO method and describes the computational details of our implementation. Section III presents the results and analyses of the applications to the four model systems. Conclusions are presented in Section IV.

F˜p,XCp,X ) Sp,XCp,Xεp,X

where S, F˜, C, and ε are the overlap, Fock, expansion coefficient, and eigenvalue matrices, respectively. The superscripts e and p denote electrons and quantum nuclei, respectively, and the superscript X indicates the monomer I or dimer IJ. The electronic and nuclear Hartree-Fock-Roothaan equations are solved iteratively to convergence. In eq 2, the electronic Fock operator is defined as

˜ e,X + Ge,X F˜e,X ) H

e,X e,X i e,X ˜ µν H ) Hµν + Vµν + B ∑ Pµν

EFMO )

N

∑ EI + ∑ (EIJ - EI - EJ) I

(1)

IJ

(12) Typically the third term in eq 5 is not necessary for quantum nuclei because the nuclear densities are so localized, but this term could be included if necessary. Thus, every fragment is immersed in the electrostatic potential arising from all classical nuclei, electrons, and quantum nuclei in the other fragments. In practice, the FMO-NEO-HF method is implemented as follows. First the monomer energies are converged selfconsistently with all other monomers at the NEO-HF level of theory if the monomer contains quantum nuclei or at the RHF level if the monomer contains no quantum nuclei. The external electrostatic potentials due to the other monomers are calculated using eqs 6 and 9. After convergence of the monomer energies, the dimer energies are calculated at the NEO-HF level of theory if the dimer contains quantum nuclei or at the RHF level if the dimer contains no quantum nuclei. The dimer calculations are performed only once in the electrostatic potential from the converged monomers and are not self-consistently converged with respect to the other dimers. Our implementation of the FMO method in the context of NEO-DFT(ee) includes electronic correlation using DFT with standard electronic functionals.10 In the FMO-NEO-DFT(ee) method, the Kohn-Sham equations for the electrons are solved iteratively with the NEO-HF equations for the quantum nuclei for each fragment. The external electrostatic potentials are defined to be the same as in the FMO-NEO-HF method, and the computational procedure is identical to that used for the FMO-NEO-HF method. Our implementation of the FMO method in the context of NEO-MP2(ee,ep) includes second-order perturbation theory corrections for both electron-electron and electron-proton correlation.9 After the monomers are converged at the NEOHF level, the NEO-MP2(ee,ep) corrections are calculated for all monomers and dimers to obtain the correlated monomer and dimer energies. The FMO-NEO-MP2(ee,ep) energy can be expressed as

EFMO-NEO-MP2 ) EFMO-NEO-HF + EFMO-corr

(10)

where N

EFMO-corr )

∑ I

N

Ecorr + I

corr - Ecorr ∑ (Ecorr IJ - EI J ) I>J

(11) and Ecorr and Ecorr I IJ are the monomer and dimer electron-electron and electron-proton NEO-MP2(ee,ep) correlation energies.

where LI is the layer to which monomer I belongs, and min(LI,LJ) is the lower level of theory for monomers I and J. The algorithm for a two-layer FMO calculation is as follows. First, all monomers are converged at the level of theory and basis set for the first layer (i.e., the lower level). Second, the monomers in the second layer are reconverged at the level of theory and basis set for the second layer, keeping the densities of the monomers belonging to the first layer fixed. The external electrostatic potential due to each monomer is calculated using the monomer density corresponding to the appropriate layer. Third, the dimer corrections are calculated at the lower level of the two monomers comprising each dimer. If desired, the monomer energies could be calculated self-consistently at the different levels of theory, thereby combining the first two steps of this algorithm. The combination of the multilayer FMO and NEO methods is particularly suitable for calculations in which only one or a few fragments contain quantum nuclei. In this case, we can implement a two-layer FMO-NEO-HF calculation, where the first layer is RHF and the second layer is NEO-HF. First all monomers are calculated at the RHF level. The quantum nuclei in the NEO fragment are treated as classical point charges for these RHF calculations. Second, the monomers in the second layer are reconverged at the NEO-HF level, treating the specified nuclei quantum mechanically in these monomers and fixing the densities of the other monomers. Third, the dimer calculations are performed at the appropriate level. Using this algorithm, many NEO-HF dimer calculations are omitted. When the second layer is comprised of only one monomer, all NEO-HF dimer calculations are avoided. We have implemented the FMO-NEO method at the NEOHF, NEO-DFT(ee), and NEO-MP2(ee,ep) levels in the GAMESS electronic structure package.32,33 In the FMO method, bonds are fractioned electrostatically, and electrons are divided in a 2:0 ratio such that both electrons in a bond are assigned to one of the two fragments. Projection operators are used to localize the electronic MOs within each fragment when fragmenting molecules at covalent bonds. A detailed description of techniques for fragmenting covalent bonds within the FMO approach is provided elsewhere.21,31 Recently a new fragmentation scheme suitable for solids has been developed,34 but we have implemented the NEO-FMO method only using the original fragmentation scheme.31 No special considerations concerning fragmenting molecules are required for the nuclear MOs because the nuclear densities are much more localized than the electronic densities, and a covalent bond with a quantum nucleus should not be fractioned (although in principle this is possible). Thus,

NEO Method within the FMO Approach

J. Phys. Chem. C, Vol. 114, No. 12, 2010 5585

Figure 1. The four systems examined in this study: (a) the water hexamer, (b) the glycine dimer, (c) the 32 water cluster, and (d) the phenol-water cluster. The fragmentation scheme used to partition the glycine dimer into three fragments is indicated in (b), where the dashed lines denote a fractioned covalent bond. In part (c), the blue oxygen identifies the water molecule with the quantum hydrogen nucleus, and the green oxygen atoms identify the two water molecules that are included as separate fragments in the second layer for 1W ML-3 and are included in the fragment containing the quantum nucleus for 3W and 3W ML-1. The quantum hydrogen nucleus in (c) is identified with an asterisk.

TABLE 1: Exponents of the DZSPDN-HD Nuclear Basis Set

bond fragmentation only impacts the quantum nuclei indirectly via the electronic density. As implemented in GAMESS, several approximations to the external electrostatic potentials and dimer energies can be used to enhance computational efficiency.35,36 In the present study, however, these options are not utilized (i.e., the GAMESS keywords RESPAP, RESPPC, and RESDIM are set to zero). The self-consistent-field and monomer convergences were set to 10-7, and a spherical harmonic basis was used. III. Results We performed calculations on four different systems to test the FMO-NEO method in a variety of situations. These four systems are depicted in Figure 1. First we examined a water hexamer and a methyl-capped glycine dimer that requires fragmenting at covalent bonds. In both of these studies, all hydrogen nuclei were treated quantum mechanically. Then we examined two larger systems, a cluster of 32 waters and a phenol molecule solvated by 16 waters, to test the multilayer FMO method. In these two studies, only one hydrogen nucleus was treated quantum mechanically. All of these calculations were performed with a modified version of GAMESS. In all cases, the 6-31G(d,p) electronic basis set was used.37-40 For the water hexamer and the glycine dimer, the DZSPDN nuclear basis set developed for hydrogen nuclei was used.5 The DZSPDN basis set contains two each of s-type, p-type, and d-type Gaussian functions. For the 32 water cluster and the phenol-water system, the DZSPDN basis set developed for hydrogen was combined with a DZSPDN basis set developed

DZSPDN-HD s-type s-type s-type p-type p-type p-type d-type d-type d-type

42.32953 30.60121 20.17583 40.25819 29.49726 18.21766 36.88654 25.51488 16.73036

for deuterium to facilitate the study of isotope effects.41 Since some of the exponents in the hydrogen and deuterium basis sets are similar, one each of the s-type, p-type, and d-type functions were eliminated, leading to a basis set with three each of s-type, p-type, and d-type functions. This nuclear basis set, denoted as DZSPDN-HD, is given in Table 1. The first system we examined was a water hexamer in the prism structure that was optimized at the RHF/6-31G(d,p) level of theory. Each fragment contained one water molecule, leading to six monomers, and all 12 hydrogen nuclei were treated quantum mechanically in the NEO calculations. These calculations were performed at the HF, DFT/B3LYP,42-44 and MP2(ee,ep) levels. Table 2 provides the values for the HF and DFT/B3LYP energies with all nuclei treated classically and with all hydrogen nuclei treated quantum mechanically, as well as ∆EQN, which corresponds to the zero-point energy associated with the quantum nuclei. At the HF level, the error in the absolute energy introduced by using the FMO method is ∼0.0025 hartree for both the classical and quantum treatment

5586

J. Phys. Chem. C, Vol. 114, No. 12, 2010

Auer et al.

TABLE 2: Energies for the Water Hexamer at the Specified Level of Theorya method

Eclb

Equc

∆EQNd

HF FMO/HF DFT B3LYP FMO/DFT B3LYP

-456.218396 -456.220940 -458.389890 -458.393542

-455.727655 -455.730412 -457.907326 -457.911221

0.490741 0.490528 0.482564 0.482321

a Energies given in atomic units. FMO results include dimer calculations. b Ecl is the energy with all protons treated as classical point charges. c Equ is the energy with all protons treated quantum mechanically using NEO. d ∆EQN is defined as the difference between energies obtained with classical protons and quantum protons.

TABLE 3: Energies for the Glycine Dimer at the Specified Level of Theorya method

Eclb

Equc

∆EQNd

HF FMO/HF DFT B3LYP FMO/DFT B3LYP

-660.666496 -660.667792 -664.182803 -664.186589

-660.106024 -660.107285 -663.627077 -663.629634

0.560472 0.560507 0.555726 0.556955

a

Energies given in atomic units. FMO results include dimer calculations. b Ecl is the energy with all protons treated as classical point charges. c Equ is the energy with all protons treated quantum mechanically using NEO. d ∆EQN is defined as the difference between energies obtained with classical protons and quantum protons.

of the nuclei. At the DFT level, the error in the absolute energy is slightly greater. The error in ∆EQN introduced by using the FMO method is ∼0.0002 hartree at both the HF and DFT levels. Thus, although a significant error is introduced in the absolute energy, the zero-point energies associated with the quantum nuclei are calculated with reasonable accuracy. We also determined the NEO-MP2(ee,ep) energy of this system by calculating the second-order corrections for the monomers and the dimers after the monomers were converged self-consistently at the HF level. The electron-electron and electron-proton correlation energies calculated with the FMO-NEO-MP2(ee,ep) method are 1.191012 and 0.051169 hartree, respectively. These values are in excellent agreement with the full NEO-MP2(ee,ep) results of 1.190993 and 0.051180 hartree, respectively. The agreement is slightly better for the electron-proton correlation. The second system we examined was the methyl-capped glycine dimer that was also optimized at the RHF/6-31G(d,p) level. All hydrogen nuclei were treated quantum mechanically for the NEO calculations. The glycine dimer was divided into three fragments at the R-carbons, as depicted in Figure 1b. As mentioned above, covalent bonds are fractioned so that one fragment receives both electrons in the bond while the other fragment receives none. For this system, both electrons associated with the bond fractioned between the CH2 and CO groups were assigned to the fragment containing the CO, as shown in Figure 1b. The same series of calculations that were performed for the water hexamer were also performed for this system. The results are presented in Table 3. Once again, the error introduced into the absolute energy by using the FMO method for either classical or quantum protons at the HF and DFT levels is on the order of 10-3 hartree. The error in ∆EQN introduced by using the FMO method is 0.00003 hartree at the HF level and 0.0012 hartree at the DFT level. In addition, we calculated the NEOMP2(ee,ep) correlation energy of this system. The electronelectron and electron-proton correlation energies calculated with the FMO-NEO-MP2(ee,ep) method are 1.962119 and 0.057635 hartree, respectively. These values are in excellent

TABLE 4: Energies for the 32 Water Cluster at the Specified Level of Theorya method HF FMO/HF FMO/HF FMO/HF FMO/HF FMO/HF

1W 3W 1W ML-1 1W ML-3 3W ML-1

Eclb

Equc

∆EQNd

-2433.128602 -2433.150290 -2433.148037

-2433.087864 -2433.109594 -2433.107322 -2433.108476 -2433.109101 -2433.106828

0.040738 0.040696 0.040715 0.041814 0.041189 0.041209

a Energies given in atomic units. FMO results include dimer calculations. b Ecl is the energy with all protons treated as classical point charges. c Equ is the energy with a single proton treated quantum mechanically using NEO. d ∆EQN is defined as the difference between energies obtained with classical protons and quantum protons.

agreement with the full NEO-MP2(ee,ep) results of 1.961845 and 0.057637 hartree, respectively. Once again, the agreement is slightly better for the electron-proton correlation. The third system we examined was a cluster of 32 water molecules. This structure was generated by Fedorov and Kitaura19 and was obtained from the official FMO Web site.45 In these calculations, only one hydrogen nucleus on a water molecule in the middle of the cluster was treated quantum mechanically, while the other nuclei were treated as classical point charges. All calculations for this cluster were performed at the HF level. Our objective was to test the multilayer FMO approach, as well as to analyze the impact of the size of the fragment containing the quantum proton on the accuracy. The results are summarized in Table 4. When each fragment contains one water molecule, denoted as FMO/HF 1W in Table 4, the error introduced by using the FMO method is ∼0.02 hartree for the absolute energies but is less than 10-4 hartree for ∆EQN. Previously the error in the FMO method for absolute energies in electronic structure calculations was found to increase approximately linearly with system size,31 but the energetic properties associated with the nuclear quantum effects do not appear to follow this trend. The accuracy is improved further when the water molecule containing the quantum nucleus is combined with the two water molecules containing the oxygen atoms closest to the quantum nucleus into a single fragment, denoted as FMO/HF 3W in Table 4. We also applied the two-layer FMO method to this system, where the first layer is RHF and the second layer is NEO-HF. At the FMO/HF 1W ML-1 level, each fragment includes one water molecule, and the second layer includes only the water molecule with the quantum nucleus. In this case, the error in ∆EQN introduced by the FMO method increases to ∼0.001 hartree. The accuracy is improved by including more water molecules in the second layer. When the two water molecules with oxygen atoms closest to the quantum nucleus are included in the second layer but are still treated as separate fragments, denoted FMO/HF 1W ML-3 in Table 4, the error in ∆EQN is ∼0.0004 hartree. A similar error in ∆EQN is obtained when the fragment with the quantum nucleus includes the three water molecules, all other fragments include only one water, and only the fragment containing the quantum nucleus is included in the second layer, denoted FMO/HF 3W ML-1 in Table 4. The errors introduced by the multilayer approximation within the FMO approach arise mainly from the calculation of the dimer energies at the lower level of theory, rather than from the multilayer procedure for calculating the monomer energies. In addition, we examined properties such as isotope effects and vibrational excitation energies using the FMO method. To

NEO Method within the FMO Approach

J. Phys. Chem. C, Vol. 114, No. 12, 2010 5587

TABLE 5: Energies for the 32 Water Cluster at the Specified Level of Theorya method HF FMO/HF FMO/HF FMO/HF FMO/HF FMO/HF

1W 3W 1W ML-1 1W ML-3 3W ML-1

TABLE 6: Energies for the Phenol-Water Cluster at the Specified Level of Theorya

Hb Egs

∆EHDc

∆E12d

-2433.087864 -2433.109594 -2433.107322 -2433.108476 -2433.109101 -2433.106828

0.011309 0.011296 0.011302 0.011647 0.011450 0.011458

0.022205 0.022184 0.022199 0.022986 0.022492 0.022511

a

Energies given in atomic units. FMO results include dimer H calculations. b Egs is the ground-state energy with a single hydrogen nucleus treated quantum mechanically using NEO. c ∆EHD is the energy difference between the systems with a quantum hydrogen and a quantum deuterium. d ∆E12 is the energy difference between the NEO first excited and ground states.

examine isotope effects, the calculations were performed by replacing the quantum hydrogen nucleus with deuterium. The results are summarized in Table 5. The error in ∆EHD, the difference between the hydrogen and deuterium vibrational ground-state energies, introduced by using the FMO method was ∼10-5 hartree when one or three water molecules were included in the fragment with the quantum nucleus. The error in ∆EHD introduced by using two-layer FMO methods that include three water molecules in the NEO layer was ∼0.00015 hartree. To examine the vibrational excitation energy ∆E12, we placed the single quantum nucleus in the second nuclear orbital (i.e., the orbital with the second lowest energy) when solving the Hartree-Fock-Roothaan equations. Note that this procedure is valid only for a single quantum nucleus, where the excited vibrational states correspond to the higher-energy nuclear orbitals. When multiple nuclei are treated quantum mechanically, other methods such as multiconfigurational approaches are required to calculate accurate vibrational excitation energies. As shown in Table 5, the errors in ∆E12 introduced by using the FMO and two-layer FMO methods are similar to the corresponding errors in ∆EHD. The fourth system we examined was a phenol molecule solvated by 16 water molecules. This structure was generated by Fedorov and Kitaura29 and was obtained from the official FMO Web site.45 In these calculations, the hydrogen nucleus on the hydroxyl group of the phenol was treated quantum mechanically, while all other nuclei were treated as classical point charges. All calculations for this system were performed at the HF level. Two different fragmentation schemes were implemented. In the first scheme, denoted 1WPH, the phenol and the water closest to the quantum nucleus were treated as a single fragment, and the other fragments each contained one water molecule. In the second scheme, denoted 1WOH, only the OH group of the phenol and the water closest to the quantum nucleus were treated as a single fragment, and the other fragments each contained one water molecule or the phenol ring. In this case, the C-O bond between the ring and the oxygen atom was fractioned so that both electrons in the C-O bond were included in the fragment containing the OH group. The results for the phenol-water system are summarized in Table 6. The error in ∆EQN introduced by using the full FMO method with either the 1WPH or 1WOH fragmentation scheme is on the order of 10-5 hartree. The error in ∆EQN introduced by using the two-layer FMO method in which only the fragment with the quantum nucleus is included in the NEO layer increases to ∼10-3 hartree for both schemes. The error is larger for the 1WOH scheme than for the 1WPH scheme because of the fractioned covalent bond. Additional errors are introduced for the 1WOH scheme because the bond is fragmented at an sp2-

method HF FMO/HF FMO/HF FMO/HF FMO/HF

1WPH 1WOH 1WPH ML 1WOH ML

Eclb

Equc

∆EQNd

-1522.001180 -1522.008593 -1522.007914

-1521.960801 -1521.968233 -1521.967588 -1521.967092 -1521.964967

0.040379 0.040360 0.040326 0.041501 0.042947

a

Energies given in atomic units. FMO results include dimer calculations. b Ecl is the energy with all protons treated as classical point charges. c Equ is the energy with the phenol hydroxyl proton treated quantum mechanically using NEO. d ∆EQN is defined as the difference between energies obtained with classical protons and quantum protons.

TABLE 7: Energies for the Phenol-Water Cluster at the Specified Level of Theorya method HF FMO/HF FMO/HF FMO/HF FMO/HF

1WPH 1WOH 1WPH ML 1WOH ML

Hb Egs

∆EHDc

∆E12d

-1521.960801 -1521.968233 -1521.967588 -1521.967092 -1521.964967

0.011198 0.011193 0.011187 0.011495 0.012390

0.021854 0.021847 0.021845 0.022359 0.026101

a Energies given in atomic units. FMO results include dimer H calculations. b Egs is the ground-state energy with the phenol hydroxyl hydrogen nucleus treated quantum mechanically using NEO. c ∆EHD is the energy difference between the systems with a quantum hydrogen and a quantum deuterium at the phenol hydroxyl position. d ∆E12 is the energy difference between the NEO first excited and ground states.

hybridized carbon atom and the OH-water fragment is unconventional. The accuracy of the multilayer approach can be improved by increasing the size of the fragment containing the quantum nucleus. For example, including the phenol and two water molecules in the NEO fragment reduces the error in ∆EQN to 0.0006 hartree. We also examined the isotope effects and vibrational excitation energies for this system. The results are summarized in Table 7. The FMO method provides a highly accurate description of ∆EHD and ∆E12. The two-layer FMO method is not as accurate, and again the error is larger for the 1WOH scheme due to fragmentation of the covalent bond. As discussed previously,11-14,46 the NEO-HF, NEO-DFT(ee), and NEO-MP2(ee,ep) methods do not provide accurate hydrogen vibrational frequencies due to inadequate treatment of electronproton correlation. Specifically, the hydrogen vibrational frequencies tend to be too high in the absence of sufficient electron-proton correlation. As a result, the values of ∆EQN, ∆EHD, and ∆E12 calculated with these methods are correspondingly high. Future work will be aimed at including sufficient electron-proton correlation in conjunction with the FMO method using the explicitly correlated Hartree-Fock (NEOXCHF) method11-13 and the multicomponent DFT(ee,ep) method14 with electron-proton functionals derived from explicitly correlated wave functions. These approaches will provide a more quantitatively accurate description of the nuclear wave functions and associated hydrogen vibrational frequencies. IV. Conclusions In this paper, we combined the NEO method, which treats selected nuclei quantum mechanically on the same level as the electrons, with the FMO method. We applied this FMO-NEO approach to four model systems at the NEO-HF, NEO-DFT(ee), and NEO-MP2(ee,ep) levels of theory using a variety of fragmentation schemes, including those that fraction covalent

5588

J. Phys. Chem. C, Vol. 114, No. 12, 2010

bonds. We found that the FMO-NEO results are in excellent agreement with full NEO results for the calculation of properties associated with the nuclear quantum effects. We also implemented a multilayer FMO method in which the NEO calculation is performed for only a specified portion of the total system. We found that the multilayer method provides reasonable results, and systematic improvements to the accuracy can be achieved by increasing the number of fragments in the NEO layer. The FMO-NEO method can be extended in a variety of directions. Most importantly, the quantitatively accurate calculation of nuclear wave functions requires an adequate treatment of electron-proton correlation. For this purpose, we will implement the FMO method in conjunction with the NEOXCHF11-13 and NEO-DFT(ee,ep)14 methods, which have been shown to include sufficient electron-proton correlation to provide accurate nuclear densities. In addition, many of the technical developments that have been applied to the electronic structure FMO method can also be applied to the FMO-NEO method. These technical developments include approximations to the external electrostatic potentials to enhance the computational efficiency35,36 and alternative fragmentation schemes to minimize the errors associated with fragmenting molecules at covalent bonds.34 Since the nuclear densities are much more localized than the electronic densities, the quantum mechanical treatment of nuclei is not expected to lead to additional errors associated with these approximations. Furthermore, multilayer FMO-NEO approaches involving the mixture of various levels of theory for both the NEO and the conventional electronic structure layers will also be investigated. These extensions will enable the FMO-NEO approach to be utilized to treat the key hydrogen nuclei quantum mechanically in electronic structure calculations on large molecular systems. Acknowledgment. We thank Anirban Hazra for helpful discussions. We gratefully acknowledge the support of AFOSR Grant FA9550-07-1-0143. References and Notes (1) Shigeta, Y.; Ozaki, Y.; Kodama, K.; Nagao, H.; Kawabe, H.; Nishikawa, K. Int. J. Quantum Chem. 1998, 69, 629. (2) Tachikawa, M.; Mori, K.; Nakai, H.; Iguchi, K. Chem. Phys. Lett. 1998, 290, 437. (3) Nakai, H.; Sodeyama, K.; Hoshino, M. Chem. Phys. Lett. 2001, 345, 118. (4) Kreibich, T.; Gross, E. K. U. Phys. ReV. Lett. 2001, 86, 2984. (5) Webb, S. P.; Iordanov, T.; Hammes-Schiffer, S. J. Chem. Phys. 2002, 117, 4106. (6) Bochevarov, A. D.; Valeev, E. F.; Sherrill, C. D. Mol. Phys. 2004, 102, 111. (7) Pak, M. V.; Hammes-Schiffer, S. Phys. ReV. Lett. 2004, 92, 103002. (8) Skone, J. H.; Pak, M. V.; Hammes-Schiffer, S. J. Chem. Phys. 2005, 123, 134108. (9) Swalina, C.; Pak, M. V.; Hammes-Schiffer, S. Chem. Phys. Lett. 2005, 404, 394. (10) Pak, M. V.; Chakraborty, A.; Hammes-Schiffer, S. J. Phys. Chem. A 2007, 111, 4522.

Auer et al. (11) Swalina, C.; Pak, M. V.; Chakraborty, A.; Hammes-Schiffer, S. J. Phys. Chem. A 2006, 110, 9983. (12) Chakraborty, A.; Pak, M. V.; Hammes-Schiffer, S. J. Chem. Phys. 2008, 129, 014101. (13) Chakraborty, A.; Hammes-Schiffer, S. J. Chem. Phys. 2008, 129, 204101. (14) Chakraborty, A.; Pak, M. V.; Hammes-Schiffer, S. Phys. ReV. Lett. 2008, 101, 153001. (15) Gao, J.; Truhlar, D. G. Quantum mechanical methods for enzyme kinetics. In Annual ReViews of Physical Chemistry; Annual Reviews, Inc.: Palo Alto, CA, 2002; Vol. 53, p 467. (16) Friesner, R. A.; Guallar, V. Annu. ReV. Phys. Chem. 2005, 56, 389. (17) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Chem. Phys. Lett. 1999, 313, 701. (18) Fedorov, D. G.; Kitaura, K. J. Phys. Chem. A 2007, 111, 6904. (19) Fedorov, D. G.; Kitaura, K. J. Chem. Phys. 2004, 120, 6832. (20) Fedorov, D. G.; Kitaura, K. Chem. Phys. Lett. 2004, 389, 129. (21) Nakano, T.; Kaminuma, T.; Sato, T.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Chem. Phys. Lett. 2000, 318, 614. (22) Shimodo, Y.; Morihashi, K.; Nakano, T. J. Mol. Struct.: THEOCHEM 2006, 770, 163. (23) Sugiki, S.-I.; Kurita, N.; Sengoku, Y.; Sekino, H. Chem. Phys. Lett. 2003, 382, 611. (24) Fedorov, D. G.; Kitaura, K. J. Chem. Phys. 2005, 123, 134103. (25) Kochizuki, Y.; Koikegami, S.; Nakano, T.; Amari, S.; Kitaura, K. Chem. Phys. Lett. 2004, 396, 473. (26) Fedorov, D. G.; Kitaura, K. J. Chem. Phys. 2004, 121, 2483. (27) Fedorov, D. G.; Ishimura, K.; Ishida, T.; Kitaura, K.; Pulay, P.; Nagase, S. J. Comput. Chem. 2007, 28, 1476. (28) Fedorov, D. G.; Ishida, T.; Kitaura, K. J. Phys. Chem. A 2005, 109, 2638. (29) Fedorov, D. G.; Kitaura, K. J. Chem. Phys. 2005, 122, 054108. (30) Ishimoto, T.; Tachikawa, M.; Nagashima, U. J. Chem. Phys. 2006, 124, 014112. (31) Fedorov, D. G.; Kitaura, K. In Modern Methods for Theoretical Physical Chemistry and Biopolymers; Starikov, E. B., Lewis, J. P., Tanaka, S., Eds.; Elsevier: Amsterdam, 2006; p 3. (32) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347. (33) Gordon, M. S.; Schmidt, M. W. Advances in electronic structure theory: GAMESS a decade later. In Theory and Applications of Computational Chemistry: The First Forty Years; Dykstra, C. E., Frenking, G., Kim, K. S., Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005. (34) Fedorov, D. G.; Jensen, J. H.; Deka, R. C.; Kitaura, K. J. Phys. Chem. A 2008, 112, 11808. (35) Nakano, T.; Kaminuma, T.; Sato, T.; Fukuzawa, K.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Chem. Phys. Lett. 2002, 351, 475. (36) Fedorov, D. G.; Kitaura, K. Chem. Phys. Lett. 2006, 433, 182. (37) Ditchfield, R.; Hehre, W. J.; Pople, J. A. J. Chem. Phys. 1971, 54, 724. (38) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257. (39) Hariharan, P. C.; Pople, J. A. Theor. Chim. Acta 1973, 28, 213. (40) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. J. Chem. Phys. 1982, 77, 3654. (41) Reyes, A.; Pak, M. V.; Hammes-Schiffer, S. J. Chem. Phys. 2005, 123, 064104. (42) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (43) Stephens, P. J.; Devlin, F. J.; Chablowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623. (44) Hertwig, R. H.; Koch, W. Chem. Phys. Lett. 1997, 268, 345. (45) http://staff.aist.go.jp/d.g.fedorov/fmo/main.html, 2009. (46) Swalina, C.; Hammes-Schiffer, S. J. Phys. Chem. A 2005, 109, 10410.

JP907193G