Communication: Adjusting charge transfer state energies for configuration interaction singles: Without any parameterization and with minimal cost

In a recent article, we showed that configuration interaction singles (CIS) has a systematic bias against charge-transfer (CT) states: CT vertical excitation energies are consistently too high (by 1-2 eV) as compared with non-CT energies [J. E. Subotnik, J. Chem. Phys. 137, 071104 (2011)]. We now show that this CIS error can be corrected approximately by performing a single Newton-Raphson step to reoptimize orbitals, thus establishing a new set of orbitals which better balances ground and excited state energies. The computational cost of this correction is exactly that of one coupled-perturbed Hartree-Fock calculation, which is effectively the cost of the CIS calculation itself. In other words, for twice the computational cost of a standard CIS calculation, or roughly the same cost as a linear-response time-dependent Hartree-Fock calculation, one can achieve a balanced, size-consistent description of CT versus non-CT energies, ideally with the accuracy of a much more expensive doubles CIS(D) calculation.


I. INTRODUCTION
Configuration interaction singles (CIS) is the simplest and most intuitive approach for constructing excited electronic states. A CIS wavefunction has the form | CIS ⟩ = ia t a i | a i ⟩ and depends on two sets of variables: (i) the choice of occupied ("i") and virtual ("a") orbitals, and (ii) the choice of amplitudes for singles excitations ("t a i "). In a standard calculation, the choice of orbitals is dictated by a Hartree-Fock (HF) calculation to ensure orthogonality to the ground-state, and the amplitudes are chosen variationally by diagonalizing the Hamiltonian in the basis of single excitations A, Although it is well known that CIS does not recover accurate vertical excitation energies from the ground state, 1 CIS is often good enough to predict accurate rates of electronic excitation transfer between non-CT excited states. 2,3 By implication, this means that CIS often does a decent job of predicting relative energies between non-CT excited states. Other attractive features of CIS include: (i) it is variational; (ii) it is computationally cheap; (iii) it recovers the correct −1/r asymptotic behavior of CT states that comes about because of the Coulombic attraction between electron attachments and detachments. 4 For all of the reasons above, over the last year, our research group has attempted to use CIS theory to consider electron transfer events between excited states, though we have had little success. As we showed in a recent publication, 5 even though CIS recovers the correct −1/r asymptotic behavior of a) Electronic mail: subotnik@sas.upenn.edu.
CT states, CIS excitation energies are highly biased against CT states, shifted usually by 1-2 eV. Thus, the relative energies between CT and non-CT states are unreliable, and quite often the very ordering of CT and non-CT excited states is incorrect with CIS. Now, the standard alternative to CIS for large systems is time-dependent density functional theory 6 (TD-DFT), a method that typically obtains better vertical excitation energies than CIS for non-CT states. Unfortunately, however, current implementations of TD-DFT fail miserably for CT states because the methods do not recover the correct −1/r asymptotic behavior, 4 which leads to CT excitation energies that are often many eV too low (and getting worse for larger systems). [7][8][9][10][11] This failure of TD-DFT stems from using approximate exchange-correlation functionals, 12,13 and Tozer and co-workers 14,15 have argued that TD-DFT errors can be correlated in general with a measure of charge-transfer (though this is not always true 16 ). To correct the CT problem in TD-DFT, cutting-edge research in quantum chemistry is creating new long-range corrected (LRC) TD-DFT functionals that add in exact Hartree-Fock exchange at long distances by partitioning the Coulomb operator. 17-20 LRC functionals are a creative approach to blend together DFT functionals (that underestimate CT state energies) with CIS theory (for which CT state energies are overestimated, but with the correct asymptotic behavior). In the future, it will be interesting to see whether LRC-TD-DFT functionals can give a correct and robust description of both CT and non-CT excited state energies. As with all DFT development, there is no systematic way to improve accuracy in general.
Rather than exploring TD-DFT, the goal of this communication is to provide a simple approach for correcting CIS energies to give a balanced description of CT versus non-CT states. While the accuracy of CIS can always be improved by using an expanded configuration interaction subspace (i.e., including doubles, à la CISD or CIS(D), 21,22 etc.), we will show below that one can find the correct balance between CT and non-CT states simply by reoptimizing orbitals. We emphasize that we do not address here the intrinsically poor vertical excitation energies of CIS, which arise from not including electron-electron correlation. We also do not address errors in the relative energies of non-CT states if doubles corrections are strong and truly unavoidable. Instead, here we intend only to improve relative excitation energies between CT and non-CT states, with the aim in mind of using CIS to model electron transfer between excited states in the near future.

