Communication: A mean field platform for excited state quantum chemistry

In a nutshell, Hartree-Fock (HF) theory¹ applies the ground state variational principle to a wave function ansatz that includes only those correlations that are absolutely necessary to produce a qualitatively correct description of the electrons in a simple molecular ground state. Indeed, HF theory’s Slater determinant hews as closely as possible to a classical mean field state, in which the particles would be completely uncorrelated, while accommodating the Pauli correlations that must be included when describing electrons. This simplicity keeps HF theory (relatively) affordable as it proceeds to minimize its variational principle and make the energy stationary with respect to changes in the molecular orbital basis. That HF theory has for decades been the central platform on which high-accuracy weak correlation treatments are built^2,3 is a reminder of how valuable a qualitatively correct minimally correlated wave function with relaxed orbitals can be. For strongly correlated ground states in which an unambiguous minimal active space exists, the story is much the same, but with complete active space self-consistent field (CASSCF)^4–7 and its descendants^8–15 replacing HF theory as the basic platform upon which weak correlation methods are built.^16–30 In this paper, we develop an affordable mean field platform for simple excited states while also providing an initial weak correlation treatment and a discussion of how the approach can be generalized to strongly correlated excitations.

Unlike the situation for ground states, even zeroth order descriptions of most excited states require more than Pauli correlations. For example, consider the open shell singlet of a simple HOMO → LUMO excitation. In this state, two opposite-spin electrons are strongly correlated with each other so as not to occupy the same orbital at the same time, implying that a minimally correlated excited state ansatz must incorporate correlations not present in a Slater determinant. An obvious candidate for this job is a configuration interaction singles (CIS) ansatz³¹ which, to stay in line with ground state mean field methods like HF and CASSCF, has had its energy made stationary with respect to relaxations of the molecular orbital basis. While this direction has of course been explored before, previous approaches have tended to approximate the orbital relaxation in order to suppress cost.^32–35 For example, Subotnik’s OO-CIS approach³² achieves its efficiency by employing an incomplete and HF-approximated Newton-Raphson optimization. Inspired by the need for fully relaxed excited state orbitals in situations like charge transfer and core spectroscopy, we seek here to achieve them at the same cost-scaling as CIS by combining recent progress in excited state variational principles^36–39 with a compound application of automatic differentiation (AD). To demonstrate the efficacy of this approach as a mean field platform on which to build excited state correlation methods, we use it as the basis for an excited state analog of second order Møller-Plesset (MP2) theory.³ 

As we are looking to describe an excited state, minimizing the energy

in pursuit of the ground state variational principle will not guarantee convergence to the desired state. Instead, we will seek to minimize the Lagrangian

in which the excited state variational principle⁴⁰ 

guarantees that a sufficiently flexible wave function will converge to the exact (excited) energy eigenstate with energy closest to ω while the Lagrange multipliers $μ→$ guarantee that an energy stationary point with respect to the wave function variables $ν→$ will be achieved even when working with an approximate ansatz. While we could have done without W if we merely wished to ensure that an energy stationary point was reached, its presence guarantees that in the limit of a sufficiently flexible wave function, the global minimum of L will be the stationary point corresponding to the desired excited state. Although we do not pursue it here, one can imagine augmenting W with functions of the dipole moment or other observables in order to differentiate between states that are energetically degenerate.

To make the minimization of this Lagrangian affordable, we must deal both with the difficulty of the Ĥ² term and with the fact that derivatives of L with respect to the variational parameters $ν→$ (e.g., the CI coefficients and the orbital rotation variables) lead to second derivatives of the energy. To avoid Ĥ², we can resolve an identity in the basis of Slater determinants

and, following the approach of coupled cluster (CC) theory,² make an approximation in which we restrict our attention (and the range of the sum) to the most chemically relevant corner of Hilbert space, in this case the span of the closed shell “ground state” determinant and the singles excitations. Note that W need not be evaluated exactly for the theory to operate correctly, as its role is merely to guide the optimization into the correct stationary point. So long as W provides a sufficiently strong nudge to get us close, the Lagrange multiplier term in Eq. (2) will ensure convergence to the stationary point. This type of “nudged” convergence to a stationary point has been seen before in the application of full configuration interaction quantum Monte Carlo to excited states,⁴¹ where a heavily approximated projection operator that could only serve to nudge the state propagation away from the lower states was sufficient to converge the imaginary time evolution onto excited states. When an approximate nudge is insufficient, our formulation allows for a systematic retreat to safety via increasingly accurate evaluations of W, but we stress that in our initial testing this has yet to prove necessary. Indeed, in all cases tested so far, the even more aggressive approximation

