Improved Complete Active Space Configuration Interaction Energies

Article pubs.acs.org/JCTC

Improved Complete Active Space Configuration Interaction Energies with a Simple Correction from Density Functional Theory Shiela Pijeau† and Edward G. Hohenstein*,†,‡ †

Department of Chemistry and Biochemistry, The City College of New York, New York, New York 10031, United States Ph.D. Program in Chemistry, The City University of New York, New York, New York 10016, United States

‡

S Supporting Information *

ABSTRACT: Recent algorithmic advances have extended the applicability of complete active space configuration interaction (CASCI) methods to molecular systems with hundreds of atoms. While this enables simulation of photochemical dynamics in the condensed phase, the underlying CASCI method has some well-known problems resulting from a severe neglect of dynamic electron correlation. Vertical excitation energies, vibrational frequencies, and reaction barriers are systematically overestimated; these errors limit the applicability of CASCI. We develop a correction for the CASCI energy using density functional theory (DFT). The DFT correction incorporates the effect of dynamic electron correlation among the core electrons into the CASCI Hamiltonian. We show that the resulting DFT-corrected CASCI approach is applicable in situations where the usual single-reference DFT methods fail, such as the description of systems with biradicaloid electronic structure and conical intersections between ground and excited electronic states. Finally, we apply this DFT-corrected CASCI approach to ultrafast excited-state proton transfer dynamics. Without the DFT correction, CASCI predicts spurious reaction barriers to these processes, and, as a result, a qualitatively correct description of the dynamics is not possible. With the DFT-corrected CASCI method, we demonstrate qualitative and quantitative agreement with both theory and experiment for two model systems for excited-state intramolecular proton transfer. Finally, we apply the DFT-corrected CASCI method to excitedstate proton transfer dynamics in a system with more than 150 atoms.

1. INTRODUCTION The efficient determination of reliable excited-state potential energy surfaces is one of the great remaining challenges in electronic structure theory. This is key to the simulation of photochemical processes where there is a dearth of generally applicable methods.1−6 Ideally, one would like to be able to perform dynamical simulations of photochemical processes in condensed phase systems. Situations where the environment surrounding a chromophore is an active participant in the photochemistry (e.g., excited-state proton transfer) are of particular interest.7,8 However, this imparts a unique set of requirements on an electronic structure method. The method must be scalable to hundreds of atoms, while remaining efficient enough to perform dynamical simulations. The method must also provide a robust description of large regions of multiple potential energy surfaces. In general, this might include chemical reactions and degeneracies between electronic states (conical intersections). Finally, for quantitative predictions of experimental observables, the method must also be able to accurately predict excitation energies and reaction barriers. While there are suitable methods available for most ground state processes, there is no one method that satisfies all of these requirements for excited-state chemistry. The most widely applied excited-state methods are single reference methods based on a single Slater determinant © 2017 American Chemical Society

reference state. This includes time-dependent density functional theory (TDDFT),9,10 configuration interaction singles (CIS),11 or coupled-cluster (CC) based approaches (in the equation of motion (EOM) or linear response formalisms).12−15 In general, single reference methods suffer from two major problems. These methods tend to perform well near equilibrium geometries; however, during the course of a photochemical process, regions of the potential surface may become relevant where the description of the ground state wave function becomes poor. Since the treatment of the excited state is intrinsically tied to the treatment of the ground state, the accuracy of the excited-state potential surface might degrade. Another problem is that many of these methods cannot describe conical intersections between electronic states.16−18 In some common single reference methods, the ground state is completely decoupled from the excited states; as a result, the topology of the potential energy surface around points of degeneracy is incorrect. One solution to these problems is to apply methods that define the reference state to have a different spin or a different number of electrons from the states of interest. These are the spin-flip19−21 or EOM ionization potential (EOM-IP) and EOM electron affinity (EOM-EA) Received: September 9, 2016 Published: February 3, 2017 1130

DOI: 10.1021/acs.jctc.6b00893 J. Chem. Theory Comput. 2017, 13, 1130−1146

Article

Journal of Chemical Theory and Computation methods.22,23 Since the ground state is not the reference state, the problems mentioned above are avoided. However, in the context of CIS and TDDFT, the usual formulation of spin-flip breaks spin-symmetry and introduces spurious mixing of singlet and triplet states, for example.24 Recently, exciting progress has been made to remedy the issues that have plagued spin-flip TDDFT.25 The EOM-CC approaches are also very promising;26 however, they remain too computationally demanding to apply in dynamical simulations of most photochemical processes. Ideally, one would apply methods with a true multiconfigurational reference state such as multireference configuration interaction (MRCI)27 or multireference perturbation theory.28 By using a multiconfigurational reference, these methods are typically well behaved over larger regions of the potential energy surface; methods that include multiple electronic states in the reference avoid bias toward one particular state, which improves the description of excited states. When multireference methods are applied in dynamical simulations of photochemical processes, the results can be quite accurate.29−36 Unfortunately, these methods are prohibitively expensive for most chemical systems, and as a result their range of application is severely limited. In the context of complete active space second order perturbation theory (CASPT2), exciting progress has been made in the development of efficient implementations of the analytic gradients that promise to extend the applicability of this method in nonadiabatic dynamics simulations.37−39 It is also possible to reduce the cost of multireference configuration interaction and multireference perturbation theory computations through the use of restricted active spaces (RAS)40,41 or general active spaces (GAS).42,43 Due to the problems incurred by single reference methods and the computational expense of multireference methods, more approximate methods such as complete active space configuration interaction (CASCI) or the complete active space self-consistent field (CASSCF) method are more commonly applied.44,45 These methods typically provide a consistent description of multiple potential surfaces at a wide range of geometries and are well behaved around conical intersections. The primary problem with methods based on a CASCI wave function is that they neglect dynamic electron correlation. Typically, this results in blue-shifted vertical excitation energies and an overestimation of reaction barriers. In many cases, it is possible to extract a robust description of photochemical processes in spite of these shortcomings.2−5,46,47 One drawback of conventional CASCI methods is the factorial scaling of the CI problem as the size of the active space increases. As a result, solving the problems mentioned above by enlarging the active space is often infeasible due to the immense computational demands. Alternatives to the CASCI wave function such as the density matrix renormalization group,48 two-electron reduced density matrix methods,49,50 or stochastic approaches to the construction of the multiconfigurational wave function51,52 offer the possibility to use much larger active spaces. A common approach to correct these issues is to use a semiempirical parametrized Hamiltonian,53−55 which can be tuned for a specific problem of interest.56−59 While the semiempirical approach can be very fruitful, the results are highly dependent on the parameters that have been chosen, and those parameters might not be transferable between systems. Recently, we have developed an algorithm that exploits sparsity in the atomic orbital basis to allow CASCI-type wave functions (and their derivatives) to be evaluated with

computational effort scaling quadratically, 6(N 2), with system size (assuming a constant number of active orbitals).60 With this approach it is easily possible to determine CASCI-type wave functions in systems that include more than 1000 atoms and 10000 basis functions. This algorithm has been applied to the gradients and nonadiabatic coupling vectors of stateaveraged CASSCF wave functions.61 Recently, this CASSCF program has been interfaced with the CI program of Fales and Levine62 to allow for the use of arbitrary active spaces in these state-averaged CASSCF gradient computations.63 We have also developed an implementation of floating occupation molecular orbital CASCI (FOMO-CASCI) that uses the same sparsityexploiting algorithm;64,65 FOMO-CASCI is a popular alternative to CASSCF for simulations of photochemical processes.56−59,66−69 We have leveraged the FOMO-CASCI program to perform molecular dynamics simulations of photoexcited aqueous methyl viologen that included between 300 and 600 atoms in the CASCI computation.70 The algorithm we have developed is ideally suited for computations of a small chromophore embedded in an insulating environment: this could be a dye molecule in solution or the chromophore of a photoactive protein. The ability to treat the surrounding environment quantum mechanically in a single, large CASCI computation is especially useful if the surroundings are an active participant in the photochemistry (i.e., excited-state proton transfer). Our CASCI algorithm allows the atomic and molecular orbitals to be efficiently manipulated; however, the method still scales factorially as the number of active orbitals are increased. As a result, it is best suited for systems that can be described by compact CASCI wave functions; fortunately, modest active spaces are often successful in simulations of photochemical processes. The case of a dye molecule in solution clearly illustrates the primary shortcoming of the CASCI approach: in order to treat a few electrons exactly, the rest of the system is treated under a meanfield approximation (i.e., Hartree−Fock-like). Our goal is to improve upon this mean-field treatment of the inactive electrons, while preserving the desirable properties of the CASCI wave function. In the present work, we address some of the shortcomings of CASCI-based methods by introducing a treatment of dynamic electron correlation into the Hamiltonian. Motivated by embedding approaches, we introduce a correction from density functional theory (DFT) that accounts for the missing electron correlation between core electrons. We demonstrate that this correction ameliorates some of the problems exhibited by CASCI methods, while not damaging the attractive features of the method. We test this DFT-corrected CASCI method for its ability to describe the reaction coordinate of two prototypical excited-state intramolecular proton transfer systems. Finally, we show that this method is both sufficiently accurate and computationally efficient to provide an accurate treatment of excited-state proton transfer dynamics for these two molecules as well as an aqueously solvated system.

2. THEORY The complete active space configuration interaction (CASCI) method partitions the molecular orbitals (indexed p, q, r, s) into core orbitals (indexed i, j, k, l), active orbitals (indexed t, u, v, w, ...), and strictly unoccupied virtual orbitals (indexed a, b, c, d). A set of Slater determinants is generated by taking all possible combinations of active electrons in active orbitals; the core 1131


Article

Journal of Chemical Theory and Computation Ecore = 2 ∑ (i|h|̂ i) +

orbitals remain occupied in all of the determinants. It is convenient to start from a single determinant describing the core electrons, viz. α Ncore

β Ncore

i

j

∏ αî † ∏

|Φcore⟩ =

i

|ΦI ⟩ =

∏

αt̂ †

t

β Nact