II. THEORY
The CIS wavefunctions are optimized with respect to amplitudes -∂E CIS /∂t a i = 0 -but CIS wavefunctions are certainly not optimized with respect to choice of orbitals. Using standard analytical gradient theory, 23,24 one can parameterize the space of orbital rotations using the anti-symmetric generator of the orthogonal group. Starting with a set of atomic orbitals {χ µ } and a fixed initial set of orthonormal molecular orbitals, {φ 0 i }, φ 0 i = µ χ µ C 0 µi , all possible choices of molecular orbitals are parameterized by orbital coefficients C, Using the standard definition of the Fock matrix F pq , the CIS energy has the form Differentiating with respect to θ pq , we find that We will show below that Y ai is much larger for CT states than for non-CT excited states. While this result in itself is not surprising, the scale of the energy difference is rather surprising: looking at Y alone is often enough to discern a CT state from a non-CT state. Given this result, one is tempted to correct standard CIS states by accounting for orbital optimization. The simplest correction is to take a Newton-Raphson step. Thus, we expand the excitation energy to second order (6) and searching for the optimal , we find Now, unfortunately, this approach has two drawbacks. First, the method requires us to invert the second-derivative matrix individually for each excited state, rather than all at once. Second, there is no guarantee that ∂ 2 E CIS ∂θ ai ∂θ bj will be or should be a positive definite matrix, and as such, the method may be unsatisfactory. After all, if the Hessian were degenerate, the method would be unstable. To that end, a reasonable solution may be to replace the second-derivative in Eqs. (6)-(8) with the HF second-derivative; our intuition here is that for larger displacements in the choice of orbitals (i.e., large ), the HF term in Eq. (4) will often dominate over the three other terms. Thus, our final expression for the orbital-optimized configuration interaction singles (OO-CIS) energy correction is which corresponds to a first-order (in ) perturbative wavefunction In Sec. III, we will show that Eq. (10) yields a strong correction for CT states that is in approximate agreement with CIS(D) 21,22 in the limit of long-range charge transfer. Conveniently, Eq. (10) can be solved using only one z-vector call 25 to invert ∂ 2 E HF ∂θ ai ∂θ bj for all excited states at once.

