Absolutely Localized Projection-Based Embedding for Excited States
Xuelan Wen Daniel S. Graham Dhabih V. Chulhai Jason D. Goodpaster


[
    Grace Hui Yang
    July 2019
==================





We present a quantum embedding method that allows for the calculation of local excited states embedded in a Kohn-Sham density functional theory (DFT) environment.  Projection-based quantum embedding methodologies provide a rigorous framework for performing DFT-in-DFT and wave function in DFT (WF-in-DFT) calculations.  The use of absolute localization, where the density of each subsystem is expanded in only the basis functions associated with the atoms of that subsystem, provide improved computationally efficiency for WF-in-DFT calculations by reducing the number of orbitals in the WF calculation.  In this work, we extend absolutely localized projection-based quantum embedding to study localized excited states using EOM-CCSD-in-DFT and TDDFT-in-DFT. The embedding results are highly accurate compared to the corresponding canonical EOM-CCSD and TDDFT results on the full system, with TDDFT-in-DFT frequently more accurate than canonical TDDFT.  The absolute localization method is shown to eliminate the spurious low-lying excitation energies for charge transfer states and prevent over delocalization of excited states. Additionally, we attempt to recover the environment response caused by the electronic excitations in the high-level subsystem using different schemes and compare their accuracy.  Finally, we apply this method to the calculation of the excited state energy of green fluorescent protein and show that we systematically converge to the full system results.  Here we demonstrate how this method can be useful in understanding excited states, specifically which chemical moieties polarize to the excitation.   This work shows absolutely localized projection-based quantum embedding can treat local electronic excitations accurately, and make computationally expensive WF methods applicable to systems beyond current computational limits. 





§ INTRODUCTION
 

Understanding photo-induced processes in chemical, biological, and material systems necessitates accurate descriptions of electronic excited states.
Linear-response time-dependent density functional theory (LR)-TDDFT<cit.>—the excited state extension to ground state density functional theory (DFT)<cit.>—is oftentimes the method of choice for excited states of medium to large molecules (up to a few hundreds of atoms) due to its reasonable cost scaling (formally 𝒪(N^4)<cit.>, or 𝒪(N^3)<cit.> depending on implementation) and accuracy.
However, implementations of TDDFT rely on approximate exchange-correlation functionals and typically utilize the adiabatic approximation.<cit.>
In general, TDDFT is found to describe valence excitations quite accurately, however, these approximations lead to well-studied shortcomings of TDDFT in accurately describing Rydberg states, long range charge-transfer excitations, conical intersections, and double excitations<cit.>.


Excited states properties may be calculated with high accuracy using ab initio correlated wave function (WF) methods.
Equation of motion coupled cluster with singles and doubles (EOM-CCSD)<cit.> has been shown to be accurate for the excited states of many systems<cit.>. 
However, the computational cost of EOM-CCSD (scaling as 𝒪(N^6)) restricts its usage to about 50 or fewer atoms in a moderate basis.
Methods that extend the applicability of EOM-CCSD to larger systems include the use of local orbitals,<cit.> restricted virtual spaces,<cit.> and multiscale approaches.<cit.>
Multiscale approaches treat different regions of the molecule with methods of varying cost and accuracy to reflect their importance.<cit.>
For example, one may use EOM-CCSD to describe the important region of the molecule while using molecular mechanics (MM), in EOM-CCSD-in-MM methods<cit.>, or DFT, in EOM-CCSD-in-DFT methods<cit.>, to describe the remainder of the molecule.
Such approaches take advantage of the high accuracy of EOM-CCSD and the low scaling of MM or DFT to go beyond the size/accuracy limits of any one method alone.

Density functional theory embedding provides a formally exact framework for performing multiscale calculations. This approach has been shown to be accurate for combining density functional theory with wave function theory.
The key choice in density function theory embedding is the treatment of the non-additive kinetic energy which can be treated using either approximate kinetic functionals<cit.>, or optimized effective potentials<cit.>.
An alternative approach is to employ projection operators to avoid the use of a non-additive kinetic potential<cit.>.  

Projection-based embedding is a density functional theory embedding methodology where one uses projection operators to ensure that the orbitals of two or more regions (or subsystems) of the molecule are mutually orthogonal. This projection operator enforces the Pauli exclusion principle between subsystems and avoids the need to use non-additive kinetic energy functionals. 
Two such projection operators have since been used: (1) the parameter dependent operator (hereinafter the μ operator);<cit.> and (2) the Huzinaga operator.<cit.>



Extension of projection-based embedding for excited states has been explored by several groups with applications to EOM-CCSD<cit.>, complete active space self-consistent field (CASSCF),<cit.> linear response TDDFT (LR-TDDFT),<cit.> real-time TDDFT (RT-TDDFT)<cit.>, and density matrix renormalization group (DMRG)<cit.>.
Of primary concern is the polarization response of the environment subsystem B to an excitation in the high-level subsystem A. Chulhai and Jensen introduced an environment response term derived from the μ operator.<cit.>
Bennie and co-workers accounted for this response by systematically including important environment orbitals in the description of the high-level subsystem,<cit.> though this increases the computational cost of the EOM-CCSD calculations.





Chulhai and Goodpaster<cit.> have presented an absolute localization projection-based embedding method based on the Huzinaga operator.<cit.>
In the absolute localization method, the embedded wave function region is restricted to the basis functions of the embedded subsystem only, thereby significantly reducing the computational cost.
In this paper, we extend the absolute localization projection-based embedding method to describe excited states using both EOM-CCSD and TDDFT.
We show that this approach is able to accurately describe localized excitations in subsystems that are divided across covalent bonds, eliminate spurious low-lying charge-transfer excitations in TDDFT, and efficiently account for environment polarization.












§ THEORY
 



 §.§ Absolute Localization Projection-Based Embedding
  

In projection-based embedding methods, the total system is divided via the electron density such that for two subsystems 


    γ_tot^0 = γ_A^0 + γ_B^0


where γ_A^0, γ_B^0, and γ_tot^0 are the ground state density matrices of subsystems A and B, and the total system, respectively.
The Fock matrix of subsystem A embedded in the environment of B is


    𝐟^A-in-B = 𝐡 + 𝐉[γ_A^0 + γ_B^0] + 𝐯_xc[γ_A^0 + γ_B^0] + 𝐏^B


where 𝐡 is the total one-electron Hamiltonian that contains the kinetic and nuclear potentials, 𝐉 is the Coulomb potential, 𝐯_xc is the exchange-correlation (XC) potential, and 𝐏^B is the projection operator that enforces inter-subsystem orbital orthogonality.

The Manby and Miller groups proposed the use of the parameter-dependent level-shift projection operator (hereinafter the μ operator) as<cit.>


    P^B =  μ( 𝐒^ABγ_B^0 𝐒^BA)


where μ is a very large parameter (usually 10^6), and 𝐒^AB(BA) are the AO overlap matrices between subsystems A and B. This operator has since seen widespread application and development,<cit.> including basis set truncation<cit.> and extensions to describe excited states.<cit.>

Alternatively, the Huzinaga operator<cit.> may also be used in DFT embedding;<cit.> it is defined as



    𝐏^B = - 1/2( 𝐅^ABγ_B^0 S^BA + S^ABγ_B^0 𝐅^BA)


where 𝐅^AB(BA) are the AB (or BA) block of the total Fock matrix.
The Huzinaga operator yields accurate results when each subsystem density is expanded in only basis functions associated with that subsystem.<cit.>
In these cases, the basis sets of A and B are mutually exclusive and the electrons of A and B are restricted to their respective subsystems only. This strategy, termed “absolute localization”, has been shown to increase accuracy through error cancellation.<cit.> This allows one to significantly reduce the size of basis sets used to describe subsystems A and B, and therefore significantly reduce the computational cost of the WF calculation.

Standard EOM-CCSD has computational scaling of the order O^2_A+BV^4_A+B, where O indicates the number of occupied orbitals and V indicates the number of virtual orbitals.
When subsystem density is expanded in the basis functions of the full system (or full-system basis, as shown later), WF-in-DFT calculations scale as O^2_AV^4_A+B.
When the subsystem density is expanded in only the basis functions associate with itself (or subsystem basis), WF-in-DFT calculations scales as O^2_AV^4_A. Therefore, the size of subsystem B could be as large as the computational limit of DFT for EOM-CCSD-in-DFT calculations.



 §.§ Absolutely-Localized Embedding for Excited States
   


In this paper, we separate the environment polarization into two parts—a ground state polarization and a polarization response—and we introduce two approaches to include the latter in absolute localization embedding for excited states.
For the ground-state polarization, we use the ground state embedding potential—obtained from a ground state absolute localization embedding method as in ref. —for excited states calculations.
That is, the embedding potential is built using the DFT ground-state densities of subsystems A and B (γ_A^0 and γ_B^0) and included through the use of a modified core Hamiltonian (𝐡^A-in-B) defined as


    𝐡^A-in-B[γ_A^0,γ_B^0] = 𝐡 + 𝐉[γ_tot^0] - 𝐉[γ_A^0] + 𝐯_xc[γ_tot^0] - 𝐯_xc[γ_A^0] + 𝐏^B


This modified core Hamiltonian is then included in all subsequent correlated excited states WF calculations.
However, as previously mentioned, this strategy ignores the polarization response of the environment subsystem B.


In order to account for the polarization response of the environment, we propose two approaches. The first is a state-averaged approach that polarizes subsystem B to an average of the ground state and excited state(s) densities of subsystem A. This state-averaged polarization is included in a self-consistent method as follows (the algorithm is include in Figure <ref>):




(1) Perform absolute localization DFT-in-DFT to obtain self consistent densities γ_A^0 and γ_B^0, and embedded core Hamiltonian 𝐡^A-in-B[γ_A^0,γ_B^0].

(2) Perform TDDFT on subsystem A to obtain excited state(s) densities γ_A^i, where i denotes the ith excited state. The transition density matrix γ_A^0 → i is defined in the eqn. 23 from ref. .

(3) Calculate the state-averaged electron density of subsystem A (γ_A^avg). This state-averaged electron density, assuming that we are interested in the first n excited states, is defined as 

    γ^avg_A = 1/n+1∑_i=0^n γ_A^i


(4) Re-optimize the ground state of subsystem B using the new embedded core Hamiltonian 𝐡^B-in-A[γ_B^0,γ^avg_A] to obtain a polarized electron density γ^0,pol_B, where 0 denotes the ground state. This step allows subsystem B to be polarized by γ^avg_A.

(5) Re-optimize the ground state of subsystem A with 𝐡^A-in-B [γ^avg_A, γ^0,pol_B] to obtain a polarized γ_A^0,pol. This step allows the electron density in A to be polarized by γ^0,pol_B.

(6) Repeat steps 2–5 using γ_A^0,pol instead of γ_A^0 until both γ^0,pol_A and γ^0,pol_B are converged; the convergence criteria is set to 10^-3 for the Frobenius norm of the change in the density matrices in this paper. 

(7) Perform excited states WF (or TDDFT) calculations using the embedded core Hamiltonian 𝐡[γ_A^0,pol, γ_B^0,pol] with the converged γ_A^0,pol and γ_B^0,pol. 


The second approach to account for the environment polarization response, only applicable to embedded excited states WF methods, is by using a correction term obtained by TDDFT on the full system, defined as


    ω̃_WF-in-DFT[Ψ^A;γ_A^0, γ_B^0,γ_KS-DFT^0]   
        =ω_TDDFT[γ_KS-DFT^0]  - ω_TDDFT-in-DFT[γ_A^0, γ_B^0]  
        +ω_WF-in-DFT[Ψ^A;γ_A^0, γ_B^0]


where ω̃_WF-in-DFT is the corrected embedded WF excitation energy, ω_WF-in-DFT is the embedded WF excitation energy (with ground-state polarization or state-average polarization), ω_TDDFT-in-DFT is the embedded TDDFT excitation energy, and ω_TDDFT is the canonical TDDFT excitation on the full system with KS-DFT ground density γ_KS-DFT^0.
This correction term (ω_TDDFT - ω_TDDFT-in-DFT) is meant to account for missing environment polarization response for WF methods at the TDDFT level and is further explored in Section <ref>.




§ COMPUTATIONAL DETAILS


All embedded TDDFT calculations were perform in our Quantum Solid State and Molecular Embedding (QSoME) code<cit.>, which utilizes the Python-based Simulations of Chemistry Framework (PySCF).<cit.>
All full-system and embedded EOM-CCSD calculations were performed using Molpro 2015.1<cit.>.
For embedded EOM-CCSD calculations, the embedded core Hamiltonian 𝐡^A-in-B is calculated and then exported from QSoME to Molpro for the EOM-CCSD calculations only.
All full-system TDDFT calculations and Natural transition orbitals (NTOs) were performed using Gaussian 16 program<cit.>. 
 








The geometries of long-chain hydrocarbons with different functional groups (Section <ref>) were obtained using the M06 functional<cit.> and the cc-pVTZ basis set<cit.> using Gaussian 16 program<cit.>; all other calculations used the cc-pVDZ basis set.<cit.>
Geometries for acrolein in water in Section <ref> were generated using molecular dynamics, and the aug-cc-pVDZ basis set<cit.> was used for all embedding calculations.
The geometry for acrolein in two water molecules in Section <ref> was taken from ref , and the aug-cc-pVDZ basis set was used for all calculations;
the state-average density includes the ground state and the first excited state.
Geometries for cis-7HQ (Section <ref>) were taken from ref  and the aug-cc-pVDZ basis set was used for all calculations.
The geometry for GFP (Section <ref>) was taken from ref ; the B3LYP functional<cit.> and def2-TZVP basis set<cit.> were used for all calculations.




§ RESULTS AND DISCUSSION






 §.§ Choice of Basis for the Subsystem
  





First, we will focus on the choice of basis functions to use in embedding calculations for excited states embedded correlated wave function methods — namely, equation of motion coupled-cluster with single and double excitations (EOM-CCSD).<cit.> EOM-CCSD scales as N_occ^2N_vir^4, where N_occ and N_vir are the number of occupied and virtual orbitals, respectively. Therefore, embedding calculations, even if performed in the full basis of the full system still benefit since the number of occupied orbitals in the high-level subsystem N^A_occ may be significantly reduced (N^A_occ≪ N_occ). One may further improve the efficiency of such methods by employing basis set truncation schemes<cit.>, thereby reducing the available virtual orbitals of the subsequent WF calculation (N^trunc_vir < N_vir).<cit.>




With the Huzinaga operator, we have shown that we can truncate the basis to only include those that describe the atoms in subsystem A.<cit.> This truncation restricts the electrons of subsystem A to that subsystem only, and we label this method as “absolute localization.” For ground state properties of both molecular and extended materials, this absolute localization strategy for WF-in-DFT provides similar or improved accuracy to other embedding schemes that also use projection operators at a significantly reduced cost. In this subsection, we will continue to explore the performance of Huzinaga and absolute localization embedding strategies for embedded excited states WF-in-DFT.



In order to compare these embedding strategies, we will look at their accuracy and efficiency at correctly reproducing the localized excited states in six long-chain organic molecules. We chose the three lowest singlet excitations in decanoic acid, decan-1-amine, dec-1-ene, and decanal, and the two lowest excitations in 1-chlorodecane and decan-1-ol.
The NTOs for these excitations are shown in Figure <ref>. As seen in this figure, these excitations are localized on the functional groups and up to five carbon atoms of the alkane chains.
The third excited states in 1-chlorodecane and decan-1-ol are delocalized (or global) excitations; these results are included in the Supporting Information Figure S1 and Figure S2.
We chose long-chain molecules in order to systematically increase the size of the high-level subsystem while testing the accuracy of the methods. Note that, in each case, we will be cutting across covalent bonds. Ding and co-workers examined similar long-chain organic molecules, albeit for embedded TDDFT calculations, in order to show that the μ strategy is systematically improvable.<cit.>



Our embedding results are shown in Figure <ref>. In panel A, we show the unsigned errors of the first excitation of decanol where the basis set for subsystem A includes all the basis functions of the system (full-system, orange) or only the basis functions associated with the atoms of subsystems A (absolute localization, green).  In panel B, we show the mean absolute errors (M.A.E.) of all 16 local excitations mentioned in Figure  <ref>, the values for all excitations are provided in the Supporting Information Table S1 and Table S2. 

Chemical accuracy is reached — that is, less than 1 kcal/mol or 0.04 eV — with four embedded carbons for the absolute localization method, compared to only three carbons for the full-system basis. 
This is because, for the HONTO → LUNTO transitions shown in Figure <ref>, the LUNTOs are generally more delocalized than the HONTOs.
For the embedding calculations that use the full-system basis, the embedded EOM-CCSD uses the virtual orbital space of the full system, while the absolute localization strategy uses the virtual orbital space of the high-level subsystem only. Although both strategies have the same size of occupied space, the larger virtual space considered in the full basis strategy makes the results slightly more accurate. We want to emphasize here that absolute localization<cit.> is not simply an extreme case of basis set truncation,<cit.> as both subsystems A and B are also solved self-consistently in the absolute localized basis at the DFT-in-DFT level, similar to “freeze-and-thaw” subsystem DFT methods.<cit.>

The computation time of the EOM-CCSD calculations, the most time-consuming step in each of these methods, are shown in Figure <ref>C and  <ref>D. Although the absolute localization strategy is marginally less accurate, the subsequent embedded EOM-CCSD calculation are still 1–3 orders of magnitude faster due to the significantly reduced number of virtual orbitals.

Typically, N^A_vir≤ N^trunc_vir < N_vir, where N^A_vir is the number of virtual orbitals for subsystem A in the absolute localized basis. This significant reduction in the number of virtual orbitals results in a more efficient EOM-CCSD calculation since EOM-CCSD scales as N^4 with respect to the number of virtual orbitals, compared to only N^2 with respect to the number of occupied orbitals.
A strategy to reduce the number of virtual orbitals in a correlated wave function treatment of excited states is crucial to cost saving, especially where large basis sets are required.
For example, Silva-Junior and co-workers found that at least a TZVP basis set are required for low-lying valence excited states when using correlated wave function methods,<cit.> while additional diffuse functions should be used for Rydberg states.<cit.> In situations where such large basis sets are required, the speedup of our absolute localization embedding strategy would be even more significant than the cases shown in this subsection, which used cc-pVDZ basis set.  Given the significant computational cost savings with a relatively minor increase in error, we will use the absolute localization embedding strategy in the remainder of this paper.










 §.§ Eliminating Spurious Charge-Transfer States in TDDFT
 

It is a well-known that TDDFT, particularly for local and semi-local exchange-correlation (XC) functionals, may result in artificial low-lying charge-transfer (CT) excited states.<cit.>
This is caused by the incorrect 1/R asymptotic behavior of approximate XC kernels.<cit.>

Long-range corrected functionals—such as CAM-B3LYP<cit.>, LRC-ωPBEh<cit.>, LRC-ωPBE<cit.>, and ωB97XD<cit.>—correct for this issue by enforcing the correct 1/R asymptotic behavior at large distances.
Embedded TDDFT<cit.> is uniquely suited to solve this problem. Using the subsystem basis restricts the electron density of each subsystem on their own spatial regime, and therefore forbids the charge-transfer excitations between subsystems.
As a result, the artificially low long-range CT excitations between different weakly interacting subsystems will be eliminated.
Similarly, absolute localization can eliminate the spurious low-lying CT states since the subsystems are described with the subsystem basis. However, global and real CT electronic excitations will also be eliminated for the same reason.



Figure <ref> shows the results for the lowest n→π^* and π→π^* transitions for s-trans acrolein surrounded by different amounts of water molecules.
In Figure <ref> (A), we observe that embedded B3LYP gives similar results to B3LYP for the n →π^* transition, while CAM-B3LYP predicts higher excitation energies.
In Figure <ref> (B), both embedded B3LYP and CAM-B3LYP predict higher excitation energies for the π→π^* transition.
In Figure <ref> (C), we give the number of excited states to be calculated to obtain the bright π→π^* transition.
For gas-phase acrolein, the first excited state is n →π^* transition, and the second excited state is a bright π→π^* transition responsible for the experimental absorption peak.
As we include more water molecules into the B3LYP and CAM-B3LYP calculations, the number of excited states that needs to be calculated before we obtain the bright π→π^* transition increases due to spurious low-lying CT states between the water molecules and arcolein; this was also observed in other similar systems, for example in uracil solution.<cit.>
The spurious CT problem is severe in B3LYP, is largely mitigated by CAM-B3LYP, and is completely avoided by embedded B3LYP.
Figure <ref> (D) lists the product of the numbers of the occupied and virtual orbitals N_occ-vir, where the scaling of TDDFT calculation is 𝒪(N^3_occ-vir). 

Based on the NTO analysis in Figure S4 and S5, the n →π^* transition is localized with both B3LYP and CAM-B3LYP.
For this reason, embedded B3LYP gives similar result to B3LYP.
As we will further explore in Section <ref>, the environmental response is also weak; therefore, the ground-state embedding potential provides an accurate description for the environmental polarization resulting the close agreement between B3LYP and embedded B3LYP.
The CAM-B3LYP predict higher excitation energy than B3LYP because it has larger HF exchange at the long distance.<cit.>
The NTO analysis in Figure S6 shows that the π→π^* transition has a small delocalized character with CAM-B3LYP and large delocalization with B3LYP.
CAM-B3LYP results in higher excitation energies than B3LYP due to the larger HF exchange at long distances, while embedded B3LYP increases the excitation energies (as compared to B3LYP) by forcing the π→π^* transition to be localized on the acrolein molecule, which increases the spatial overlap between HONTO and LUNTO<cit.>.
As seen in Figures S3–S5, the NTOs of embedded B3LYP are fully localized on acrolein, and the calcualted oscillator strengths are much larger than those of B3LYP.







The example of acrolein in water suggests that absolute localization can simultaneously treat local excitation accurately, eliminate the spurious low-lying charge-transfer excitation effectively, and speed up the calculation significantly. In addition, it enforces the localization of π→π^* that are fictitiously delocalized by B3LYP.  Therefore, embedded TDDFT appears to be an excellent tool to study the excited states of solvated molecules.  










 §.§ Accounting for Environment Polarization
  



For excited states calculations, there are two types of polarization to consider: (1) The polarization of the environment to the ground-state of the high-level subsystem A; and (2) the polarization response of the environment to excitations in the high-level subsystem A. The former is exactly captured within projection-based embedding methods, at least when using the full system basis. For the latter in a quantum embedding framework, there has so far been three general strategies:
(i) accounting for the coupled response of the environment within a TDDFT framework;<cit.>
(ii) including important environment orbitals in the description of the high-level subsystem A within correlated wave function methods like EOM-CCSD;<cit.> or
(iii) using a state specific embedding potential.<cit.>
Within the embedded TDDFT framework, the response of the environment is often accounted for using a coupled frozen-density embedding term added to the TDDFT kernel.<cit.> This coupled response approach have also been extended to projection-based TDDFT embedding with the μ projection operator.<cit.> Within conventional projection-based excited states embedding,<cit.> all of the environment virtual orbitals, unless some basis set truncation scheme<cit.> is used, are included in the subsequent correlated WF calculation of the high-level subsystem. The polarization response of the environment therefore only requires inclusion of selected important environment occupied orbitals into the high-level subsystem.<cit.> However, as far as we are aware, there is no a priori method to determine which environment orbitals are important, and current methods<cit.> rely on a guess-and-check scheme.



We examine the effects of various types of polarization on the lowest n→π^* transition in an acrolein+2water system, which was previously studied using μ projection-based excited states embedding.<cit.> In Figure <ref>,
We will first consider the effects of ground-state polarization on this transition. For embedded TD-HF-in-HF (or TD-PBE0-in-PBE0), we find that ground-state polarization (purple bars) — that is, absolute localization embedding without consideration for the environment polarization response — already accounts for ∼90% of the missing polarization in the canonical full system results (blue bars) when compared to the isolated acrolein results (“no polarization;” green bars).  For embedded EOM-CCSD-in-HF(or PBE0), ground-state polarization similarly captures ∼90% of the polarization, regardless of the quality of the low-level method (HF or PBE0 in the case). We will discuss the effect of the lower level method in the next subsection. The small fraction of  missing polarization — around 10% in the case of the hydrogen-bonded acrolein+2water examined here — is due to a missing environment response. This is because there are no orbitals, neither occupied nor virtual, from environment that are included in the high-level subsystem A in the absolute localization method.



In order to include the polarization response of the environment, we will use two approaches: (1) Using a state-averaged approach; and (2) Using a TDDFT correction. In the state-averaged approach (pink bars in Figure <ref>), subsystem B is polarized to the average density of the ground-state and first excited state of subsystem A. However, we find that this approach lowers the excitation energy, which in this case leads to increased error.    

Possible reasons for the behaviour of state-average embedding potential are:
1. Huzinaga operator shifts the occupied orbitals of the environment to positive energy levels, but does not removed them, and therefore creates fictitious virtual orbitals that might be involved in the subsequent excited-state calculations. 

This unrealistic relaxation on excitation energies might be exaggerated by the optimizing the subsystem densities using state-average embedding potential self-consistently. 
2. The errors associated with embedding are more complex than simply the lack of environmental response.  Recent work has shown that the accuracy of different embedding strategies for excitation energies are largely similar with or without environmental response.<cit.>  It is possible that there are errors associated with the embedding potential that get exaggerated when the state-averaged embedding potential is used; therefore, although the electrostatic potential is more accurate due to environmental response, the overall potential is of similar quality.  Regardless of the exact reason, this suggests that including the response of the environment does not always lead to better accuracy.  



















The second approach to account for the environment polarization response is by using a correction term obtained by TDDFT; this correction is defined in eqn. <ref>. We explore a similar correction term for ground state energies in ref. , where the correction accounted for the errors introduced by using the absolutely localized basis — errors introduced, for example, by preventing electrons from moving between subsystems. Analogously, by performing an embedded excitation calculation of subsystem A in the absolutely localized basis, we lose the contributions to the excitation from coupling to subsystem B. In ref. , this missing coupling is described in terms of orbital transitions — transitions involving both occupied and virtual orbitals of subsystem B are missing — while in ref. , this missing coupling is described in terms of a coupled environment response — changes to the density of subsystem B when subsystem A is excited are missing. The correction presented in eqn. <ref> corrects for this missing coupling (or environment polarization response) using the canonical TDDFT excitation on the full system. Since we are interested in cases where canonical EOM-CCSD on the full system is prohibitively expensive and the embedded EOM-CCSD is the most time consuming step of our methodology, there are systems in which a full system TDDFT calculation is an insignificant cost to pay to include the environment response; of course, for some systems this cost could become intractable.  

The results with this TDDFT correction are also shown in Figure <ref> and included in the Supporting Information Table S3. We observe that this correction result in similar accuracy, regardless of whether we add the correction to the ground-state polarization or state-average polarization approaches. With this correction, we also observe a small dependence on the quality of the lower level method, with the TDHF corrected EOM-CCSD-in-HF results 0.03 eV higher than the benchmark canonical EOM-CCSD excitation energy, while the TDDFT corrected EOM-CCSD-in-DFT results exactly, errors less than 0.01 eV, reproduced the benchmark. These corrected errors are small, given that canonical TDHF and TD-PBE0 have errors of 0.76 and 0.30 eV, respectively, compared to canonical EOM-CCSD. Our results are similar in accuracy to the those in ref. for a smaller embedded EOM-CCSD orbital space. 

We will continue to examine the accuracy of this correction approach and its dependence on the quality of the lower level method in next subsection.

In summary, our results suggest that ground-state polarization of the environment provides adequate accuracy.  Additional schemes to include the polarization response of the environment do not always increase the accuracy.  However, the error associated with ground-state polarization for EOM-CCSD-in-PBE0 or EOM-CCSD-in-HF is only 0.015 eV and 0.019 eV, respectively.  Therefore, this absolute localization scheme with ground-state polarization is very accurate for WF-in-DFT embedding, and the lack of need of polarization response of the environment means the algorithm remains highly computationally efficient.  



 §.§ Quality of Low Level Method
  









To examine the effects of the choice exchange-correlation functional on the quality of the EOM-CCSD-in-DFT embedding, we looked at the solvent shifts of cis-7-hydroxyquinoline (cis-7HQ) due to eight hydrogen-bonded complexes: water, ammonia, methanol, formic acid, two water molecules, ammonia-water-water trimer, ammonia-water-ammonia trimer, and ammonia-ammonia-water trimer; solvent shifts are measured with respect to the gas phase 7HQ. These results are shown in Figure <ref> (and tabulated in Supporting Information Table S4), where we plot the solvation energy shifts of the lowest π→π^* excitation (Δω_π→π^*) compared to the canonical EOM-CCSD results.
In this figure, we examine two exchange-correlation functionals, LDA and CAM-B3LYP, and show the results for canonical TDDFT, embedded TDDFT, embedded EOM-CCSD with ground-state polarization, and embedded EOM-CCSD with ground-state polarization and the TDDFT correction in eqn. <ref>.

Of the methods explored, canonical TDDFT is the least accurate, with MAEs of 930 and 406 cm-1 for LDA and CAM-B3LYP, respectively. These errors also appear to be dependent on the size of the solvent molecule(s).
The solvent molecule(s) from smallest to largest solvation shifts (right to left in Figure <ref>), are CH3OH, H2O, HCOOH, NH3, 2H2O, NH3-H2O-NH3, NH3-H2O-H2O, and NH3-NH3-H2O.
Fradelos and co-workers examined the same systems<cit.> and suggested that the poor performance of TDDFT is most likely caused by the incorrect 1/R asymptotic behavior of the LDA potentials; they found that the statistical average of orbital potentials (SAOP)<cit.>, an asymptotically correct exchange-correlation potential, gives better results than the LDA and PW91 functionals in both embedded and canonical TDDFT calculations.
Therefore, this size-depedency of the canonical TDDFT errors can be attributed to erroneously delocalizing—and therefore lowering the energy of—the excited state over the entire system<cit.>.



Canonical TD-CAM-B3LYP, which has a factional 1/R dependence,<cit.> provides significantly more accurate excitations energies compared to TD-LDA.


Embedded TDDFT-in-DFT, with ground-state polarization only, performs better than canonical TDDFT, with MAEs of 539 and 153 cm-1 for LDA and CAM-B3LYP, respectively.  As explored in Figure <ref>, this is due to restricting the transition to the high-level subsystem only through absolute localization, and preventing the delocalization of the π^* state over the full system. This restriction results in more accurate localized excitations, regardless of the asymptotic behavior of the exchange-correlation functionals. Therefore, for these types of localized solvated excitations, embedded TDDFT can be more accurate and significantly cheaper than canonical TDDFT.  



For the embedded EOM-CCSD-in-DFT calculations with ground-state polarization only, we find that the choice of low-level functional—either LDA or CAM-B3LYP, shown in Figure <ref>A—does not affect the results significantly; this was also observed for acrolein in Figure <ref>.
The EOM-CCSD-in-LDA and EOM-CCSD-in-CAMB3LYP have MAEs of 183 and 204 cm-1, respectively.

Given that embedded EOM-CCSD-in-DFT calculations do not significantly depend on the choice of the DFT exchange-correlation functional, it suggests that the reason the embedded TDDFT-in-DFT calculations do depend on the exchange-correlation functional is due to the TDDFT functional and the not the environment DFT functional.  Therefore, it is not surprising that CAM-B3LYP provides significant better accuracy compared to LDA.  We can conclude that the low level functional does not significant influence our embedded excitation energies.  This is consistent with our previous results, which found that the choice of exchange-correlation functionals have a minimal effect on the ground state properties of absolutely localized embedded WF-in-DFT.<cit.>







Including state-average polarization (results not shown) only marginally improves (∼ 20 cm-1) on the ground-state polarization results. The additional cost of multiple embedded TDDFT-in-DFT calculations for the state-average approach therefore makes it less desirable. However, for the TDDFT correction in eqn. <ref>, we find that the results of EOM-CCSD-in-CAMB3LYP are significantly improved (MAE of 50 cm-1), while the EOM-CCSD-in-LDA results remain the same (MAE of 209 cm-1).
Because this correction depends on the canonical TDDFT on the full system, the corrected EOM-CCSD-in-DFT are therefore dependent on the quality of the XC functional used.
With this TDDFT correction, we also observe that there is no size-dependence in the solvation energy shifts, unlike what was observed for ground-state polarization EOM-CCSD-in-DFT.
We ascribe this to the larger polarization response of the larger solvation molecules; this response is not accounted for by the ground-state embedding but is accounted for by the TDDFT correction.




§ A SYSTEMATICALLY IMPROVABLE METHOD: GREEN FLUORESCENT PROTEIN
  




In this subsection, we will show how the absolute localization method is systematically improvable by examining a 161-atom model of green fluorescent protein (GFP). The geometry is taken from ref , which contains the neutral p-hydroxybenzylidene-imidazolinone chromophore (chro), nine amino acid residues including T62, Q69, Q94, R96, H148, V150, T203, S205, and E222, and four water molecules.



We systematically increased the size of the subsystem A in two ways: (i) by successively adding the H2O, S205, E222-/R96+, and 3H2O residues to the chromophore; or (ii) by successively adding the same moieties but in the reverse order.
We performed two sets of calculations: (1) isolated TDDFT calculations (no embedding) where the size of the system included only the chromophore and additional moieties and (2) embedded TDDFT-in-DFT (ground-state embedding) where the high-level system was the same size as the isolated TDDFT calcualtions but embedded in the the remainder of the 161-atom model as the environment subsystem B.  We note that these calculations are of similar computational cost since the TDDFT cost is identical.
These results are shown in Figure <ref>, where Figure <ref>A shows the 161-atom model with the chromophore (balls and sticks model), residues that are successively added to the high-level subsystem (opaque tubes model), and residues in the environment (translucent tubes model). The transition of interest, shown in Figure <ref>B, was identified by examining
the calculated oscillator strengths and shapes of the natural transition orbital (NTOs).<cit.>  Our canonical TDDFT using B3LYP on the full 161-atom model showed that the brightest low-lying state (oscillator strength f=0.76 a.u.) was 3.12 eV, which is in good agreement with experimental results of 3.12-3.14 eV for the neutral A-form GFP.<cit.>

The results of successively increasing the size of the high-level subsystem are shown in Figures <ref>C and D for both the cluster TDDFT and embedded TDDFT-in-DFT results. The black line is the canonical TDDFT results on the full 161-atom model.
For the cluster (no embedding) calculations, we observe that the results are erratic; increasing the size of the cluster does not always translate to a more accurate excitation energy.
However, for the embedded calculations, we always obtain a more accurate results with a larger subsystem A.
The fundamental difference between these two (cluster and embedded) models is the presence of the ground-state polarization of the environment.
For embedded TDDFT-in-DFT, the ground-state polarization — the most important contribution to solvent shifts as explored in Figure <ref> — is always included, and adding residues to the high-level subsystem systematically increases the amount of the system which is responding to the excited state.
This systematic improvability of projection-based embedding was previously observed for ground state properties,<cit.> and we show here that it also applies to the excited states properties.

Additionally, embedding calculations allow for a systematic study of the relative importance of the water molecules and residues in the excitation energy.  A series of embedding calculations allow for the decoupling of the environment effects between ground state polarization and excited state polarization.  By including the entire system for all calculations, the environment polarization (in the ground state) is included and by expanding the size of the TDDFT calculation, one can then determine how much each neighboring water molecule or residue polarizes to the excited state.   This allows one to understand the role and relative importance of the environment to the excited state of the chromophore.  

Embedding calculations are robust irrespective of which other water molecules or residues are included in the excited state calculation.   For instance, starting with the chromophore (panel C), inclusion of the 3H_2O water molecules lowers the excitation energy by -0.032 eV.  For the pathway followed in panel B, the inclusion of the 3H_2O water molecules occurs last, and lowers the excitation energy by -0.033 eV, nearly identical.  Therefore, regardless of when that moiety is included, embedding calculations can be used to determine the importance of that moiety in excited state polarization.  

These results, and the results from Figure <ref>, demonstrate that systematic improvability is retained when using projection based embedding for excited states. Which is to say that larger subsystems for the excited state regions will lead to more accurate calculations.  Additionally, our calculations on GFP show that these calculations can be used to determine the relative importance of chemical moieties to the excitation energy in a robust manner.    



§ SUMMARY AND CONCLUSIONS


We explored the behavior of absolute localization projection-based embedding method for excited states. We show that, for localized excited states with embedded EOM-CCSD, the absolutely localized method is marginally less accurate but significantly cheaper than embedding using the basis set for the full system.  Additionally, we show how this method can accurately eliminate spurious low-lying excited states in embedded TDDFT.
We found that the ground-state polarization of the environment, naturally included with embedded excited states calculations, accounts for ∼90% of the excitation energy shifts due to solvents. The remainder was ascribed to environment polarization response.  We proposed two approaches to include this polarization response — a state-average approach and a TDDFT correction approach — and found that the TDDFT correction accurately accounts for the environment response; this TDDFT correction is only applicable to embedded EOM-CCSD, where the full system TDDFT cost is small compared to the cost of the subsequent EOM-CCSD calculation. We showed that this embedding strategy is systematically improvable, and we are always guaranteed more accurate results for larger regions included in the high-level subsystem.  Finally, we showed that embedding could be a useful tool for understanding the relative importance of chemical moieties to excited states.  

Our results show correlated wave function level of accuracy at a significantly reduced cost. Absolutely localized basis<cit.> reduces computational cost by aggressively truncating the basis functions to only those on atoms involved in the excitation. This truncation results in far fewer virtual orbitals, and thus lower computational cost.  We used EOM-CCSD as the excited states wave function method throughout this paper, however, we see no limitations in applying this absolutely localization embedding strategy to other excited states wave function methods as well, like the complete active space with second-order perturbation theory (CASPT2).<cit.> We are currently exploring such embedding strategies.



                                                                                      
This research was carried out within the Nanoporous Materials Genome Center, which is supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences under Award DE-FG02-17ER16362. The authors acknowledge the Minnesota Supercomputing Institute (MSI) at the University of Minnesota and the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, for providing resources that contributed to the results reported within this paper.





NTO analysis and embedding results for the third excited states of chlorodecane and decanol; numerical data for Figure 3; The NTO analysis for acrolein in different number of water molecules with embedded B3LYP, canonical B3LYP, and CAM-B3LYP discussed in Figure <ref>; numerical data for Δω_π→π^* discussed in Figure <ref>; geometries used in this paper.