leads the minimization of L to converge on to the same stationary point as when truncating the sum in Eq. (4). If this equivalence is maintained after more extensive testing in the future, there would be a strong simplicity argument in favor of employing Eq. (5).

Unlike the Ĥ² terms, the challenge of second energy derivatives can be overcome by a compound use of automatic differentiation (AD),⁴² which ensures that all of the derivatives of a many-input-single-output function can be evaluated for a small constant multiple of the cost of evaluating the function itself. As AD’s use in quantum chemistry is still in its infancy,^43–45 let us briefly explain the principles so that its usefulness for our purposes is clear. First, consider that many complicated functions, such as the CIS energy, can be written as a graph in which each node is one of the four basic binary operations shown in Table I in which two input quantities a and b go in and the output quantity f that comes out may then become one of the inputs for one or more of the other nodes in the graph. One can compute the overall function value $g(x→)$ by traversing the graph, starting from the dangling edges that are the inputs $x→$ and moving forward through all the nodes until the final output $g(x→)$ is reached. Now, if one can afford to store the outputs of each node in the graph, it becomes possible to evaluate all of g’s first derivatives with respect to the elements of $x→$ for a cost that is a small constant multiple of the evaluation cost of g via a sort of reverse traversal of the graph.⁴² Crucially, if one considers a given node and assumes that they already know the partial derivative of g with respect to the output f of that node, then the chain rule and the simple derivative formulas in Table I ensure that the partial derivatives of g with respect to the node’s inputs a and b can be evaluated via three or fewer binary arithmetic operations. Starting at the final node and working backwards, one finds that the number of operations required to get all the derivatives of g with respect to all the intermediates and all the elements of $x→$ is not worse than four times the number of nodes in the graph, and so the cost to get all the derivatives $∂g/∂x→$ is a small constant multiple of the cost to evaluate g itself.

This approach, known as reverse accumulation, provides a straightforward if tedious recipe for constructing a low cost implementation of analytic derivatives, and we may, for example, apply it to the CIS energy E to obtain an efficient function for the energy first derivatives that appear in Eq. (2). Folding this logic over on itself, we recognize that thanks to the dot product of these efficient derivatives with $μ→$⁠, L itself can be implemented as a many-input-single-output function of the variables $ν→$ and $μ→$ whose cost is a small constant multiple of that of E. By a second, compound application of reverse accumulation, we may thus arrive at an implementation that delivers analytic derivatives of L with respect to both $ν→$ and $μ→$ for a cost that is also a small constant multiple of that of a CIS energy evaluation and with an additional memory requirement for intermediate storage that scales only as the square of the system size. Note that, in practice, setting up reverse accumulation is tedious work that can and has been automated in many software packages. In the present study, we have leveraged the machine learning community’s rapid progress in reverse accumulation software by employing the TensorFlow⁴⁶ framework to evaluate our Lagrangian’s derivative vector, the norm of which we minimize by a quasi-Newton approach.⁴⁷ As the cost of a CIS energy evaluation scales as the cost of a Fock matrix build, so does evaluation of L and the necessary derivatives, leading our excited state mean field (ESMF) approach to have the same cost scaling as CIS.

While we have motivated this Lagrangian approach with the prospects of relaxing the orbitals in CIS, the logic supporting its construction is much more general, and indeed L can in principle be employed efficiently with any ansatz for which E and a reasonable approximation to W can be efficiently evaluated. L could, for example, be used as a more rigorous alternative to maximum overlap methods⁴⁸ when optimizing a CASSCF wave function for an individual excited state. Although we are quite curious about this possibility, we do not pursue strongly correlated excited state treatments in this study. Instead, we focus on delivering a fully orbital-relaxed CIS wave function and testing its ability to act as a platform for excited state correlation treatments in the same way HF theory does for ground states.

To this end, we employ the CIS-like ansatz

in which excitations are labeled by alpha (i, a) or beta (ī, ā) indices and the closed shell determinant |0⟩ is included to help the orbital-relaxed excited state better maintain orthogonality to the RHF ground state. The vector of variables $ν→$ that we optimize via L includes the coefficients c₀, c_ia, and c_īā, as well as the elements of the matrix X that defines the orbital rotation operator,

which for the present study we constrain so as to keep the orbitals spin-restricted. Although we do not explore the possibility in the present study, this approach could be generalized to work with a CASSCF wave function by replacing the CIS expansion in Eq. (6) with the CASSCF CI expansion. Given the much greater size of such CI expansions, it may in that case be more effective to use L only for the orbital rotation optimization and to instead rely on modern CI solvers to keep the energy stationary with respect to the CI coefficients.