̂ |q) = (p|h|̂ q) + (p|hcore

(1)

The Hamiltonian is diagonalized in the basis of these determinants producing the CASCI wave function as a linear combination of Slater determinants |ΨCI⟩ =

∑ cI|ΦI ⟩

(3)

I

where cI is an eigenvector of the Hamiltonian. The CASCI method treats electron−electron interactions exactly within the active space and in a mean-field picture outside of the active space. In the limit where the active space includes all molecular orbitals, CASCI is equivalent to full CI. As a result of the exact treatment of the electrons within the active space, CASCI performs well in particularly challenging situations. For example, CASCI is able to capture the multiconfigurational character of wave functions that is present during bond dissociation and at intersections between electronic states (i.e., conical intersections). This allows CASCI to provide a consistent description of multiple electronic states, simultaneously, over large regions of their potential energy surfaces. From a theoretical perspective, CASCI has the attractive property that the energy is invariant to orbital rotations within an orbital space (i.e., core−core, active−active, and virtual−virtual rotations). However, the mean-field treatment of the core electrons results in a nearcomplete neglect of dynamic electron correlation (when only a few electrons are active). The resulting errors manifest in many situations; particularly troubling is the systematic overestimation of vertical excitation energies and, often, reaction barriers. Although CASCI wave functions can provide accurate descriptions of a variety of chemical problems, these systematic errors limit the applicability of the method and necessitate thorough benchmarking against more reliable methods. We will start from a general expression for the electronic Hamiltonian in terms of second-quantized operators Ĥ =

∑ (p|h|̂ q)Epq̂ + pq

1 2

E DFT = 2 ∑ (i|h|̂ i) + 2 ∑ (ii|jj) + Exc(ρ) i

(4)

where Ê is the usual singlet excitation operator. Once a complete active space partitioning of the molecular orbitals is introduced, the Hamiltonian may be normal-ordered with respect to the determinant describing the core electrons, |Φcore⟩. Ĥ = Ecore +

̂ |u)E ̂ + ∑ (t |hcore tu tu

1 2

ij

(8)

for a closed-shell singlet. Here, Exc(ρ) is the exchangecorrelation energy, which is interpreted generally to include both DFT and exact exchange contributions. If the indices i, j are allowed to include all electrons, as is usually done, then EDFT is an expression for the KS-DFT energy. Typically, this energy is obtained in a self-consistent manner as in HF; however, it can also be obtained as a nonself-consistent postprocessing step once a suitable electron density has been obtained. This nonself-consistent approach has been used for many years and originates in the work of Lie and Clementi.79,80 More recently, this approach has been suggested by Burke and co-workers as a remedy in cases where self-consistent iteration exacerbates the errors of the exchange-correlation functional.81−85 If the sums over i and j in eq 8 are restricted to include only the core electrons (as defined by the CAS partitioning of the molecular orbitals), then the energy becomes the nonself-consistent DFT energy of a cationic electron density. We replace the mean-field energy of the

∑ (pq|rs)(Epq̂ Erŝ − δqrEpŝ ) pqrs

(7)

Hartree−Fock (HF) theory can be obtained from this formalism if no electrons are active, and the molecular orbitals variationally minimize the energy given in eq 6. In that case, the ̂ |q) (eq 7), becomes the effective one-electron operator, (p|hcore usual Fock operator. Generally, CAS methods can be viewed as an operator-level embedding71 approach: an exact treatment of the active electrons embedded within a mean-field treatment of the core electrons (eq 5). The embedding potential in a CASCI method is simply the mean-field potential of the electrons. As a result there is an imbalance in the handling of the electron−electron interactions: electrons in the active orbitals feel the exact potential of the other active electrons but only the mean-field potential of the electrons in the core orbitals. The electrons in the core orbitals feel both a mean-field potential of the other electrons in core orbitals and the electrons in the active orbitals. This view of CASCI as a type of HF embedding that relies on the orthogonality of the molecular orbitals to partition the system has been previously recognized by Manby and coworkers.72 Embedding can be performed more robustly by using optimized effective potential (OEP) methods to improve on the mean-field potential.73−77 It is also possible to exploit orbital orthogonality to embed two DFT subsystems exactly.72 To obtain exact DFT-in-DFT embedding, the embedding potential necessarily includes exchange-correlation effects. Use of orbital orthogonality also allows wave function subsystems to be embedded in a DFT subsystem;78 however, these methods typically partition systems based on spatial groupings of atoms rather than by defining orbital subspaces. In the present work, we take inspiration from this perspective of CAS-type methods as a type of embedding and attempt to introduce an account of dynamic electron correlation into the CASCI energy. In Kohn−Sham (KS) DFT, the electronic energy is obtained as a function of the one-particle density, ρ, viz.

(2)

u

∑ [2(pq|ii) − (pi|qi)] i

†

∏ βû |Φcore⟩

(6)

ij

and (p|ĥcore|q) is the sum of the bare one-electron operator and the effective one electron potential of the core electrons, viz.

† βĵ |⟩

A complete N-electron Slater determinant (indexed I, J) can be generated by acting Nactive creation operators on the reference determinant. α Nact

∑ [2(ii|jj) − (ij|ij)]

∑ (tu|vw)(Etû Evŵ − δuvEtŵ ) tuvw

(5)

In the above equation, Ecore is the energy of the core electrons, viz. 1132


Article

Journal of Chemical Theory and Computation

CASCI energy are largely unchanged. The new terms that appear in the gradient expressions have a negligible impact on the computational cost of the gradient; when pure density functionals are applied, the cost even decreases slightly relative to the canonical CASCI gradient. Here, we will describe the modifications required for the gradient of the DFT-corrected CASCI energy. Analytic formulations of the FOMO-CASCI and CISNO-CASCI derivative are presented elsewhere.64,91 We begin by noting that the DFT-corrected CASCI energy can be written in terms of the DFT energy of the core electrons and the one- and two-particle density matrices of the active electrons (γ and Γ, respectively).

reference state, eq 1, in the CAS partitioned Hamiltonian with the DFT energy from eq 8. The resulting method constructs CASCI wave functions using a modified Hamiltonian that includes the effect of dynamic electron correlation between core electrons. Ĥ eff = EDFT +

̂ |u)E ̂ + ∑ (t |hcore tu tu

1 2

∑ (tu|vw)(Etû Evŵ − δuvEtŵ ) tuvw

(9)

This attempts to correct the energy of the core electrons for dynamic correlation, while leaving the interactions of active electrons unchanged. The mean-field interaction between core and active electrons also remains. When Exc is the exact (meaning HF) exchange functional with no correlation functional, eq 8 becomes equivalent to eq 6; and EDFT is simply the mean-field energy of the core electron density, and the usual electronic Hamiltonian is recovered. The molecular orbitals used in this approach could be obtained in any manner deemed appropriate for a particular application. In the present work, we will use CASSCF orbitals, floating occupation molecular orbitals (FOMO),86,87 and configuration interaction singles natural orbitals (CISNOs).88 While it is certainly possible to optimize the orbitals in a self-consistent manner with the CASCI wave function while including the DFT correction (as in CASSCF), we will not pursue that possibility here. It should be noted that in the two limits of active space size, the DFT-corrected CASCI method reduces to an established approach. In the limit where no electrons are active, the method is equivalent to nonself-consistent DFT. In the limit where all electrons and orbitals are active, the method is equivalent to full CI.89 The approach described above has several advantages and disadvantages as a result of its simplicity. For a given set of nuclear coordinates, the CASCI vertical excitation energies remain unchanged. This is an advantage when intersections between electronic states are important. CASCI properly describes the topology of conical intersections; this is not damaged by introducing the DFT correction. On the other hand, the CASCI vertical excitation energies are well-known to be systematically overestimated, and our DFT correction offers no improvement in this regard. By including only the portion of the dynamic electron correlation that is identical for every electronic state, treatments of Rydberg excitations, for example, will not be improved. Essentially, the CASCI wave function must provide a reasonable starting point for the description of the electronic states of interest for the DFT correction to be potentially useful. Since the DFT correction is added entirely within a single invariant orbital subspace, the orbital invariance properties of CASCI are also retained (although other formal invariance properties are lost).90 Similarly, by adding the DFT correction only to the core electrons, there is no double counting of correlation effects; of course, this is at the expense of neglecting all correlation involving active electrons that is not included in the CASCI wave function. Finally, the computational expense due to the DFT correction is negligible. This allows the DFT-corrected CASCI method to be applied to systems with hundreds of atoms that are far too large for more rigorous multireference methods. Obviously, the applicability of this method remains limited by the size of the active space; if one cannot afford to form the CASCI wave function, the method cannot be applied. A final advantage of the present approach is that the equations defining analytic derivatives of the DFT-corrected

ECASCI(DFT) = E DFT +

∑ γpqhpq + ∑ Γpqrs(pq|rs) pq

pqrs

(10)

Two different derivatives of EDFT appear in the gradient expressions. First, and most obviously, is the derivative of EDFT with respect to the nuclear coordinates. ξ E DFT = 2 ∑ (i|h|̂ i)ξ + 2 ∑ (ii|jj)ξ + Excξ (ρ) i

(11)

ij

Here, the superscript, ξ, implies that the derivative of these terms has been taken with respect to the nuclear coordinates, and only the so-called skeleton terms92 are included (changes in the wave function are included elsewhere). The other type of derivative that is needed is the derivative of EDFT with respect to rotations of molecular orbitals. These derivatives of the usual core energy can be constructed in terms of a Fock-like matrix. f pi = (p|h|̂ i) +

∑ [2(pi|jj) − (pj|ij)] (12)

j

When the DFT correction is included, this Fock-like matrix needs to be consistent with the definition of EDFT that has been used. f piDFT = (p|h|̂ i) + 2 ∑ (pi|jj) + [Vxc(ρ)]pi j

(13)

Note that only first derivatives of the exchange-correlation energy are needed: derivatives with respect to changes in the nuclear coordinates and derivatives with respect to changes in the core density. Finally, we note that analytic expressions for nonadiabatic derivative coupling vectors remain unchanged. This can be seen from the Hellmann−Feynman expression for the nonadiabatic coupling vector, viz. ⟨Ψ

∂ ⟨Ψ(m)| ∂ξ Ĥ |Ψ(n)⟩ d (n) | Ψ ⟩≈ dξ E(n) − E(m)

(m)

(14)

By the orthogonality of CASCI wave functions, derivatives of the core energy, Ecore, do not contribute to the coupling vectors; it follows that the DFT corrections we introduce also do not contribute to the coupling vectors. There have been other attempts to improve CI methods by introducing correction terms taken from DFT.79,80,93−97 Of particular relevance is a technique developed by McDouall that uses the DFT correlation energy to improve CASSCF (or, more generally, multiconfiguration self-consistent field (MCSCF)) energies.98−101 Although the intention is the same, McDouall’s formulation is slightly different than the one presented here. We directly compute the exchangecorrelation energy of the core electron density, whereas McDouall determines the DFT correlation energy of the 1133


Article


the functional itself). Since density functional methods are generally transferable between chemical systems, it stands to reason that the DFT-corrected CASCI approach will be more transferable than other semiempirical CASCI methods, which often require extensive tuning before they can be applied to new systems.