III. RESULTS: 2-(4-(PROPAN-2-YLIDENE) CYCLOHEXYLIDENE)MALONONITRILE (PYCM)
To test the theory above, we have studied the PYCM molecule from Ref. 5 (shown in the inset of Fig. 1) using a developers' version of the Q-Chem software package. 26 For a set of 500 different nuclear geometries, we have computed the first 12 excited states, amongst which there is almost always at least one CT state: the electron donor is the dimethyl alkene group and the electron acceptor is the dicyano group. Electronic absorption experiments have shown that the CT state should be the lowest lying excited state, 27 but CIS calculations in vacuum drastically overestimate the vertical excitation energy of this CT state, ranking the CT state always between the third and seventh excited states, in disagreement with experiment. Of course, in solution, a CT state would be stabilized by the surrounding solvent molecules -which we ignore in our calculation -but solvent effects are not large enough to account for the discrepancy. In Ref. 5, we showed that CIS(D) gives a strong correction, lowering the energy of the CT state (making it the first or second excited state), and that correction is to a good approximation proportional to the excited state relative dipole moment. If OO-CIS is a valid theory, we expect that it should behave similarly. We now present three pieces of evidence in the method's favor.
As our first piece of evidence, in Fig. 1, we show that the orbitals are indeed far less optimized for CT-states as compared with non-CT states. In this figure, we have plotted a histogram of the trace of Y T Y (Eq. (5)), separated according to the relative dipole moment of the excited states | ⃗ µ rel | = | ⃗ µ − ⃗ µ HF |. Non-CT states are colored blue, and CT states are colored red. Of our 6000 calculations (500 × 12), we identify 498 CT states, all of which have relatively large norms for Y ai .
Next, in Fig. 2, we present a scatter plot of both E CIS(D) = E CIS(D) − E CIS and E OO-CIS = E OO-CIS − E CIS versus the magnitude of the relative dipole moment of each excited state | ⃗ µ rel |. For non-CT states (on the left-hand side), E CIS(D) follows no obvious pattern and can be positive or negative. For OO-CIS, the energy correction is always negative, because the HF Hessian is positive definite, but the energy correction is very small for non-CT states. By contrast, for CT states (on the right-hand side of Fig. 2), both CIS(D) and OO-CIS are proportional to the relative dipole moment, with nearly the same slope! Thus, our OO-CIS approach approximately recovers the CIS(D) correction for strong CT states, up to a constant shift in energy; according to OO-CIS, the CT state is usually the first, second, or third excited states.
Finally, to strengthen our argument, we present in Fig. 3  ) are colored red. While the correction energies appear uncorrelated for non-CT states, they become highly correlated for CT states. In fact, a linear fit of these points shows a slope of 0.79 for the red and green points, and a slope roughly 0.96, much closer to 1, for the red dots alone. In total, this data indicates that, for at least one molecule with CT excited states, OO-CIS gives a very meaningful energy correction, quite comparable to CIS(D).

IV. DISCUSSION AND CONCLUSIONS
The data above demonstrate that orbital optimization is crucial for correcting the relative energies of CT versus non-CT excited states. This fact suggests many novel avenues for exploration. First, it will be crucial in the future to implement and analyze the more rigorous corrections in Eqs. (7) and (8), where the step size in orbital space is determined by the second derivative of the CIS energy (rather than the HF energy). Preliminary evidence suggests that Eq. (8) may yield larger corrections (in absolute value) for CT state energies as compared with Eq. (10), but that the CIS Hessian may not always be positive definite. Is it reasonable to expect that excited states must correspond to a local minima or are saddle points also physical? Second, one may wonder if similar treatments of orbital optimization can be applied to TD-DFT (in the Tamm-Dancoff approximation); the CT problem is not yet completely solved in the framework of TD-DFT, and orbital optimization may yield new insight. Third and finally, OO-CIS wavefunctions are of the form in Eq. (11) and are thus orthonormal only to first order, and orthogonal to the groundstate only to zeroth order. Future research must analyze the properties of these wavefunctions beyond energetics. Moreover, do the doubles corrections found in Eq. (11) match explicitly the doubles corrections found in CIS(D), which are clearly necessary to describe CT states 5 ? In summary, we have shown that a simple one-step orbital optimization (Eq. (10)) yields a meaningful correction to CIS excitation energies for CT states, comparable in fact to a CIS(D) correction. The approach is clearly size-consistent and the computational cost of this correction is minimal: inverting ∂ 2 E HF ∂θ ai ∂θ bj has the same cost (approximately) as the CIS calculation itself, which is an order of magnitude cheaper than CIS(D). Thus, for only twice the computational cost of standard CIS, or exactly the same cost as TDHF, the OO-CIS approach rebalances the relative excitation energies of CT versus non-CT CIS states. Because this cost is so minimal, we expect Eq. (10) will likely become a standard component of all future CIS calculations.