While the excited state mean field (ESMF) ansatz |Φ⟩ is more flexible than CIS and might therefore be expected to be more accurate, one should remember that, due to the significant effect of weak correlation on energetics, HF theory itself is quite poor quantitatively even when it is qualitatively a good zeroth order wave function. By the same reasoning, total energies and energy differences from ESMF wave functions (which revert to RHF when ω targets the ground state) are not likely to be competitive in accuracy with methods that incorporate correlation effects, such as equation of motion coupled cluster singles and doubles (EOM-CCSD).⁴⁹ 

As in the ground state, we are interested in the ESMF wave function not as a destination in and of itself, but rather as a reliable platform upon which to construct correlation treatments that can reasonably hope to achieve more quantitative accuracy. Considering first the prospects for a coupled cluster theory, we note that due to the natural termination of the Baker-Campbell-Hausdorff (BCH) expansion in traditional, similarity-transformed CC singles and doubles (CCSD),² the usual approach of projecting the CC eigenvalue equation

into the space of low order (in this case internally contracted) excitations from the reference will, as in the ground state theory, lead to a polynomially complex system of equations for the cluster amplitudes. However, due to the fact that bare triples excitations would be present within the internally contracted doubles, the cost scaling of such an approach would, although still polynomial, be substantially higher than in ground state CCSD. More enticing is the prospect of leveraging the fact that the potentially long-range nature of the excitation and the orbital relaxations it induces should already be accounted for in the ESMF reference state, allowing lower-cost local CC approaches^50–53 to focus on what they do best: treating short-ranged weak correlation.

As intriguing as CC methods may be, the principle of Occam’s razor suggests that an analog of the much simpler MP2 theory would be a wiser starting point for investigations into post-ESMF correlation methods. We can arrive at just such a method, which we will denote as ESMP2, by applying standard Rayleigh-Schrödinger perturbation theory to the zeroth order Hamiltonian

where $F^$ is the Fock operator with respect to the ESMF one-body density matrix, $R^=|ΦΦ|$⁠, $P^$ is the projector to the span of |0⟩ and the singles excitations in the ESMF orbital basis, and $Q^=1−P^$⁠. Note that $P^Φ=Φ and Q^Φ=0.$ This choice of Ĥ₀ leads to the MP2-like zeroth order relationship,

and indeed we see that when we set ω so as to target the ground state, |Φ⟩ becomes the RHF state and ESMP2 simplifies to MP2. In the excited state case, as in traditional MP2, the first order wave function contains no determinants with fewer than two excitations, which is a consequence of including the $P^ĤP^$ term in Ĥ₀. We also maintain the relationship E_ESMF = E₀ + E₁ in direct analogy to the MP2 relationship E_HF = E₀ + E₁. As the Fock operator is not diagonal for excited states, the first order amplitudes on doubles and triples are found by inverting $F^−E0$ via the minimal residual Krylov subspace method, for which the MP2-style denominators are an excellent preconditioner. Note that this approach produces a fully excited-state-specific first order wave function, as opposed to the CIS(D) method where the triples are a product of the CIS coefficients and the ground state MP2 amplitudes.^56,57

While we have, for the sake of simplifying development and testing, written our pilot ESMP2 implementation in a fully uncontracted form whose triples part has an N⁷ cost scaling, multiple avenues exist for recovering the N⁵ scaling of traditional MP2. On the one hand, we could exploit sparsity in the ESMF coefficient matrix, which we observe to have only a handful of elements that are not small. Indeed, we have tested this idea by setting all but three of the elements to zero when solving the ESMP2 linear equation and evaluating its second order energy for the molecules discussed below and found that excitation energy predictions are not strongly affected. On the other hand, we could follow the internally contracted approach of Evangelista and co-workers.⁵⁸ 

Table II reveals that, like its HF cousin in the ground state, the ESMF approximation does not confer quantitative accuracy. This behavior can be understood as a direct consequence of our design goal of hewing as closely as possible to a classical mean field theory. By including only those correlations that are absolutely necessary to realize a fermionic excited state (namely, Pauli correlations and the open shell correlation), ESMF is missing all weak correlations and so, like HF, does not produce quantitative energies. The fact that ESMF tends to underestimate excitation energies can also be understood in terms of what correlations are included. Indeed, in creating their open-shell arrangements, ESMF gives each excited state roughly one pair’s worth of electron correlation, and so these states’ energies are biased low compared to that of the closed shell ground state. Thus, in direct analogy to HF theory, the simplicity and mean field nature of ESMF prevents it from delivering accurate energetics on its own, but this of course was not the intention. What is more important is the question of whether ESMF can match HF theory’s ability to act as a foundation for correlation methods, a question that our early results appear to answer strongly in the affirmative.