complete electron density and removes the DFT correlation energy of the active electron density. As a result, it is unclear how this approach should be extended to excited states. If the DFT correction is evaluated for each state individually, the topology of potential surfaces around intersections would be destroyed. If this approach was applied instead to a stateaveraged electron density, then the topology of conical intersections would be preserved. Our approach avoids this ambiguity by operating only on the core electron density. Correction through the exclusive use of correlation functionals appears as a subset of the possible approaches described in the present work. The fundamental challenge of developing methods that attempt to combine DFT with multiconfigurational wave functions is that it is not known how DFT affects the coupling between two arbitrary Slater determinants (an off-diagonal element of the Hamiltonian). One possible solution is to introduce a phenomenological definition for the Hamiltonian that includes contributions from DFT. This is done in Grimme’s MRCI/DFT method.96,97 The advantage of this approach is that the electronic states are obtained from diagonalization of a single Hamiltonian and, therefore, should exhibit the proper behavior in the vicinity of conical intersections. An alternative is to develop postprocessing corrections to multiconfigurational wave functions. In these methods, typically, the multiconfigurational wave function is constructed in the usual manner, and the energy of each state is evaluated using DFT. Recently, these methods have been formulated in terms of the on-top pair density. Gusarov et al. as well as Cremer and co-workers developed a CAS-DFT method based on postprocessing corrections to CASSCF wave functions in terms of both the total density and the on-top pair density.95,102,103 A similar idea has been recently pursued by Gagliardi and co-workers, who are beginning to use their multiconfigurational pair-density functional theory to treat excited states.104 The preliminary results from their method for vertical excitation energies are quite promising.105,106 Related approaches have also been pursued by Savin and coworkers.107,108 It is also possible to formulate these corrections in the context of spin projected Hartree−Fock.109 The inherent disadvantage of postprocessing-type DFT corrections is that they necessarily destroy the topology of intersections between electronic states. One may draw an analogy to state-specific CASPT2 (which suffers from the same problem) as compared to multistate CASPT2 (which properly describes conical intersections). In certain contexts, such as the prediction of vertical excitation energies, there is no issue applying statespecific methods. In the context of photochemical dynamics where conical intersections play a critical role, it is crucial that the topology near seams of intersection be described correctly. A conservative view of the DFT-corrected CASCI approach is that it is simply a semiempirical CASCI method. Indeed, the present approach shares some similarities with the usual semiempirical CASCI methods, since the CASCI wave function remains unchanged and the Hamiltonian is parametrized to account for the neglected dynamic electron correlation. Semiempirical CASCI methods have been successful in applications to photochemical dynamics.56−59,66−69 In contrast to the usual semiempirical methods where the one- and twoelectron integrals appearing in the Hamiltonian are parametrized, in the present work, the only adjustable parameter is the choice of functional with which to evaluate the energy of the core electrons (and any adjustable parameters contained in

3. COMPUTATIONAL DETAILS The DFT-corrected CASCI method has been implemented within a developmental version of the graphical processing unit (GPU) accelerated quantum chemistry package, TeraChem.110−112 All CASCI computations reported in the present work were performed with this program. Complete active space second-order perturbation theory (CASPT2)28 computations were performed with the Molpro package.113 Second-order approximate coupled cluster singles and doubles (CC2) computations were performed with the density-fitted implementation114 in Turbomole.115 Throughout the manuscript, CAS active spaces will be specified as m/n, where m is the number of active electrons, and n is the number of active orbitals. In principle, the molecular orbitals defining the core and active spaces can be obtained by any technique. Three different schemes were used to obtain the molecular orbitals: CASSCF, FOMO,87 and CISNO.88 In the case of FOMO, a temperature parameter, β, of 0.35 au and Gaussian broadening of orbital energy levels is used in this work. The CISNOs are always constructed by including the first singlet excited state in the one-particle density matrix. Of the schemes used to obtain molecular orbitals in this work, CASSCF (or state-averaged CASSCF) is the most robust in terms of flexibility in the definition of the active space. CASSCF orbitals are also the most computationally demanding to obtain and may become unstable in certain regions of the potential energy surface. The FOMO-CASCI method uses fractional occupation number Hartree−Fock (FON-HF) orbitals, which are equivalent in cost to a usual Hartree− Fock computation and exhibit enhanced stability near regions of intersection between electronic states. In our implementation of FOMO-CASCI, we allow only the orbitals that will be active in the CASCI wave function to be fractionally occupied in the FON-HF procedure. This allows all of the active orbitals to enter the FON-HF energy expression and often provides a good set of molecular orbitals with which to construct the CASCI wave function. Previously, we have shown that the cost of a FOMO-CASCI gradient computation is similar to the cost of a TDDFT gradient (for a small active space with less than six active orbitals).64 We have found FOMO-CASCI to be a practical and efficient alternative to state-averaged CASSCF for simulations of photochemical processes. CISNO-CASCI is another alternative to CASSCF that works well when CIS provides a reasonable description of the electronic states of interest. The natural orbitals of the CIS wave function are used to define the active space for the CASCI wave function. Since these orbitals are generated from excited-state wave function(s), they can be quite useful for describing excited-state processes. This method is roughly intermediate in computational expense between that of FOMO-CASCI and CASSCF. Additional details of the FOMO-CASCI and CISNO-CASCI approaches can be found in refs 87 and 88, respectively. To study the proton transfer dynamics of 10-hydroxybenzo[h]quinoline (HBQ) and 2-(2′-hydroxyphenyl)benzothiazole (HBT), classical molecular dynamics simulations are performed on the first singlet excited states computed with the CASCI and 1134


Article


these ESIPT processes. These simulations demonstrate that the DFT-correction is required to obtain a qualitatively correct description of these reactions from a CASCI wave function. We also obtain the isotope effect on the proton transfer time, demonstrating that our approach provides a physically reasonable description of the ESIPT reactions. In these dynamics simulations, we exhaustively sample the initial conditions to show the efficiency of the DFT-corrected CASCI approach. Finally, we apply our approach in a QM/ MM simulation of an excited-state proton transfer reaction that occurs in aqueous solution. The QM region in this simulation contains more than 150 atoms and is offered as a demonstration of the scalability of the DFT-corrected CASCI approach. A simple system that illustrates some of the problems with single reference methods is the torsion of ethylene. In its planar geometry, ethylene is well described by a single Slater determinant, but when it is twisted 90°, the π and π* orbitals become degenerate and the wave function gains biradical character. This system is shown in Figure 1. Here, CASPT2

DFT-corrected CASCI methods. CISNOs as described above were used in the molecular dynamics simulations. Initial conditions are sampled from a Wigner distribution constructed from harmonic vibrational frequencies obtained with HF (for CASCI) or DFT (for DFT-corrected CASCI methods, using the same functional). Trajectories are integrated for 200 fs in 0.5 fs time steps. There were issues with energy conservation in a small fraction of the trajectories. These were attributed to crossings between the S1 and S2 states computed in the CIS method; this led to discontinuities in the CISNOs within the CASCI active space. Trajectories where such discontinuities exist were discarded and not included in the final averaging. Although quite common throughout the literature, the removal of trajectories in this manner will systematically bias the results away from a particular outcome (i.e., the outcome associated with the discontinuity) and should be avoided whenever possible. In the present case, the discontinuities occur well after the proton transfer event, and the average proton transfer time is unaffected by the inclusion or exclusion of the problematic trajectories. Note that the issue arises from the choice of molecular orbitals, not the DFT correction (which is the focus of the present work). If these simulations were to be extended to longer times, a different choice for the molecular orbitals would be required to avoid this issue. Each DFT-corrected CASCI simulation of excited-state proton transfer was averaged over 500 trajectories, and the canonical CASCI simulations were averaged over 100 trajectories. We also study excited-state proton transfer in aqueous 1Hpyrrolo[3,2-h]quinoline (HPQ) via classical excited-state molecular dynamics simulations. Here, the effect of solvent is modeled with a quantum mechanics/molecular mechanics (QM/MM) approach. The HPQ molecule and 50 water molecules are treated with the DFT-corrected CASCI method and FOMO molecular orbitals. A 25 Å sphere of TIP3P water surrounds the HPQ molecule. A soft spherical harmonic restraint (with a force constant of 10 kcal mol−1 Å−2 and radius of 25 Å) is applied to the system in order to prevent evaporation. The initial conditions for the molecular dynamics trajectories were taken from a classical molecular dynamics simulation using the general Amber force field (GAFF)116,117 for HPQ and TIP3P118 for the water molecules. This simulation was equilibrated in an NVT ensemble at 300 K for 500 ps. After equilibration, the simulation was continued for 1 ns; the initial positions and momenta come from snapshots taken every 100 ps from this trajectory. The subsequent excited-state molecular dynamics trajectories are integrated for 600 fs using 0.5 fs time steps in an NVE ensemble.

Figure 1. Rigid torsion of ethylene computed at four levels of theory B3LYP/cc-pVDZ, FOMO-CASCI/cc-pVDZ, CASPT2/cc-pVDZ, and FOMO-CASCI/cc-pVDZ with a B3LYP correction for core correlation. The CASCI methods use a 2/2 active space; the CASPT method uses a 12/12 active space. The geometries with a twist of 0 degrees are obtained via optimization at each level of theory.

(using a 12/12 active space) should provide an accurate description of the torsional potential. B3LYP, on the other hand, cannot capture an exact degeneracy between occupied and unoccupied orbitals; it also cannot describe an open-shell singlet electronic state. As a result, the torsional potential exhibits a cusp. A common solution to this problem is to break spin symmetry and use unrestricted molecular orbitals. While this removes the cusp from the torsional potential, it also introduces spurious mixing between singlet and triplet states; this is unacceptable for simulations of photochemical processes. FOMO-CASCI can remove this cusp while preserving spin symmetry since the underlying FON-HF allows the π and π* orbitals to have equal fractional occupation; the CASCI expansion captures the biradical character of the wave function. However, despite producing a smooth potential, the barrier is somewhat overestimated relative to CASPT2. When a B3LYP correction is added to FOMO-CASCI, a smooth potential is obtained, the spin symmetry is preserved, and the barrier is

4. RESULTS In this section, we present a series of computations chosen to illustrate the utility of our DFT-corrected CASCI approach. First, we examine the torsion of ethylene to demonstrate that biradicaloid electronic structure can be described with this method. We test the performance of the DFT-corrected CASCI approach for bond dissociation where full configuration interaction (FCI) reference data is available. Then, we show the topology of the ground and first excited state in the vicinity of a conical intersection to show that the seam space exhibits the proper dimensionality. Next, we apply our approach to excited-state intramolecular proton transfer (ESIPT) reactions. Here, we demonstrate that the DFT-correction is capable of greatly improving reaction barriers as compared to the standard CASCI treatment. Next, we perform dynamical simulations of 1135


Article


Figure 2. Rigid torsion of ethylene computed with FOMO-CASCI/cc-pVDZ using a variety of active spaces and a B3LYP correction for core correlation. The CASPT method uses a 12/12 active space. The upper panels show the energy relative to the planar geometry, and the lower panels show the error in the relative energy computed with the CASCI methods with respect to CASPT2. The left panels show symmetric active spaces, and the right panels show active spaces including one π* orbital. Note that a few points have been omitted due to discontinuities in the definition of the FOMO-CASCI molecular orbitals.

active spaces is likely a result of fortuitous error cancellation. Note that, when a symmetric active space is applied, the best results are obtained from the largest active spaces. With a 12/12 active space, the effect of the DFT correction on the torsional potential is minimal. Using the nonsymmetric active spaces, the largest active space does not produce the most desirable results. The 8/5 and 8/8 active spaces exhibit the worst performance; this is likely a result of errors due to the mean-field embedding potential. The DFT-corrected CASCI torsional potentials are well-behaved with all of the active spaces, although some perform better than others. We do not regard the active space dependence of the results as a deficiency of the method. Methods that utilize active spaces are generally not “black box” approaches; the DFT-corrected CASCI method we introduce here is no different. The application of this method to a chemical problem will require testing and calibration to ensure a satisfactory active space has been chosen. To test the DFT-corrected CASCI for bond breaking reactions, we compute bond dissociation potential energy curves for the HF and CH4 molecules (see Figure 3). In this case, CASSCF orbitals and a 2/2 active space are used. The DFT correction to the CASSCF energy is computed with the B3LYP functional. Note that the CASSCF orbitals are not optimized in the presence of the DFT correction; we will consider self-consistent optimization in future work. Benchmark potential energy curves computed with FCI from ref 119

greatly improved relative to canonical B3LYP, although slightly underestimated relative to CASPT2. Our assumption is that this form of the DFT correction is best suited for application with compact CASCI wave functions where minimal active spaces are employed. Because we only treat the core electrons with DFT, as the active space is enlarged, one begins to trade the DFT treatment of the electron correlation for the CASCI treatment. It also should be expected that our use of a mean-field embedding potential will introduce larger errors as the number of active electrons increases (assuming the total number of electrons is much greater than the number of active electrons). In Figure 2, we examine the torsional potential of ethylene using FOMOCASCI wave functions with DFT corrections from B3LYP. We consider active spaces that are expanded symmetrically: 2/2, 4/ 4, 6/6, etc. as well as active spaces that are expanded only in the occupied space: 2/2, 4/3, 6/4, etc. In the former case, there is a trade off between the dynamic electron correlation captured in the active space and the correlation captured by the DFT correction. With the nonsymmetric active spaces, the total amount of dynamic electron correlation decreases as the active space is expanded. We find that good results are obtained with the 4/3 and 4/4 active spaces. The small, 2/2 active space underestimates the torsional barrier when the DFT correction is applied; larger active spaces overestimate the torsional barrier. Consequently, the performance of the 4/3 and 4/4 1136


Article


Figure 3. Bond dissociation of HF (left panels) and one bond of CH4 (right panels) computed with FCI, CASSCF, CAS(B3LYP)-SCF, and B3LYP. CASSCF computations use a 2/2 active space consisting of one σ orbital and one σ* orbital. Potential energy curves are shown in the upper panels and are reported relative to the minimum energy on the curve. Errors relative to FCI are shown in the lower panels. Computations on HF use a 631G** basis, and computations on CH4 use a 6-31G* basis. FCI and B3LYP energies are taken from ref 119.

are used as a reference. We find that the DFT-correction offers only a modest improvement over CASSCF in the case of HF and very little improvement in the case of CH4. For HF, the largest change is an improvement in the reaction energy seen at the dissociation limit (of about 5 kcal mol−1). This is likely a result of the fact that the largest changes in the electron density are captured by the CASCI wave function and the core electron density (where the DFT correction is applied) remains relatively constant throughout the bond breaking reaction. This implies that the DFT correction will be more useful for treating reactions involving larger molecules where changes in the core electron density are more pronounced. Importantly, CASSCF with our DFT correction easily outperforms canonical B3LYP. Our approach inherits the qualitative description of the bond breaking reaction from the CASCI wave function. This allows for a reasonable description of bond dissociation while preserving spin symmetry. We offer this test to demonstrate that the DFT correction as formulated here is well-behaved and does not spoil the description of bond dissociation offered by the CASCI wave function. A theoretical analysis of the DFT-corrected CASCI method indicates that it should be able to describe the topology of the GH branching plane in the region around a minimum energy conical intersection (MECI).16 We test this hypothesis computationally by optimizing an MECI of the anionic chromophore of the green fluorescent protein. The GH branching plane in the area of that MECI is shown in Figure

4. Indeed, the FOMO-CASCI method with a B3LYP correction for core correlation is capable of describing the potential energy surfaces in the vicinity of the conical intersection. Both the MECI geometries and the potential energy surfaces are in good agreement with previous work on this compound using FOMO-CASCI and SA-CASSCF. It should be noted that DFT-corrected CASCI should not be expected to provide significant improvement upon the usual CASCI methods for the prediction of MECI geometries. Changes in the vertical excitation energies typically introduce larger changes in the MECI geometry than would the introduction of core electron correlation. Generally speaking, DFT-corrected CASCI should be expected to provide a description of MECI geometries that is similar to the underlying CASCI methodology. In the cases of the torsion of ethylene and the prediction of MECI geometries, CASCI performs well without corrections from DFT. The next question is whether or not the DFTcorrection offers improvement on the underlying CASCI. Here, we will examine the proton transfer reaction coordinate on the S1 state of HBQ120−125 and HBT.126−140 These are two model systems for ultrafast ESIPT reactions.141,142 The proton transfer reaction in HBQ is described as an active transfer mode, where, after excitation, the proton transfer precedes other structural changes.122,142,143 Recent estimates of the proton transfer time of HBQ are 13 ± 5 fs,141 25 ± 15 fs,121 and between 24 and 45 fs (depending on the choice of probe wavelength).122 Proton transfer times in HBT are slightly longer, although still ultrafast. 1137


Article


Figure 4. Topology of the GH branching plane around the minimum energy conical intersection of the anionic chromophore of the green fluorescent protein computed at the FOMO-CAS(2/2)-CI/6-31G** level of theory with a B3LYP correction for core correlation is shown ∂ here. Here, the G vector (shown in red) is g ⃗ = ∂R⃗ (ES0 − ES1), and the ∂

H vector (shown in blue) is h ⃗ = (ES1 − ES0)⟨ΨS0 ⃗ ΨS1⟩. Further, the ∂R H vector is Gram-Schmidt orthogonalized against the G vector.

Proton transfer in HBT is described as passive; that is, bond rearrangement precedes proton transfer and acts as a gating mechanism for the reaction.142,143 Experiments on HBT put the proton transfer time in the 50 fs143 to 62 fs142 range. Theoretical studies of these molecules also predict ultrafast proton transfer and a barrierless or nearly barrierless excitedstate reaction coordinate. The extensive experimental and theoretical investigations of these molecules make them ideal test cases for the DFT-corrected CASCI method. To begin, we examine the excited-state potential surface along the proton transfer coordinate of HBQ and HBT. The proton transfer coordinate is defined as the difference between the OH distance and the NH distance, i.e. RPT = R OH − RNH

Figure 5. Minimum energy proton transfer path of HBQ (upper panel) and HBT (lower panel) on the S1 state computed at the CC2/ cc-pVDZ, CISNO-CAS(2/2)-CI/6-31G*, and PBE-corrected CISNOCAS(2/2)-CI/6-31G* levels of theory.

been established, previously, that the S1 state of both HBQ and HBT retains its ππ* character throughout the proton transfer process.122,137,144,145 While limiting the active space to include π orbitals is certainly insufficient for ESIPT in general,146 it is a reasonable approximation for HBQ and HBT (and, as a first approximation, for molecules containing extended π systems). For both HBQ and HBT, CC2 predicts barrierless proton transfer. CASCI provides a reasonable description of the ESIPT reaction coordinate for HBQ although a small barrier (1.3 kcal mol−1) is present. In contrast, CASCI is extremely poor for the ESIPT potential of HBT; a large (9 kcal mol−1) barrier to proton transfer is present. The addition of a DFT correction to CASCI improves the excited-state potential substantially. Here, we use the PBE exchange-correlation functional to describe the correlation of the core electrons. In the case of HBQ, the DFTcorrected CASCI is in remarkably good agreement with CC2. For HBT, the DFT-corrected CASCI removes the barrier that is present on the CASCI potential; however, it slightly underestimates the magnitude of the reaction energy relative to CC2. For both prototypical ESIPT systems, the DFTcorrection improves CASCI and is in good agreement with CC2. CC2 scales with the fifth power of system size, 6(N 5),147,148 whereas the DFT-corrected CASCI method scales with the fourth power of system size, 6(N 4). In practice, our implementation of DFT-corrected CASCI exhibits quadratic, 6(N 2), scaling due to the exploitation of sparsity

(15)

The proton transfer coordinate is constructed by tracing the minimum energy path from the enol Franck−Condon point to the keto minimum on the S1 potential surface. These potential curves are reported for HBQ and HBT in Figure 5. To facilitate comparison between different methods, we introduce a relative proton transfer coordinate RPT − RPT,enol RPT,rel = RPT,keto − RPT,enol (16) where the enol geometry is constrained to have the same OH bond length as the S0 enol minimum and all other coordinates relaxed. The relaxed geometries are consistent with the level of theory used to compute the energetics. We use CC2 to provide a reference potential for HBQ and HBT; CC2 is expected to perform well for these systems and has previously been applied successfully to these and similar molecules.144 The CASCI computations are performed with CISNO-CAS(2/2)-CI/631G* and include one π and one π* in the active space; it has 1138


Article

Journal of Chemical Theory and Computation in the atomic orbital basis.60,64 When the great disparity in computational expense between CC2 and DFT-corrected CASCI is considered, the level of agreement between these two methods is quite remarkable. As a rough point of comparison, for HBT, the CC2/cc-pVDZ gradient evaluation is about 32× more expensive than the corresponding calculation with DFT-corrected CISNO-CAS(2/2)-CI/ccpVDZ and about 55× more expensive than DFT-corrected FOMO-CAS(2/2)-CI/cc-pVDZ.149 Note that the discrepancy in computational cost will grow cubically with system size and that HBQ and HBT are still small molecules (24 and 25 atoms, respectively). It is interesting to consider the functional dependence of these reaction paths using DFT-corrected CASCI. As formulated, this approach is compatible to all usual density functionals; it is possible to apply pure density functionals, hybrid density functionals, and range-separated hybrid functionals. In Figure 6, we have computed the excited-state proton

functionals, PBE0 and B3LYP, predict small barriers to the proton transfer reaction of between 1 and 2 kcal mol−1. Even these small barriers are sufficient to slow the proton transfer time significantly and are not consistent with the experimental observations. The best agreement with CC2 is obtained with the pure PBE functional. With this DFT correction, no energetic barrier to proton transfer appears, and the reaction energy is improved (although underestimated by about 2 kcal mol−1). The range-separated hybrid, ωPBE, again significantly overestimates the reaction energy (by about 6 kcal mol−1); however, it predicts a barrierless reaction path. In that way, it agrees better with the experimental observations than either PBE0 or B3LYP. Perhaps counterintuitively, the behavior of ωPBE may be the most desirable in general. The overestimated reaction energies are likely inherited from the underlying CASCI wave function. However, it is still able to remove the spurious reaction barriers from the potential surface. One might view ωPBE as introducing the smallest perturbation to the CASCI results; this may lead to greater stability of the computations, since they are performed nonself-consistently. We plan to address this question in more detail in future work. For the remainder of the present work, we will use the DFT correction computed with PBE to study the excited-state proton transfer dynamics of HBQ and HBT. The origin of spurious reaction barriers to ESIPT has been the subject of discussion in previous theoretical studies of similar systems.32,150 As a result of the minimal 2/2 active space we employ, there is essentially no dynamic electron correlation included in the CASCI wave function. In the case of the DFTcorrected CASCI, it follows that the only account of dynamic correlation is through the DFT treatment of the core electrons. Since the barriers to ESIPT are removed through the introduction of the DFT correction, we conclude, a posteriori, that these spurious barriers result from the neglect of core electron correlation. The next question is why these errors were exaggerated in HBT as compared to HBQ. This seems to be related to the passive nature of the proton transfer in HBT, the bond rearrangement that precedes proton transfer in HBT. The lack of core electron correlation creates a barrier to this rearrangement. In the case of HBQ (active ESIPT), the primary driving force of ESIPT is the electronic reorganization that occurs upon excitation. Contributions from bond rearrangement are secondary, and thus only small errors result from the lack of core correlation (the gross electronic reorganization is captured within the minimal CAS active space). The size of the active space also affects the description of the potential surface. One might reasonably expect that inclusion of σ and σ* orbitals in the active space would significantly improve the treatment of the reaction path and possibly remove the spurious barrier. To test this hypothesis, we have performed CASSCF computations using the minimal 2/2 active space and a larger 6/5 active space that additionally includes two σ orbitals and a σ* orbital. The results of these computations for HBQ and HBT are contained in Figure 7. The active molecular orbitals from these computations are shown in Figures 8 and 9. We find very little difference between the 2/2 and 6/5 active spaces. For HBQ, the two potential energy curves are nearly identical, although there is a very slight (less than 1 kcal mol−1) improvement in the reaction energy. For HBT, the results are again very similar for the two active spaces. The inclusion of the σ and σ* orbitals do not significantly affect the spurious reaction barrier, and they do not significantly improve the reaction energy. For comparison, results from CISNO-CASCI

Figure 6. Minimum energy proton transfer path of HBQ (upper panel) and HBT (lower panel) on the S1 state computed at the CC2/ cc-pVDZ and CISNO-CAS(2/2)-CI/6-31G* levels of theory. CASCI computations include corrections from PBE, PBE0, B3LYP, and ωPBE.

transfer reaction path using CISNO-CAS(2/2)-CI/6-31G* with DFT corrections from PBE, PBE0, B3LYP, and ωPBE. For HBQ, the potential energy curves obtained with PBE, PBE0, and B3LYP are very similar and are all in good agreement with our CC2 benchmark. The range-separated hybrid, ωPBE, overestimates the reaction energy by about 5 kcal mol−1 relative to CC2. For HBT, however, both hybrid 1139


Article


Figure 8. Molecular orbitals of HBQ defining a 6/5 active space obtained from a CAS(6/5)-SCF/6-31G* computation (canonical CASSCF orbitals). (a) σ orbital corresponding to the OH bond, (b) σ orbital corresponding to the NH bond, (c) π orbital, (d) π* orbital, and (e) σ* orbital. Figure 7. Minimum energy proton transfer path of HBQ (upper panel) and HBT (lower panel) on the S1 state computed at the CC2/ cc-pVDZ, CAS(2/2)-SCF/6-31G*, CAS(6/5)-SCF/6-31G*, CISNOCAS(2/2)-CI/6-31G*, and FOMO-CAS(2/2)-CI/6-31G* levels of theory. The 2/2 active spaces include one π orbital and one π* orbital; the 6/5 active spaces additionally include two σ orbitals and one σ* orbital (see Figures 8 and 9).

and FOMO-CASCI with the 2/2 active space are also included. The CISNO-CASCI potential energy curves are in the worst agreement with the CC2 reference (which speaks to great improvement offered by the DFT correction). For HBQ, the FOMO-CASCI potential energy curve is in the best agreement (of the methods that neglect dynamic electron correlation); for HBT, however, FOMO-CASCI exhibits a discontinuity near the Franck−Condon point and, therefore, cannot be applied in this case. Increasing the size of the active space is not a tractable approach to improve the description of the excited-state proton transfer reactions by CASCI wave functions.150 In passing, it is also worth noting that inclusion of the σ and σ* orbitals in the active space during a molecular dynamics simulation would be technically difficult. Next, we will explore the effect of the DFT correction on the dynamic evolution of the model ESIPT systems. The results of this analysis are shown in Figure 10 for HBQ and HBT. In these figures, we show the time evolution of the average OH and NH distances involved in the ESIPT process. It is common to define the proton transfer time as the point where the OH and NH distances are equal; we will use this definition in the present work. However, it should be noted that this is not the experimental observable used to determine proton transfer times. When comparing with experiments that determine rate

Figure 9. Molecular orbitals of HBT defining a 6/5 active space obtained from a CAS(6/5)-SCF/6-31G* computation (canonical CASSCF orbitals). (a) σ orbital corresponding to the OH bond, (b) σ orbital corresponding to the NH bond, (c) π orbital, (d) π* orbital, and (e) σ* orbital.

constants from the spectroscopic signature of the keto form of HBT or HBQ, the definition of proton transfer described above will likely underestimate the experimental proton transfer time. For both molecules, we compare classical molecular dynamics simulations on the S1 electronic state with and without the DFT correction for the dynamic correlation of the core electrons. In the case of HBT, the proton transfer time is about 9 fs with DFT correction as opposed to 70 fs without DFT correction. The small barrier that is present on the HBQ 1140


Article


Figure 10. Time-evolution of the oxygen−hydrogen and nitrogen− hydrogen interatomic distances of HBQ (upper panel) and HBT (lower panel) during the first 200 fs following photoexcitation to the S1 electronic state computed at the CAS(2/2)-CI/6-31G* and the PBE-corrected CAS(2/2)-CI/6-31G* levels of theory.

Figure 11. Time-evolution of the oxygen−hydrogen/deuterium and nitrogen−hydrogen/deuterium interatomic distances of HBQ (upper panel) and HBT (lower panel) during the first 200 fs following photoexcitation to the S1 electronic state. All computations are performed at the PBE-corrected CAS(2/2)-CI/6-31G* level of theory.

reaction coordinate (Figure 5) slows the progress of the reaction by about a factor of 8. Only the proton transfer time obtained with DFT-corrected CASCI is within the experimental uncertainty for this process.141 In view of the agreement between the DFT-corrected CASCI and CC2, we believe that 9 fs is a reliable theoretical estimate of the proton transfer time. For HBT, the estimated proton transfer time is 36 fs with DFT correction; in contrast, proton transfer is not observed within the first 200 fs on the canonical CASCI potential surface. This is attributable to the 9 kcal mol−1 barrier to proton transfer that is present. Clearly, without DFT correction CASCI cannot be successfully applied to the photochemistry of HBQ and HBT. Next, we determine the effect of isotopic substitution on the proton transfer times in HBQ and HBT using the DFTcorrected CASCI potential surface. We substitute the hydroxyl hydrogen with a deuterium; no other isotopic substitutions are made. For an active ESIPT process, as is the case for HBQ, the isotope effect is predicted to be √2 due to the fact that the proton transfer time is limited by the velocity of the proton. In the case of HBT, a passive ESIPT process, proton transfer is limited by bond rearrangement, and little to no isotope effect is predicted. Indeed these predictions are supported by experiment.142 In Figure 11, proton transfer times of 8.7 and 12.1 fs are observed for HBQ and deuterated HBQ, respectively. This is in good agreement with the expected isotope effect of √2. For HBT (Figure 11), proton transfer times of 32.1 and 40.3 fs

are observed. This isotope effect (of approximately 1.25) is more pronounced than what is observed in experiment.142 However, we note that our results do agree to within the experimental uncertainty of the measured proton (and deuteron) transfer time for HBT, ±10 fs.141 The description of ESIPT in HBT and HBQ by DFTcorrected CASCI is quite promising for the application of this method to other ESIPT processes. The efficiency of the DFTcorrected CASCI method makes it feasible to study these processes with methods for quantum nuclear dynamics that are usually bottlenecked by the computational expense of the electronic structure methods.3,6,30 This is especially interesting in the case of HBT, which has a gas-phase excited-state lifetime of 2.6 ps. Deactivation of the S1 state occurs via twisting about the central carbon−carbon bond. Direct simulation of the nonadiabatic dynamics of HBT has proven elusive, since an electronic structure method must be able to satisfy the following three criteria simultaneously for successful application to this problem. First, the method must be able to accurately describe the ESIPT process (this rules out conventional CAStype methods, see Figures 10 and 5). Next, the method must be able to describe the biradicaloid electronic structure of HBT near a seam of conical intersections between S0 and S1 (this rules out single reference methods such as TDDFT or CC2). Finally, the method must be efficient enough to simulate many 1141


Article


separated ωPBE exchange-correlation functional to evaluate the DFT correction (ω = 0.3 au). The active orbitals are localized on the HPQ molecule (Figure 12). The active space is comprised of one π orbital and one π* orbital. It has been previously established that this process proceeds along the S1 electronic state; this state retains its ππ* character throughout the process.157 However, the charge-transfer character of this state does increase after proton transfer occurs (Figure 12). In that way, one might think about the tautomerization of HPQ as a type of photoinduced proton-coupled electron transfer. We followed 10 independent trajectories of aqueous HPQ for 600 fs (additional computational details were given previously). While this is not a sufficient number of trajectories to obtain converged statistics (as was done for HBT and HBQ), it is sufficient to gain a qualitative understanding of the ESPT process. The purpose of these simulations is 2-fold: to explore the effect of bulk solvent on ESPT in aqueous HPQ and to demonstrate that the DFT-corrected CASCI method can be applied to larger molecular systems. Of the 10 trajectories, proton transfer reactions were observed in 6 of them, and, in 5 trajectories, the tautomerization process was completed within 600 fs. The time-evolution of the six reactive trajectories is shown in Figure 13. To simplify the analysis of

trajectories for several picoseconds (this rules out multireference methods such as CASPT2). The DFT-corrected CASCI, however, satisfies all three of these requirements; work is currently underway to directly simulate the nonadiabatic dynamics of HBT with this method. Finally, we examine the solvent-assisted tautomerization of aqueous 1H-pyrrolo[3,2-h]quinoline (HPQ) via excited-state proton transfer (ESPT). The ESPT process of HPQ has been studied both experimentally151−154 and theoretically.155−157 Experiments have been performed on gas-phase clusters (HPQ and a few solvent molecules) and HPQ in bulk solvent, while theoretical studies have focused on the small gas-phase clusters. Here, we will apply the DFT-corrected CASCI method to ESPT in aqueous HPQ. In order to simulate bulk solvation, we construct a 25 Å sphere of water around a single HPQ molecule (Figure 12). The 50 water molecules nearest the nitrogen

Figure 13. Time evolution of the HPQ reaction coordinate (as defined in eq 17) is shown for the six reactive molecular dynamics trajectories for the first 600 fs following excitation to S1.

the molecular dynamics trajectories, we define a reaction coordinate where we use the error function as a switching function to distinguish between the steps of the proton transfer process.

Figure 12. Representative geometries from the simulation of aqueous HPQ: (a) the QM region including HPQ and 50 water molecules and (b) the total QM/MM system including a sphere of water molecules with a radius of 50 Å. The active orbitals of HPQ and HPQ tautomer: (c) the π* orbital of HPQ, (d) the π* orbital of the HPQ tautomer, (e) the π orbital of HPQ, and (f) the π orbital of the HPQ tautomer.

RPT =

erf(α(R OH2 − R NH2)) − erf(α(R OH1 − R NH1)) 2 (17)

This reaction coordinate takes a value of −1 when the normal structure of HPQ is present. After the first proton transfer event, the reaction coordinate takes a value of 0. Finally, the tautomer of HPQ corresponds to a reaction coordinate of 1. The interatomic distances are defined in Figure 13, and units of Ångströms are used. The constant, α, is chosen arbitrarily to define the sharpness of the switching function; here, we take α = 4. This simple three-state model is sufficient to describe the reactions observed in our simulations: only tautomerization involving a single water molecule (double proton transfer) is observed. The same reaction mechanism is observed in each case (Figure 13). The first step in the tautomerization is a proton transfer from water to HPQ generating OH−. The

atoms of HPQ are treated quantum mechanically with the DFT-corrected CASCI method; the remaining water molecules are described with the TIP3P force field. We note that the QM region of these computations includes 171 atoms and 1161 basis functions. The energy and gradient evaluation at each time step took, on average, 170 s on a six-core Intel Core i75820K CPU clocked at 3.30 GHz and four Nvidia GeForce GTX 970 GPUs. On this workstation, one 600 fs trajectory takes approximately 2.5 days to complete. The electronic structure of the HPQ−(H2O)50 complex is treated at the FOMO-CAS(2/2)-CI/6-31G* level of theory using the range1142


Article


method depends entirely on the particular questions one wishes to address. Our approach is not well well suited to applications that require accurate predictions of vertical excitation energies. If accurate energetics on a single (ground- or excited-state) potential energy surface are of primary importance, then the DFT correction can offer significant improvement. As a result of these considerations, perhaps the most obvious application of this method is in the generation of potential energy surfaces for dynamical simulations of photochemical processes. Work is currently underway to interface DFT-corrected CASCI with the ab inito multiple spawning (AIMS) approach to nuclear quantum dynamics.30 Preliminary applications of the DFTcorrected CASCI approach in nonadiabatic molecular dynamics simulations of photochemical processes have been very encouraging. At present, we view the DFT-corrected CASCI approach as an alternative to semiempirical configuration interaction (CI) methods. It should find applications in simulations of photochemical processes where standard CAS-type methods are not suitably accurate, single-reference methods are not applicable, and multireference methods (MRCI, CASPT2, etc.) are too computationally demanding. The principal advantage of DFT-corrected CASCI compared to semiempirical CI is that there is no need for reparametrization of the Hamiltonian. However, both approaches do require validation against higher levels of theory. As a result, we certainly do not view the DFTcorrected CASCI approach as a replacement for existing excited-state methods. Rather, we view it as complementary to the more accurate but more expensive methods. This is especially true in the context of dynamical simulations where high computational efficiency is particularly important. The quadratic-scaling GPU-based implementation of the DFTcorrected CASCI approach makes dynamical simulations of systems with less than 50 atoms routine with even modest computational resources; dynamical simulations of larger systems (up to 500 atoms) are possible but do require significant computational effort. Much work is still needed to properly test and calibrate the DFT-corrected CASCI approach. At present, we treat the core electrons with popular density functionals (PBE, B3LYP, and ωPBE). In general, density functionals are not optimized to treat cationic densities, nor are they designed for nonselfconsistent application. We are optimistic that the performance of the DFT-corrected CASCI approach can be improved through the careful selection and optimization of a density functional. Another potential improvement is in the choice of embedding potential. Here, we use the mean-field potential, which neglects dynamic correlation between the core and active subsystems. It is not immediately obvious how to improve the embedding potential while still retaining the correct topology of intersections between electronic states, but this is a promising direction for future research. The DFT-corrected CASCI approach presented here is quite promising as a highly efficient method for treating multiple electronic states, and more work is warranted to further refine this method.

second step is proton transfer from the protonated HPQ to the OH−. The second step completes the cycle, regenerating a water molecule and producing the tautomer of HPQ. This mechanism is consistent with dynamics simulations of 7azaindole, a related molecule, performed in the condensed phase using a smaller number of QM water molecules (one or two), but a polarizable description of the MM waters.158 The most important observation from these simulations is that tautomerization typically occurs in a sequential (asynchronous) manner. For a gas-phase cluster, the tautomerization process completes in under 10 fs (from the first proton transfer event to the generation of the tautomer).157 In bulk, the time between proton transfer events varies greatly. In one case (Figure 13, blue curve), the tautomerization completes in a single, concerted step. In four other simulations, the process takes between 50 and 300 fs to complete. In the final reactive trajectory, the initial proton transfer event occurs within 100 fs of excitation but does not complete within 600 fs. The difference between gas-phase and solution-phase behavior can be partially explained by the fact that the intermediate geometries (the protonated HPQ and OH− or a deprotonated HPQ and H3O+) are unstable in the gas phase. In a gas-phase complex of HPQ and one or two solvent molecules, the transient species formed after the first proton transfer event cannot be stabilized; this forces the reaction to proceed quickly. In contrast, in aqueous solvent, the hydrogen bonding network of the surrounding water molecules can stabilize these intermediates and extend their lifetimes by (at least) an order of magnitude. Preorganization of the solvent also plays a role in determining the time scale of the process. In the gas phase, the solvent molecules are in the optimal location to facilitate concerted tautomerization. In bulk, it is possible not to have any solvent molecules in a position to undergo ESPT (as evidenced by the four nonreactive trajectories). It is also possible for the solvent to be in a suitable position for the first proton transfer event but not to be in a position to rapidly complete the tautomerization process. In the gas phase, the solvent molecules that mediate tautomerization are aptly described as a proton-conducting “water wire”. The molecules that make up the wire are optimally placed, and the energetics of the system prevent the proton from becoming trapped before tautomerization is completed. In contrast, we do not observe the formation of water wires in our simulations of aqueous HPQ. The hydrogen bonds between water molecules are sufficiently strong to both stabilize the transient species and to disrupt the interactions between water molecules and HPQ. These interactions between solvent molecules lead to a stepwise tautomerization mechanism in the condensed phase that can take hundreds of femtoseconds to complete.

5. CONCLUSIONS We have presented an efficient density-based correction to CASCI energies. By using an effective Hamiltonian that incorporates the effect of the dynamic electron correlation within the core orbitals, spurious reaction barriers that appear on CASCI potential surfaces can be removed. This formulation preserves the ability of CASCI to properly describe the topology of conical intersections. Coupled with the recent algorithmic advances in CASCI methods, the DFT-corrected CASCI approach can be extended to systems with hundreds of atoms or molecular dynamics simulations with extensive sampling of initial conditions. Whether or not this simple correction offers improvement over the canonical CASCI

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00893. Energies and optimized geometries of the GFP chromophore, HBQ, and HBT (PDF) 1143


Article


■

(30) Ben-Nun, M.; Quenneville, J.; Martínez, T. J. J. Phys. Chem. A 2000, 104, 5161. (31) Barbatti, M.; Aquino, A. J. A.; Lischka, H. Mol. Phys. 2006, 104, 1053. (32) Coe, J. D.; Levine, B. G.; Martínez, T. J. J. Phys. Chem. A 2007, 111, 11302−11310. (33) Tao, H.; Allison, T. K.; Wright, T. W.; Stooke, A. M.; Khurmi, C.; van Tilborg, J.; Liu, Y.; Falcone, R. W.; Belkacem, A.; Martínez, T. J. J. Chem. Phys. 2011, 134, 244306. (34) Kuhlman, T. S.; Glover, W. J.; Mori, T.; Møller, K. B.; Martínez, T. J. Faraday Discuss. 2012, 157, 193. (35) Mori, T.; Glover, W. J.; Schuurman, M. S.; Martínez, T. J. J. Phys. Chem. A 2012, 116, 2808. (36) Sellner, B.; Barbatti, M.; Muller, T.; Domcke, W.; Lischka, H. Mol. Phys. 2013, 111, 2439. (37) Győ rffy, W.; Shiozaki, T.; Knizia, G.; Werner, H.-J. J. Chem. Phys. 2013, 138, 104104. (38) MacLeod, M. K.; Shiozaki, T. J. Chem. Phys. 2015, 142, 051103. (39) Vlaisavljevich, B.; Shiozaki, T. J. Chem. Theory Comput. 2016, 12, 3781−3787. (40) Olsen, J.; Roos, B. O.; Jorgensen, P.; Jensen, H. J. A. J. Chem. Phys. 1988, 89, 2185−2192. (41) Malmqvist, P.-A.; Rendell, A.; Roos, B. O. J. Phys. Chem. 1990, 94, 5477−5482. (42) Ma, D.; Li Manni, G.; Gagliardi, L. J. Chem. Phys. 2011, 135, 044128. (43) Ma, D.; Li Manni, G.; Olsen, J.; Gagliardi, L. J. Chem. Theory Comput. 2016, 12, 3208−3213. (44) Roos, B. O. Adv. Chem. Phys. 1987, 69, 399−445. (45) Shepard, R. Adv. Chem. Phys. 1987, 69, 63. (46) Coe, J. D.; Martínez, T. J. J. Am. Chem. Soc. 2005, 127, 4560. (47) Levine, B. G.; Martínez, T. J. Annu. Rev. Phys. Chem. 2007, 58, 613. (48) Chan, G. K.-L.; Sharma, S. Annu. Rev. Phys. Chem. 2011, 62, 465. (49) Mazziotti, D. A. Chem. Rev. 2012, 112, 244. (50) Fosso-Tande, J.; Nguyen, T.-S.; Gidofalvi, G.; DePrince, A. E. J. Chem. Theory Comput. 2016, 12, 2260−2271. (51) Thomas, R. E. J. Chem. Theory Comput. 2015, 11, 5316−5325. (52) Li Manni, G.; Smart, S. D.; Alavi, A. J. Chem. Theory Comput. 2016, 12, 1245−1258. (53) Pople, J. A.; Santry, D. P.; Segal, G. A. J. Chem. Phys. 1965, 43, S129. (54) Segal, G. A. Semiempirical methods of electronic structure calculation; Plenum Press: New York, 1977. (55) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902. (56) Granucci, G.; Persico, M.; Toniolo, A. J. Chem. Phys. 2001, 114, 10608. (57) Toniolo, A.; Granucci, G.; Martínez, T. J. J. Phys. Chem. A 2003, 107, 3822. (58) Toniolo, A.; Olsen, S.; Manohar, L.; Martínez, T. J. Faraday Discuss. 2004, 127, 149−163. (59) Ciminelli, C.; Granucci, G.; Persico, M. Chem. - Eur. J. 2004, 10, 2327−2341. (60) Hohenstein, E. G.; Luehr, N.; Ufimtsev, I. S.; Martínez, T. J. J. Chem. Phys. 2015, 142, 224103. (61) Snyder, J. W.; Hohenstein, E. G.; Luehr, N.; Martínez, T. J. J. Chem. Phys. 2015, 143, 154107. (62) Fales, B. S.; Levine, B. G. J. Chem. Theory Comput. 2015, 11, 4708−4716. (63) Snyder, J. W.; Fales, B. S.; Hohenstein, E. G.; Levine, B. G.; Martínez, T. J. J. Chem. Phys., submitted for publication. (64) Hohenstein, E. G.; Bouduban, M. E. F.; Song, C.; Luehr, N.; Ufimtsev, I. S.; Martínez, T. J. J. Chem. Phys. 2015, 143, 014111. (65) Hohenstein, E. G. J. Chem. Phys. 2016, 145, 174110. (66) Toniolo, A.; Ben-Nun, M.; Martínez, T. J. J. Phys. Chem. A 2002, 106, 4679−4689. (67) Punwong, C.; Martínez, T. J.; Hannongbua, S. Chem. Phys. Lett. 2014, 610−611, 213−218.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Edward G. Hohenstein: 0000-0002-2119-2959 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

■

REFERENCES

Support for this project was provided by the Martin & Michele Cohen Fund for Science and PSC-CUNY Award #68875-00 46, jointly funded by The Professional Staff Congress and The City University of New York. Computational resources were provided through a Research Cluster Grant from Silicon Mechanics: award number SM-2015-289297. E.G.H. would like to thank Prof. Todd Martı ́nez and Dr. Robert Parrish for helpful discussions and a critical reading of the manuscript.

(1) Olivucci, M.; Ragazos, I. N.; Bernardi, F.; Robb, M. A. J. Am. Chem. Soc. 1993, 115, 3710−3721. (2) Garavelli, M.; Bernardi, F.; Olivucci, M.; Vreven, T.; Klein, S.; Celani, P.; Robb, M. A. Faraday Discuss. 1998, 110, 51−70. (3) Worth, G. A.; Robb, M. A. Adv. Chem. Phys. 2002, 124, 355−431. (4) Levine, B. G.; Martínez, T. J. J. Phys. Chem. A 2009, 113, 12815− 12824. (5) Virshup, A. M.; Punwong, C.; Pogorelov, T. V.; Lindquist, B. A.; Ko, C.; Martínez, T. J. J. Phys. Chem. B 2009, 113, 3280−3291. (6) Plasser, F.; Barbatti, M.; Aquino, A. J. A.; Lischka, H. Theor. Chem. Acc. 2012, 131, 1073. (7) Chattoraj, M.; King, B. A.; Bublitz, G. U.; Boxer, S. G. Proc. Natl. Acad. Sci. U. S. A. 1996, 93, 8362−8367. (8) Tolbert, L. M.; Solntsev, K. M. Acc. Chem. Res. 2002, 35, 19−27. (9) Runge, E.; Gross, E. K. U. Phys. Rev. Lett. 1984, 52, 997−1000. (10) Casida, M. E. In Recent Advances in Density Functional Methods; Chong, D. P., Ed.; World Scientific: Singapore, 1995; Chapter TimeDependent Density Functional Response Theory for Molecules, p 155. (11) Foresman, J. B.; Head-Gordon, M.; Pople, J. A.; Frisch, M. J. J. Phys. Chem. 1992, 96, 135. (12) Geertsen, J.; Rittby, M.; Bartlett, R. J. Chem. Phys. Lett. 1989, 164, 57−62. (13) Koch, H.; Jørgensen, P. J. Chem. Phys. 1990, 93, 3333−3344. (14) Stanton, J. F.; Bartlett, R. J. J. Chem. Phys. 1993, 98, 7029−7039. (15) Shavitt, I.; Bartlett, R. J. Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory; Cambridge University Press: Cambridge, 2009. (16) Yarkony, D. R. Acc. Chem. Res. 1998, 31, 511. (17) Levine, B. G.; Ko, C.; Quenneville, J.; Martínez, T. J. Mol. Phys. 2006, 104, 1039. (18) Levine, B. G.; Coe, J. D.; Martínez, T. J. J. Phys. Chem. B 2008, 112, 405. (19) Shao, Y.; Head-Gordon, M.; Krylov, A. I. J. Chem. Phys. 2003, 118, 4807−4818. (20) Minezawa, N.; Gordon, M. S. J. Phys. Chem. A 2009, 113, 12749−12753. (21) Zhang, X.; Herbert, J. M. J. Chem. Phys. 2014, 141, 064104. (22) Stanton, J. F.; Gauss, J. J. Chem. Phys. 1994, 101, 8938−8944. (23) Nooijen, M.; Bartlett, R. J. J. Chem. Phys. 1995, 102, 3629−3647. (24) Sears, J. S.; Sherrill, C. D.; Krylov, A. I. J. Chem. Phys. 2003, 118, 9084−9094. (25) Zhang, X.; Herbert, J. M. J. Chem. Phys. 2015, 143, 234107. (26) Krylov, A. I. Annu. Rev. Phys. Chem. 2008, 59, 433−462. (27) Werner, H.-J. Adv. Chem. Phys. 1987, 69, 1. (28) Roos, B. O. Acc. Chem. Res. 1999, 32, 137−144. (29) Ben-Nun, M.; Martínez, T. J. J. Am. Chem. Soc. 2000, 122, 6299. 1144


Article

Journal of Chemical Theory and Computation (68) Goyal, P.; Hammes-Schiffer, S. J. Phys. Chem. Lett. 2015, 6, 3515−3520. (69) Goyal, P.; Schwerdtfeger, C. A.; Soudackov, A. V.; HammesSchiffer, S. J. Phys. Chem. B 2016, 120, 2407−2417. (70) Hohenstein, E. G. J. Am. Chem. Soc. 2016, 138, 1868−1876. (71) Wesolowski, T. A.; Warshel, A. J. Phys. Chem. 1993, 97, 8050− 8053. (72) Manby, F. R.; Stella, M.; Goodpaster, J. D.; Miller, T. F. J. Chem. Theory Comput. 2012, 8, 2564−2568. (73) Fux, S.; Jacob, C. R.; Neugebauer, J.; Visscher, L.; Reiher, M. J. Chem. Phys. 2010, 132, 164101. (74) Goodpaster, J. D.; Ananth, N.; Manby, F. R.; Miller, T. F. J. Chem. Phys. 2010, 133, 084103. (75) Goodpaster, J. D.; Barnes, T. A.; Miller, T. F. J. Chem. Phys. 2011, 134, 164108. (76) Huang, C.; Pavone, M.; Carter, E. A. J. Chem. Phys. 2011, 134, 154110. (77) Nafziger, J.; Wu, Q.; Wasserman, A. J. Chem. Phys. 2011, 135, 234101. (78) Goodpaster, J. D.; Barnes, T. A.; Manby, F. R.; Miller, T. F. J. Chem. Phys. 2012, 137, 224113. (79) Lie, G. C.; Clementi, E. J. Chem. Phys. 1974, 60, 1275. (80) Lie, G. C.; Clementi, E. J. Chem. Phys. 1974, 60, 1288. (81) Lee, D.; Furche, F.; Burke, K. J. Phys. Chem. Lett. 2010, 1, 2124− 2129. (82) Lee, D.; Burke, K. Mol. Phys. 2010, 108, 2687−2701. (83) Kim, M.-C.; Sim, E.; Burke, K. J. Chem. Phys. 2011, 134, 171103. (84) Kim, M.-C.; Sim, E.; Burke, K. Phys. Rev. Lett. 2013, 111, 073003. (85) Kim, M.-C.; Sim, E.; Burke, K. J. Chem. Phys. 2014, 140, 18A528. (86) Granucci, G.; Toniolo, A. Chem. Phys. Lett. 2000, 325, 79−85. (87) Slavíček, P.; Martínez, T. J. J. Chem. Phys. 2010, 132, 234102. (88) Shu, Y.; Hohenstein, E. G.; Levine, B. G. J. Chem. Phys. 2015, 142, 024102. (89) The variational theorem does not apply to comparisions between energy expectation values correspondng to wave functions with different numbers of active electrons and orbitals once the DFT correction has been introduced. (90) In the usual CASCI method, many definitions of the reference state are possible that will lead to the same final wave function. The typical starting point occupies all the core orbitals and adds the active electrons. It is also possible to pick the reference state such that the occupied and active orbitals are fully occupied. The excess electrons would then be annihilated to obtain the desired Slater determinants. These starting points do not produce equivalent results once the DFT correction is added. (91) Fales, B. S.; Shu, Y.; Levine, B. G.; Hohenstein, E. G. Manuscript in preparation. (92) Yamaguchi, Y.; Osamura, Y.; Goddard, J. D.; Schaefer, H. F. A New Dimension to Quantum Chemistry: Analytic Derivative Methods in Ab Initio Molecular Electronic Structure Theory; Oxford University Press: Oxford, 1994. (93) Moscardó, F.; Muñoz-Fraile, F.; Pérez-Jiménez, A. J.; PérezJordá, J. M.; San-Fabián, E. J. Phys. Chem. A 1998, 102, 10900−10902. (94) Abia, L. P.; Pérez-Jordá, J. M.; San-Fabián, E. J. Mol. Struct.: THEOCHEM 2000, 528, 59−64. (95) Gräfenstein, J.; Cremer, D. Mol. Phys. 2005, 103, 279−308. (96) Grimme, S. Chem. Phys. Lett. 1996, 259, 128−137. (97) Grimme, S.; Waletzke, M. J. Chem. Phys. 1999, 111, 5645. (98) Malcolm, N. O. J.; McDouall, J. J. W. J. Phys. Chem. 1996, 100, 10131−10134. (99) Malcolm, N. O. J.; McDouall, J. J. W. J. Phys. Chem. A 1997, 101, 8119−8122. (100) Malcolm, N. O. J.; McDouall, J. J. W. Chem. Phys. Lett. 1998, 282, 121−127. (101) McDouall, J. J. W. Mol. Phys. 2003, 101, 361−371. (102) Gusarov, S.; Malmqvist, P.-A.; Lindh, R.; Roos, B. O. Theor. Chem. Acc. 2004, 112, 84−94.

(103) Gusarov, S.; Malmqvist, P.-A.; Lindh, R. Mol. Phys. 2004, 102, 2207−2216. (104) Li Manni, G.; Carlson, R. K.; Luo, S.; Ma, D.; Olsen, J.; Truhlar, D. G.; Gagliardi, L. J. Chem. Theory Comput. 2014, 10, 3669− 3680. (105) Ghosh, S.; Sonnenberger, A. L.; Hoyer, C. E.; Truhlar, D. G.; Gagliardi, L. J. Chem. Theory Comput. 2015, 11, 3643−3649. (106) Hoyer, C. E.; Ghosh, S.; Truhlar, D. G.; Gagliardi, L. J. Phys. Chem. Lett. 2016, 7, 586−591. (107) Gutlé, C.; Savin, A. Phys. Rev. A: At., Mol., Opt. Phys. 2007, 75, 032519. (108) Miehlich, B.; Stoll, H.; Savin, A. Mol. Phys. 2010, 102, 527− 536. (109) Garza, A. J.; Jiménez-Hoyos, C. A.; Scuseria, G. E. J. Chem. Phys. 2014, 140, 244102. (110) Ufimtsev, I. S.; Martínez, T. J. J. Chem. Theory Comput. 2008, 4, 222−231. (111) Ufimtsev, I. S.; Martínez, T. J. J. Chem. Theory Comput. 2009, 5, 1004−1015. (112) Ufimtsev, I. S.; Martínez, T. J. J. Chem. Theory Comput. 2009, 5, 2619−2628. (113) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; Shamasundar, K. R.; Adler, T. B.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hesselmann, A.; Hetzer, G.; Hrenar, T.; Jansen, G.; Köppl, C.; Liu, Y.; Lloyd, A. W.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklass, A.; O’Neill, D. P.; Palmieri, P.; Peng, D.; Pflüger, K.; Pitzer, R.; Reiher, M.; Shiozaki, T.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M. MOLPRO, version 2012.1, a package of ab initio programs. 2012. See http://www.molpro.net (accessed Feb 8, 2017). (114) Hättig, C.; Weigend, F. J. Chem. Phys. 2000, 113, 5154. (115) Furche, F.; Ahlrichs, R.; Hättig, C.; Klopper, W.; Sierka, M.; Weigend, F. WIRES: Comput. Mol. Sci. 2014, 4, 91−100. (116) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157. (117) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. J. Mol. Graphics Modell. 2006, 25, 247−260. (118) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (119) Dutta, A.; Sherrill, C. D. J. Chem. Phys. 2003, 118, 1610. (120) Chou, P. T.; Chen, Y. C.; Yu, W. S.; Chou, Y. H.; Wei, C. Y.; Cheng, Y. M. J. Phys. Chem. A 2001, 105, 1731. (121) Takeuchi, S.; Tahara, T. J. Phys. Chem. A 2005, 109, 10199. (122) Schriever, C.; Barbatti, M.; Stock, K.; Aquino, A. J. A.; Tunega, D.; Lochbrunner, S.; Riedle, E.; de Vivie-Riedle, R.; Lischka, H. Chem. Phys. 2008, 347, 446. (123) Higashi, M.; Saito, S. J. Phys. Chem. Lett. 2011, 2, 2366. (124) Lee, J.; Joo, T. Bull. Korean Chem. Soc. 2014, 35, 881−885. (125) Zhou, M.; Zhao, J.; Cui, Y.; Wang, Q.; Dai, Y.; Song, P.; Xia, L. J. Lumin. 2015, 161, 1−6. (126) Barbara, P. F.; Brus, L. E.; Rentzepis, P. M. J. Am. Chem. Soc. 1980, 102, 5631−5635. (127) Laermer, F.; Elsaesser, T.; Kaiser, W. Chem. Phys. Lett. 1988, 148, 119−124. (128) Brewer, W. E.; Martinez, M. L.; Chou, P.-T. J. Phys. Chem. 1990, 94, 1915−1918. (129) Frey, W.; Laermer, F.; Elsaesser, T. J. Phys. Chem. 1991, 95, 10391. (130) Potter, C. A. S.; Brown, R. G.; Vollmer, F.; Rettig, W. J. Chem. Soc., Faraday Trans. 1994, 90, 59−67. (131) Lochbrunner, S.; Wurzer, A. J.; Riedle, E. J. Chem. Phys. 2000, 112, 10699. (132) de Vivie-Riedle, R.; De Waele, V.; Kurtz, L.; Riedle, E. J. Phys. Chem. A 2003, 107, 10591. (133) Lochbrunner, S.; Wurzer, A. J.; Riedle, E. J. Phys. Chem. A 2003, 107, 10580. 1145


Article

Journal of Chemical Theory and Computation (134) Chang, S. M.; Tzeng, Y. J.; Wu, S. Y.; Li, K. Y.; Hsueh, K. L. Thin Solid Films 2005, 477, 38−41. (135) Chang, S. M.; Hsueh, K. L.; Huang, B. K.; Wu, J. H.; Liao, C. C.; Lin, K. C. Surf. Coat. Technol. 2006, 200, 3278−3282. (136) Sun, D.; Fang, J.; Yu, G.; Ma, F. J. J. Mol. Struct.: THEOCHEM 2007, 806, 105−112. (137) Barbatti, M.; Aquino, A. J. A.; Lischka, H.; Schriever, C.; Lochbrunner, S.; Riedle, E. Phys. Chem. Chem. Phys. 2009, 11, 1406− 1415. (138) Mohammed, O. F.; Luber, S.; Batista, V. S.; Nibbering, E. T. J. J. Phys. Chem. A 2011, 115, 7550−7558. (139) Kungwan, N.; Plasser, F.; Aquino, A. J. A.; Barbatti, M.; Wolschann, P.; Lischka, H. Phys. Chem. Chem. Phys. 2012, 14, 9016− 9025. (140) Luber, S.; Adamczyk, K.; Nibbering, E. T. J.; Batista, V. S. J. Phys. Chem. A 2013, 117, 5269−5279. (141) Kim, C. H.; Joo, T. Phys. Chem. Chem. Phys. 2009, 11, 10266. (142) Lee, J.; Kim, C. H.; Joo, T. J. Phys. Chem. A 2013, 117, 1400− 1405. (143) Schriever, C.; Lochbrunner, S.; Ofial, A. R.; Riedle, E. Chem. Phys. Lett. 2011, 503, 61. (144) Aquino, A. J. A.; Lischka, H.; Hättig, C. J. Phys. Chem. A 2005, 109, 3201−3208. (145) Aquino, A. J. A.; Plasser, F.; Barbatti, M.; Lischka, H. Croat. Chem. Acta 2009, 82, 105−114. (146) Coe, J. D.; Martínez, T. J. Mol. Phys. 2008, 106, 537−545. (147) Christiansen, O.; Koch, H.; Jørgensen, P. Chem. Phys. Lett. 1995, 243, 409. (148) Most implementations of CC2 use either conventional or density-fitted representations of the electron repulsion integrals scale and, as a result, scale as 6(N 5). Further reductions in the scaling could be realized through the introduction of localized molecular orbitals or the use of sparsity in the atomic orbital basis. (149) CC2 computations were performed using Turbomole with 8 threads on two 4-core AMD Opteron 2376 processors clocked at 2.3 GHz. CASCI computations were performed with TeraChem using 8 threads on two 6-core Intel Xeon E5-2640 processors clocked at 2.5 GHz and 8 Nvidia GeForce GTX 680 GPUs. (150) Coe, J. D.; Martínez, T. J. J. Phys. Chem. A 2006, 110, 618− 630. (151) Kyrychenko, A.; Herbich, J.; Izydorzak, M.; Wu, F.; Thummel, R. P.; Waluk, J. J. Am. Chem. Soc. 1999, 121, 11179−11188. (152) Kijak, M.; Zielińska, A.; Thummel, R. P.; Herbich, J.; Waluk, J. Chem. Phys. Lett. 2002, 366, 329−335. (153) Nosenko, Y.; Kyrychenko, A.; Thummel, R. P.; Waluk, J.; Brutschy, B.; Herbich, J. Phys. Chem. Chem. Phys. 2007, 9, 3276−3285. (154) Nosenko, Y.; Kunitski, M.; Riehn, C.; Thummel, R. P.; Kyrychenko, A.; Herbich, J.; Waluk, J.; Brutschy, B. J. Phys. Chem. A 2008, 112, 1150−1156. (155) Kyrychenko, A.; Stepanenko, Y.; Waluk, J. J. Phys. Chem. A 2000, 104, 9542−9555. (156) Kyrychenko, A.; Waluk, J. J. Phys. Chem. A 2006, 110, 11958− 11967. (157) Kungwan, N.; Daengngern, R.; Piansawan, T.; Hannongbua, S.; Barbatti, M. Theor. Chem. Acc. 2013, 132, 1397. (158) Kina, D.; Nakayama, A.; Noro, T.; Taketsugu, T.; Gordon, M. S. J. Phys. Chem. A 2008, 112, 9675−9683.

1146


Improved Complete Active Space Configuration Interaction Energies

Recommend Documents