The data reveal that in water, formaldehyde, stretched lithium fluoride, and neon, ESMP2 rivals EOM-CCSD in accuracy. Its errors are typically at least a factor of five smaller than CIS, and it substantially outperforms CIS(D) in the charge transfer states of stretched LiF. Our primary explanation for this success is that the reference wave function’s mean-field quality of having fully relaxed orbitals places it at a similar “distance” from the correct wave function as for HF in the ground state, with the subtle effects of weak electron correlation being all that is missing. By contrast, CIS(D), which is also inspired by and closely entwined with MP2 theory, is asking its perturbation to capture both correlation effects and orbital relaxation. While the former are typically small in systems that are not strongly correlated, the size of the latter is much more system dependent. Table II shows that in the low lying transitions of water, in which the overall spatial distribution of electrons is not greatly changed, both CIS(D) and ESMP2 are highly accurate. In stretched LiF, however, the transitions convert between ions and neutral atoms, and as one would expect these large charge density changes lead to strong orbital relaxation effects, as revealed by comparing the CIS and ESMF energies. These relaxations are much more difficult to treat perturbatively and, to make matters worse, cause the closed-shell-state MP2 amplitudes that all CIS(D) states rely on to be less appropriate for the open-shell states. Looked at from this perspective, it is not surprising that ESMP2, thanks to its orbital-relaxed reference and fully state-specific perturbation, delivers more uniform accuracy across charge-transfer, valence, and Rydberg states alike. Indeed, when compared to the high-level benchmark of EOM-CC(2,3),⁵⁹ the maximum unsigned error of ESMP2 is less than the mean unsigned error for CIS(D).

We can also compare our results with more recent attempts to provide orbital relaxations for CIS. The OO-CIS method, for example, provides relaxations via a single Newton-Raphson orbital optimization step in which the HF Hessian is used as an approximation to the CIS orbital Hessian.³² As OO-CIS lacks a treatment for weak correlation, it is not surprising that in H₂O, CH₂O, and Ne, its accuracy is for most states better than CIS (which has neither orbital relaxation nor a weak correlation treatment) but worse than CIS(D), EOM-CCSD, and ESMP2, as shown in Table II. The large errors that OO-CIS makes in LiF can be understood as a consequence of its Hessian approximation, as LiF’s smallest HF Hessian eigenvalues are about an order of magnitude smaller than those in the other molecules, and so when this Hessian is inverted to get the OO-CIS orbital relaxation, the resulting Newton-Raphson step is much too large. In this system at least, it appears that the HF Hessian is not an effective surrogate for the CIS Hessian. Unlike OO-CIS and our approach, the variationally orbital-adapted CIS (VOA-CIS) method seeks to provide orbital relaxations through specially chosen expansions of the configuration interaction space.³⁴ However, there is at present no single prescription for defining the expansion (the developers explore at least eight possibilities in their initial paper³⁴) and so instead of making an extensive and hard-to-interpret direct comparison, we will point out the developers’ conclusion that “most of the time, VOA-CIS closely follows CIS(D)” and their data that shows that it is not unusual for VOA-CIS to be in error by between 0.5 eV and 1 eV for single excitation energies.³⁴ In summary, our preliminary testing shows ESMP2, with its inclusion of both orbital relaxation and state-specific correlation, to be closer in its behavior to EOM-CCSD than to previous attempts at augmenting the CIS wave function.

We have presented an excited state mean field theory and investigated its potential as a platform on which to build excited state correlation treatments. Like HF theory, the method relies on making a minimally correlated wave function’s energy stationary with respect to orbital relaxations. While HF does this for a Slater determinant, we do so for a CIS-like wave function in order to accommodate the basic structure of simple excitations, and so the cost of our mean field optimization has the same scaling with system size as CIS. Unlike HF theory, our approach incorporates an excited state variational principle into its Lagrangian so that ground and excited states are treated equally. In exploring the method’s potential as a platform for correlation treatments, we have constructed an excited state analog of MP2 theory and found that, in initial tests, its accuracy rivals that of EOM-CCSD.

Looking to the future, there are a number of important questions to consider about this excited state mean field approach. First, how will it fare when applied to larger and more complicated charge transfer systems, core excitations, and Rydberg states? Second, will its advantages be maintained when paired with more complex and expensive methods such as coupled cluster and the complete active space self-consistent field? Finally, can this minimally correlated, excited-state-specific wave function be usefully employed as a replacement for the Slater determinant in an excited-state-specific generalization of the Kohn-Sham approximation? We look forward to exploring these and other questions that will undoubtedly arise in the context of excited state mean field theory.

We thank Joe Subotnik for many helpful discussions, and we acknowledge funding through the Early Career Research Program of the Office of Science, Office of Basic Energy Sciences, the U.S. Department of Energy, Grant No. DE-SC0017869. Calculations were performed both on our own desktop computers and at the National Energy Research Scientific Computing Center, 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.