A General Approach for Multireference Ground and Excited States using Non-Orthogonal Configuration Interaction

TL;DR

This paper introduces a general method combining non-orthogonal configuration interaction with holomorphic Hartree--Fock solutions to accurately describe multireference ground and excited states across various molecular geometries.

Contribution

It presents a novel protocol using active space SCF metadynamics and complex plane topology to identify and trace multiple Hartree--Fock solutions for multireference state calculations.

Findings

Successfully applied to fluorine dimer dissociation

Demonstrated effectiveness on cyclobutadiene distortion

Provides a scalable approach for multireference states

Abstract

A balanced description of ground and excited states is essential for the description of many chemical processes. However, few methods can handle cases where static correlation is present, and often these scale very unfavourably with system size. Recently, multiple Hartree--Fock (HF) solutions have been proposed as a basis for non-orthogonal configuration interaction (NOCI) to provide multireference ground and excited state energies, although applications across multiple geometries have been limited by the coalescence of HF solutions. Holomorphic HF (h-HF) theory allows solutions to be analytically continued beyond the Coulson--Fischer points at which they vanish but, until now, this has only been demonstrated for small model systems. In this work, we propose a general protocol for computing NOCI ground and excited state energies using multiple HF solutions. To do so, we outline an active space variation of SCF metadynamics that allows a chemically relevant set of HF states to be identified, and describe how these states can be routinely traced across all molecular geometries by exploiting the topology of h-HF solutions in the complex plane. Finally, we illustrate our approach using the dissociation of the fluorine dimer and the pseudo-Jahn--Teller distortion of cyclobutadiene, demonstrating its applicability for multireference ground and excited states.

Tables

Table 1. Table 1: Total energy E tot subscript 𝐸 tot E_{\mathrm{tot}} of the singlet ground state, and first three vertical excitation energies for cyclobutadiene.

	$E_{\mathrm{tot}}{(}^1{B}_{1g})$	${}_{}{}^{3}A_{2g}^{}$	${}_{}{}^{1}A_{1g}^{}$	${}_{}{}^{1}B_{2g}^{}$
NOCI	$-153.72075$	$0.01630$	$0.06987$
CIS	$-153.60216$	$-0.06703$	$0.01840$	$0.22862$
CASSCF $(4,4)$	$-153.70956$	$0.01657$	—	—
sa-CASSCF $(4,4)$	$-153.70721$	$0.01465$	$0.08204$	$0.14176$
EOM-CCSD	$-154.18974$	$-0.02062$	$0.24389$	$0.25355$
EOM-CCSDT	$-154.23741$	$0.00271$	—	—
ex-FCI	$-154.234\left(1\right)$	$0.013\left(3\right)$	$0.058\left(2\right)$	$0.078\left(1\right)$

Equations

$$|\Psi ⟩=\sum _x^m|{}^x\Phi ⟩a_x,$$

$$|{}^x{\varphi }_i⟩=\sum _{\mu }^n{}^x{C}_{\cdot i}^{\mu \cdot }|{\chi }_{\mu }⟩.$$

$$\sum _x^m\left(H^{wx}-E{S}^{wx}\right)a_x=0,$$

$${}^{wx}{P}_i^{\mu \nu }=({}^x\tilde{C}{)}_{\cdot i}^{\mu \cdot }{(}^w{\tilde{C}}^{\ast }{)}_{i\cdot }^{\cdot \nu },$$

$$E=\frac{⟨{\Phi }^{\ast }|\hat{H}|\Phi ⟩}{⟨{\Phi }^{\ast }|\Phi ⟩}.$$

$$\hat{H}=V_{\text{N}}+\sum _i^N\hat{h}\text{quantity}(i)+\sum _{i<j}^N\frac{1}{r_{ij}},$$

$$P^{\mu \nu }=\sum _i^N{C}_{\cdot i}^{\mu \cdot }{C}_{i\cdot }^{\cdot \nu }$$

$$F_{\mu \nu }=h_{\mu \nu }+\sum _{\sigma \tau }^n⟨\mu \sigma ||\nu \tau ⟩P^{\tau \sigma },$$

$$\sum _{\mu \nu }^n{C}_{i\cdot }^{\cdot \mu }{S}_{\mu \nu }{C}_{\cdot j}^{\nu \cdot }={\delta }_{ij},$$

$${\hat{H}}_{\lambda }=V_{\text{N}}+\sum _i^N\hat{h}\text{quantity}(i)+\lambda \sum _{i<j}^N\frac{1}{r_{ij}},$$

$${\varphi }_{\alpha }$$

$${\varphi }_{\beta }$$

Full text

I Introduction:

A balanced treatment of ground and excited states is essential for the description of a wide range of physical processes, including singlet fission,Coto et al. (2015); Mayhall (2016); Zimmerman et al. (2010) electron transfer,Jensen et al. (2018) and primary mechanisms of vision.Barca et al. (2018); Gozem et al. (2014, 2012) However, the presence of strong static correlation in many of these systems — where the single determinant self-consistent field (SCF) description breaks down — presents a challenge for many conventional electronic structure techniques. In principle, full configuration interaction (FCI) can be used to compute the exact energy spectrum for a given basis set. Practically, however, the exponential scaling of FCI limits its applicability to a minority of relatively small chemical systems.

To circumvent the exponential scaling of FCI, numerous approaches have been developed that provide approximate ground state energies. In cases dominated by a single reference determinant, these methods range from the mean-field Hartree–Fock (HF) approximation,Attila Szabo and Neil S. Ostlund (1989) to truncated Configuration Interaction (CI), Coupled Cluster (CC) and Perturbation Theory (PT) approaches.Helgaker et al. (2000); Shavitt and Bartlett (2009) Meanwhile, for excited states the most common single-reference techniques include Time-Dependent Density Functional Theory (TD-DFT),Runge and Gross (1984) uncorrelated CI singles (CIS) Foresman et al. (1992) and Equation of Motion CC (EOM-CC).Stanton and Bartlett (1993); Watts and Bartlett (1994) Of these, EOM-CC provides the most accurate excitation energies, but is usually limited to only single and double excitations (EOM-CCSD) by the prohibitive cost of including full triple excitations.Bartlett and Watts (1995)

In contrast, the range of techniques available for computing ground and excited states in systems with strong multireference character remains much more limited. Multi-configurational SCFHelgaker et al. (2000) (MCSCF) methods provide the prevailing family of techniques, with the Complete Active Space SCF (CASSCF) method being the most widely used approach.Malmqvist and Roos (1989) In CASSCF, an FCI expansion is solved within a pre-defined active orbital space while the reference orbitals are simultaneously optimised. To compute multiple excited states and prevent variational collapse, CASSCF must usually be applied using a state-averaged formalism (sa-CASSCF) where the weighted energy of multiple states is optimised rather than a single target state.Werner and Meyer (1981) Ultimately, however, CASSCF scales exponentially with the size of the active space and remains a challenge for larger systems, despite recent advances using stochasticManni et al. (2016) and selected CI (sCI)Smith et al. (2017) approaches.

In addition to conventional techniques, the multiple mean-field solutions to the HF equations have themselves been proposed as approximations to physical excited states.Barca et al. (2018); Gilbert et al. (2008); Besley et al. (2009) As a non-linear method, HF theory can yield a large number of multiple solutionsFukutome (1975); Stanton (1968); Burton et al. (2018) each representing a Slater determinant built from a bespoke set of optimal HF molecular orbitals. Since the development of SCF metadynamicsThom and Head-Gordon (2008) and the Maximum Overlap Method (MOM),Besley et al. (2009) the identification of multiple SCF solutions has become relatively routine. In many cases, however, multiple HF states can can break the symmetries of the true Hamiltonian, leading to wave functions that do not satisfy the molecular point group symmetry, or are not eigenfunctions of the spin operators $\hat{\mathcal{S}}^{2}$ or $\hat{\mathcal{S}}_{\text{z}}$.Stuber and Paldus (2003)

At the HF level, the onset of strong static correlation is indicated by the break-down of the single determinant approximation. The restricted HF (RHF) ground state in \ceH2, for example, overestimates the binding energy by incorrectly dissociating into a linear sum of the physical “radical” and “ionic” states.Coulson and Fischer (1949) In contrast, the unrestricted HF (UHF) approach provides the correct behaviour at dissociation by allowing electrons of opposite spin to occupy different spatials orbitals and localise on opposing hydrogen atoms.Attila Szabo and Neil S. Ostlund (1989) Moreover, the signature of multireference character — namely several significant reference determinants — manifests directly in the degeneracy of the HF solutions.

To exploit the natural description of static correlation provided by multiple HF solutions, recent research has focussed on using these states as a multireference basis for CI calculations.Thom and Head-Gordon (2009); Sundstrom and Head-Gordon (2014); Yost et al. (2013); Mayhall et al. (2014); Burton and Thom (2016); Oosterbaan et al. (2018); Jensen et al. (2018) Since multiple HF solutions are not required to be mutually orthogonal, this CI takes the form of a non-orthogonal CI (NOCI).Thom and Head-Gordon (2009) By constructing a wave function as a linear combination of non-orthogonal HF solutions, NOCI is able to capture electron correlation with a scaling of $\mathcal{O}(m^{2}\ \text{max}(N^{3},n^{2}))$, where $m$ is the number of determinants used in the expansion, $N$ is the number of electrons, and $n$ is the number of basis functions.Thom and Head-Gordon (2009) In addition, including all degenerate symmetry-broken HF states in the NOCI basis enables the restoration of broken symmetries in a similar style to Projected HF (PHF)Jiménez-Hoyos et al. (2012); Löwdin (1955); Scuseria et al. (2011) and Half-Projected HF (HPHF)Smeyers and Doreste-Suarex (1973); Smeyers and Delgado-Barrio (1974); Cox and Wood (1976) approaches. However, in contrast to the variation-after-projection style of PHF methods, NOCI provides a projection-after-variation formalism where symmetry is partially restored in a single CI expansion. As a result, NOCI provides only approximate symmetry restoration for infinite symmetry groups (e.g. spin), but avoids the challenging non-linearity and self-consistency of PHF methods.

As an inherently multireference approach, NOCI is well-suited to capturing static correlation and can be made size-consistent by ensuring a suitable set of determinants is used.Sundstrom and Head-Gordon (2014) Since determinants are constructed from different sets of orbitals — each optimised individually at the SCF level — NOCI can provide a more balanced treatment of ground and excited states, leading to its application for multi-electron excitations,Sundstrom and Head-Gordon (2014) core excitations,Oosterbaan et al. (2018, 2019) and charge transfer processes.Jensen et al. (2018) Moreover, the non-orthogonality of the NOCI basis can lead to more efficient and compact multideterminantal expansions that provide useful initial guess orbitals for active space techniques,Krausbeck et al. (2014) or trial nodal surfaces for quantum Monte-Carlo techniques.Pathak and Wagner (2018); Landinez Borda et al. (2019) The combination of relatively low-order polynomial scaling, compact multideterminant expansions, and the ability to ensure size-consistency presents NOCI as a promising alternative to CASSCF for treating multireference ground and excited states.

Despite the progress of NOCI methods, their application using symmetry-broken HF states over a range of molecular geometries has, until recently, been limited by the disappearance of HF solutions at Coulson–Fischer points.Coulson and Fischer (1949) At such points, there is a sudden reduction in the size of the NOCI basis set leading to discontinuities in the NOCI energy. To prevent these discontinuities, the holomorphic Hartree–Fock (h-HF) approach has been developed as a method for analytically continuing real HF solutions beyond the points at which they vanish.Hiscock and Thom (2014); Burton and Thom (2016); Burton et al. (2018). In h-HF theory, a new energy function is defined by removing the complex-conjugation of orbital coefficients from the conventional HF equations. This transformation can be seen as a complex analytic continuation of real HF theory and results in a non-Hermitian theory in which the energy is a complex-analytic polynomial of the orbital coefficients. Significantly, the stationary points of the h-HF energy appear to exist across all molecular geometries,Burton et al. (2018) extending with complex orbital coefficients beyond the Coulson–Fischer points where their real counterparts vanish. As a result, h-HF stationary states provide a continuous basis for NOCI calculations and yield continuous potential energy surfaces across all geometries.Burton and Thom (2016); Burton et al. (2018)

Although the combination of h-HF and NOCI has been shown to provide accurate approximations to FCI for small systems,Burton and Thom (2016); Burton et al. (2018) its applicability in more chemically relevant examples is yet to be demonstrated and some practical issues have persisted. In particular, since (in the absence of a magnetic field) the h-HF energy is symmetric along the real orbital coefficient axis, it can be difficult to identify the correct complex direction in which to follow h-HF solutions when their real counterparts vanish. Moreover, as the SCF solution space can become very large,Burton et al. (2018) identifying a suitable basis set of chemically relevant SCF states can prove difficult.

In the current Paper we outline a general approach for applying NOCI in combination with h-HF theory. We describe how, by applying SCF metadynamics in an active orbital space, one can identify a suitable set of chemically relevant determinants. We then demonstrate a simple approach that allows h-HF solutions to be followed around a Coulson–Fischer point and into the complex orbital coefficient plane. Finally, we illustrate our general approach by considering to the dissociation of \ceF2 and the pseudo-Jahn–Teller distortion of cyclobutadiene.

II Theory

II.1 Non-Orthogonal Configuration Interaction

The NOCI wave function is constructed as a linear combination of $m$ mutually non-orthogonal basis states $\{\ket{{}^{x}\Phi}\}$ as

[TABLE]

where we employ the non-orthogonal tensor notation of Head-Gordon et al.Head-Gordon et al. (1998) Each state $\ket{{}^{x}\Phi}$ corresponds to a single Slater determinant constructed from a bespoke set of $N$ occupied molecular orbitals (MOs), $\{\ket{{}^{x}\phi_{i}}\}$, which themselves are formed from a linear combination of $n$ (non-orthogonal) atomic orbitals (AOs), $\{\ket{\chi_{\mu}}\}$, as

[TABLE]

The NOCI eigenstates are identified by solving the generalised eigenvalue problem

[TABLE]

where $H^{wx}=\langle^{w}\Phi|\hat{H}|^{x}\Phi\rangle$ and $S^{wx}=\langle^{w}\Phi|^{x}\Phi\rangle$ are the Hamiltonian and overlap matrix elements in the non-orthogonal basis. Following the approach detailed in Ref. Thom and Head-Gordon, 2009, $H^{wx}$ and $S^{wx}$ are computed by constructing a biorthogonal set of orbitals using Löwdin’s pairing approachLöwdin (1962); Amos and Hall (1961) and then applying the generalised Slater–Condon rules.Mayer (2003) However, for complex-valued MOs we note that the correct form of the unweighted codensity matrices outlined in Ref. Thom and Head-Gordon, 2009 is

[TABLE]

in contrast to Eq. (3) provided therein.

In comparison to orthogonal CI, selecting a relevant set of determinants to form the NOCI basis is less trivial. In particular, the lack of a universal set of MOs across the non-orthogonal determinants removes the intuitive concept of excitation levels (i.e. singles, doubles) that pervades orthogonal approaches. Moreover, the lack of orthogonality between the MOs of multiple non-orthogonal determinants allows excitations of any order to couple in the NOCI expansion. As a result, more Hamiltonian and overlap matrix elements usually need explicit computation, although non-orthogonality can also lead to shorter multi-determinant expansions.

A number of approaches have been proposed for constructing non-orthogonal basis sets, ranging from the use of SCF metadynamics to identify multiple HF solutions,Thom and Head-Gordon (2008) to spin-flipMayhall et al. (2014) (SF-NOCI) and CISOosterbaan et al. (2018, 2019) (NOCIS) inspired approaches. In SF-NOCI, non-orthogonal determinants are constructed by independently relaxing all possible spin-flip excited determinants from a restricted high-spin reference using a frozen active space SCF procedure.Mayhall et al. (2014) In contrast, the NOCIS approach builds the non-orthogonal basis by first individually optimising a set of restricted high-spin reference determinants, each corresponding to the removal of a different core electron, and then reattaching the excited electron to all possible virtual orbitals.Oosterbaan et al. (2018)

Significantly, in both SF-NOCI and NOCIS, the basis states are computed by partially optimising multiple determinants constructed as excitations from a symmetry pure reference. In constrast, by building the NOCI basis from multiple HF solutions identified using SCF metadynamics, one can supplement symmetry-pure determinants with symmetry-broken HF states and capitalise on the electron correlation symmetry-broken states provide at the mean-field level. However, to ensure the existence of symmetry-broken HF states across all geometries — and thus guarantee continuous NOCI energies — HF solutions must be analytically continued beyond the Coulson–Fischer points at which they vanish. To do so, we construct the basis for NOCI using the stationary points of the h-HF equations.Hiscock and Thom (2014); Burton and Thom (2016)

II.2 Holomorphic Hartree–Fock Theory

Using the h-HF approach, real HF solutions can be analytically continued across all geometries by identifying the stationary points of the h-HF energy functionHiscock and Thom (2014)

[TABLE]

Here $\hat{H}$ is the conventional electronic Hamiltonian

[TABLE]

where $V_{\text{N}}$ is the nuclear repulsion, $\hat{h}\quantity(i)$ are the one-electron operators and $r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|$ defines the distance between electrons $i$ and $j$. In constrast to the complex-Hermitian extension of HF,Ostlund (1972) the h-HF equations form a complex-analytic continuation of real HF theory by allowing the MO coefficients to become complex without introducing complex conjugation into the energy function. As a result, we find a constant number of h-HF stationary points to Eq. (5) across all geometries and, where real HF solutions vanish, their h-HF counterparts extend with complex MO coefficients.Hiscock and Thom (2014); Burton et al. (2018) In turn, the h-HF solutions provide a continuous basis for NOCI.Burton and Thom (2016); Burton et al. (2018)

Notably, since Eq. (5) satisfies the Cauchy–Riemann conditions,Fischer and Lieb (2012) it is a complex-analytic function with no concept of minima and maxima. However, stationary points are still well-defined as points where the magnitude of the gradient becomes zero. Identifying optimal h-HF states corresponding to the stationary points of Eq. (5) requires only minor modifications to the conventional SCF procedure.Burton and Thom (2016) In particular, removing the complex-conjugation of the MO coefficients leads to complex-symmetric (cf. Hermitian in conventional complex HF) density, $\mathbf{P}$, and Fock, $\mathbf{F}$, matrices defined as

[TABLE]

and

[TABLE]

where $h_{\mu\nu}$ and $\langle\mu\sigma||\nu\tau\rangle$ are the one-electron and anti-symmetrised two-electron integrals respectively.Attila Szabo and Neil S. Ostlund (1989) Since the optimal molecular orbitals are now given by the eigenvectors of a complex-symmetric matrix, the MO coefficients must form a complex-orthogonal setCraven (1969) (cf. unitary eigenvectors for Hermitian matrices) satisfying

[TABLE]

where $S_{\mu\nu}=\langle\chi_{\mu}|\chi_{\nu}\rangle$ is the (real) AO overlap matrix. Finally, the h-HF energy and orbital energies are complex in general and the Aufbau ordering principleAttila Szabo and Neil S. Ostlund (1989) cannot be applied to select the occupied orbitals on each SCF iteration. Instead, we have found a complex-symmetric analogue of MOM to provide an effective alternative approach.Burton and Thom (2016)

Since all real HF states are also solutions the h-HF equations, from hereon we use the term h-HF to refer only to the stationary points of Eq. (5) with complex orbital coefficients.

III A General Approach

III.1 Active Space SCF Metadynamics

Identifying HF states to form a non-orthogonal basis is an integral, yet challenging, component of NOCI calculations. SCF metadynamics provides an effectively black-box approach for locating multiple stationary points.Thom and Head-Gordon (2008) Starting from an optimised solution, SCF metadynamics generates new initial guess determinants by randomly mixing occupied and virtual MOs, and then optimises these determinants in the presence of a biasing potential to prevent re-converging onto a previously known stationary point. However, even for modest systems, the number of states located may become very large and identifying a suitable set of chemically relevant determinants can be difficult.

Fortunately, in many cases, a dominant subset of “active” MOs can be identified that strongly influence the characteristics of the potential energy surface. The most relevant HF determinants are usually related to different active orbital occupations, along with symmetry-broken states formed from mixing these MOs. We can therefore define an active space SCF metadynamics approach that uses these active orbitals to identify a suitable subset of HF states as follows. Starting from an initial symmetry-pure reference determinant, e.g. the RHF ground state, a metadynamics calculation is run in which only the active MOs are allowed to mix and where the SCF optimisation proceeds only in this active space — i.e. the inactive orbitals remain frozen throughout. This process leads to a set of determinants that differ only the composition of their active MOs, but are not themselves fully optimised HF stationary points. Subsequently relaxing the inactive orbitals by optimising each determinant in the full orbital space then yields true HF solutions that form the basis for NOCI.

Significantly, the number of partially optimised HF states with frozen inactive orbtials is controlled by the size of the active space and is much smaller than the number of states in the full unfrozen HF space. As a result, active space SCF metadynamics provides a more manageable approach to identifying a chemically relevant basis of HF states for NOCI.

III.2 Moving Past the Coulson–Fischer Point

While h-HF theory allows real HF states to be analytically continued into the complex plane, identifying the correct complex direction to follow solutions in the MO coefficient space has previously proved difficult. In particular, since the h-HF energy is symmetric about the real MO coefficient axis, SCF calculations starting from a real guess (for example a real HF solution from a previous geometry) show no preference to move towards any particular complex direction. Additionally, the coalescence of symmetry-broken states often coincides exactly with a real symmetry-pure solution at a cusp catastrophe,Burton et al. (2018) and thus attempts to trace symmetry-broken states into the complex plane often remain stuck on real symmetry-pure solutions instead.

Routinely following h-HF states past the Coulson–Fischer point into the complex requires an understanding of the complex topology of h-HF solutions. Recently, by scaling the electron-electron interaction using the complex parameter $\lambda$, i.e. creating the perturbed Hamiltonian

[TABLE]

we have shown that multiple h-HF states form a continuous interconnected manifold in the complex $\lambda$-plane.Burton et al. (2019) Each individual h-HF state is given by a different branch of a Riemann surface, with Coulson–Fischer points forming isolated exceptional points corresponding to the branch points of the Riemann surface.Burton et al. (2019) Significantly, we can use this continuous topology of h-HF states to follow solutions around a Coulson–Fischer point by tracing a suitable $\lambda$-trajectory in the complex plane.

We illustrate this idea by considering the multiple HF solutions of \ceH2 in a minimal basis set (STO-3G) using two orbitals $\phi_{\alpha}$ and $\phi_{\beta}$ parameterised by the complex angle $\theta$,

[TABLE]

Stationary points correspond to the critical values $\theta_{\text{c}}$. Additionally to the RHF $\upsigma_{\text{g}}^{2}$ state (red solid line in Fig. 1), in the unperturbed case ($\lambda=1$) a doubly degenerate pair of real spatial symmetry-broken UHF (sb-UHF) solutions exist in the dissociation limit (blue solid line in Fig. 1). These correspond to the diradical configurations \ce^$\upharpoonleft$H\bond…H^$\downharpoonright$ and \ce^$\downharpoonleft$H\bond…H^$\upharpoonright$.Coulson and Fischer (1949); Attila Szabo and Neil S. Ostlund (1989); Hiscock and Thom (2014) As the bond length is shortened, the sb-UHF solutions coalesce with the RHF state at the Coulson–Fischer point, before continuing into the complex-$\theta$ plane as h-UHF solutionsHiscock and Thom (2014) (left panel in Fig. 1).

To move smoothly from the real sb-UHF branch, past the Coulson–Fischer point and onto the complex h-UHF branch, we first perturb the real solutions into the complex plane by taking $\lambda=\text{e}^{\text{i}\pi/20}$ (dashed lines in Fig. 1). Although the $\theta_{\text{c}}$ for the RHF state is by fixed spatial symmetry and remains unchanged, the corresponding $\theta_{\text{c}}$ values for the sb-UHF state become complex even for long bond lengths, smoothly connecting the real sb-UHF and complex h-UHF regimes without ever passing through the Coulson–Fischer point. Relaxing these pertubed stationary points back to $\lambda=1$ at each geometry then leads to the unperturbed states required for NOCI.

III.3 Combined Protocol

We are now in a position to define a general protocol that combines active space SCF metadynamics, h-HF and NOCI to compute the ground and excited states of molecular systems. In most systems, the multiple symmetry-broken HF states of interest occur at geometries where strong static correlation is present, for example dissociaton or transition states. Therefore, to ensure every important solution is captured, the initial active space SCF metadynamics calculation is run at all geometries of particular interest. The states identified at each specific structure can then be connected by tracing solutions across the intermediate geometries, including identifying h-HF extensions for any real solutions that coalesce and vanish.

Our combined approach can be summarised as follows:

Identify real HF solutions at the geometries of interest using active space SCF metadynamics. 2. 2.

Perturb states off the real axis using a complex $\lambda$ scaling. 3. 3.

Trace the perturbed solutions across all geometries, identifying any corresponding complex h-HF solutions required. 4. 4.

Relax all states back onto the real axis. 5. 5.

Compute NOCI energies using the resulting basis of multiple h-HF solutions.

We have implemented this approach, comprising active space SCF metadynamics, h-HF and NOCI, in a new dedicated LIBNOCI library available in Q-Chem 5.2.Shao et al. (2015)

IV Results and Discussion

We now illustrate our general approach by considering the ground-state dissociation of \ceF2 and the ground and excited states in the pseudo-Jahn–Teller distortion of cyclobutadiene. Both systems exhibit a challenging electronic structure including multireference character and strong correlation effects. To assess the performance of combining h-HF and NOCI, we compare our results to currently available methods including CIS, CASSCF, EOM-CC and FCI.

IV.1 Computational Details

All h-HF and NOCI energies, along with geometry optimisations (performed at the RHF level) and CIS excitation energies, are calculated using Q-Chem 5.2.Shao et al. (2015) Where h-HF energies become complex, only the real part is plotted. CASSCF energies are computed using the ORCA quantum chemistry packageNeese (2012) and CC calculations, including EOM-CC,Kállay and Gauss (2004) are run using MRCC.111MRCC, a quantum chemical program suite written by M. Kállay, P. R. Nagy, Z. Rolik, D. Mester, G. Samu, J. Csontos, J. Csóka, B. P. Szabó, L. Gyevi-Nagy, I. Ladjánszki, L. Szegedy, B. Ladóczki, K. Petrov, M. Farkas, P. D. Mezei, and B. Hégely. See also Z. Rolik, L. Szegedy, I. Ladjánszki, B. Ladóczki, and M. Kállay, J. Chem. Phys. 139, 094105 (2013), as well as: www.mrcc.hu To obtain approximate FCI comparisons for cyclobutadiene, we use a sCI methodHolmes et al. (2016); Evangelisti et al. (1983); Huron et al. (1973); Giner et al. (2013) shown to provide near FCI accuracy for both ground and excited states.Holmes et al. (2017); Loos et al. (2019, 2018); Chien and Zimmerman (2017); Scemama et al. (2018) In particular, we use the CIPSI (CI using a perturbative selection made iteratively) algorithmEvangelisti et al. (1983); Huron et al. (1973); Giner et al. (2013) implemented in QUANTUM PACKAGE 2.0Garniron et al. (2019) using the frozen core approximation. Futher information on our sCI calculations, including detailed results, is provided in the Supporting Information.

The cc-pVDZ basis set is used throughout, and all energies are provided in units of Hartrees, $\mathrm{E_{h}}$.

IV.2 Dissociation of \ceF2

Due to the combination of strong static and dynamic correlation effects, the ground-state binding curve for \ceF2 is notoriously difficult to compute.Sears et al. (2003); Yang et al. (2013); Gordon and Truhlar (1987); Krylov and Sherrill (2002); Laidig et al. (1987); Krylov (2000); Kowalski and Piecuch (2001) The particularly challenging nature of \ceF2 is present even at the HF level, where RHF vastly overestimates the binding energy (in common with most single bonds)Kowalski and Piecuch (2001) and UHF predicts a completely unbound potential.Gordon and Truhlar (1987); Laidig et al. (1987) Moreover, the RHF solution provides a poor reference wave function at both equilibrium and dissociation geometries,Yang et al. (2013) posing further difficulties for post-HF methods. For example, CCSD overbinds the molecule by almost a factor of two,Kowalski and Piecuch (2001) while the “gold standard” CCSD(T) fails completely at large bond lengths.Kowalski and Piecuch (2001) Even conventional multireference methods struggle to describe the dissociation energy correctly, with full valence CASSCF $(14,8)$ underestimating the potential well depth by around a factor of half.Gordon and Truhlar (1987) In contrast, CCSDT — known to describe single bond breaking wellKowalski and Piecuch (2001) — provides a remarkably close approximation to the exact FCI potential energy surface computed by Bytautas and Ruedenberg,Bytautas et al. (2007) demonstrating the importance of triple excitations for capturing electron correlation.

Taking an active space for SCF metadynamics comprising the valence $3\upsigma_{\text{g}}$ bonding and $3\upsigma_{\text{u}}$ antibonding molecular orbitals (leaving the $\uppi$ orbitals frozen), we identify eight real UHF states in the dissociation limit that directly mirror the minimal basis states of \ceH2.Burton and Thom (2016); Burton et al. (2018) Folowing relaxation in the full orbital space, these states correspond to two spatially symmetry-broken UHF radical states (\ce*↿F\bond…F^$\downharpoonright$ and \ce⇃*F\bond…F^$\upharpoonright$), the bonding $\upsigma_{\text{g}}^{2}$ and antibonding $\upsigma_{\text{u}}^{2}$ RHF states, two non-bonding $\upsigma_{\text{g}}\upsigma_{\text{u}}$ UHF solutions, and two spatially symmetry-broken RHF states resembling the ionic configurations \ceF^+\bond…F^- and \ceF^-\bond…F^+, as shown in Fig. 2. Similarly to the multiple HF states of \ceH2, as the bond length is shortened we find two distinct Coulson–Fischer points involving the coalescence of the radical UHF and ionic RHF solutions with the $\upsigma_{\text{g}}^{2}$ and $\upsigma_{\text{u}}^{2}$ states respectively. For shorter bond lengths, the corresponding h-RHF (dotted red) and h-UHF (dotted blue) solutions continue to exist with complex orbital coefficients.

We consider first the NOCI basis including only the $\upsigma_{\text{g}}^{2}$ RHF ground state and the two radical sb-UHF states along with their h-UHF counterparts, which we denote NOCI(3). Using this minimal NOCI basis yields a binding curve that closely matches the CASSCF $(14,8)$, suggesting that it captures a similar amount of the static correlation as CASSCF $(14,8)$, as shown in Fig. 2. Significantly, the variational flexibility offered by these three determinants appears sufficient to overcome the deficiency of the unbound UHF approximation, leading to a qualitatively correct bound potential while retaining size consistency in dissocation. Like the full valence CASSCF, however, NOCI(3) underestimates the binding energy and overestimates the equilbrium bond length, indicating that NOCI(3) is still unable to capture any dynamic correlation that dominates at shorter geometries.

In contrast, when all eight h-HF solutions are included in the NOCI basis, denoted NOCI(8), the ground state energy is lowered further in comparison to the CASSCF $(14,8)$, particularly in the equilibrium region. The additional multiple Hartree–Fock solutions appear to capture additional dynamic correlation at shorter bond lengths and, although NOCI(8) now overbinds the molecule, the absolute error is reduced in comparison to both NOCI(3) and CCSD. Moreover, the equilibrium bond length predicted using NOCI(8) shows a very promising correspondence to the FCI and CCSDT, as shown in the bottom panel of Fig. 2. Overall we see that, even with such a small CI basis, NOCI provides a good reproduction of the relative binding curve of \ceF2 across the full range of geometries, including where real HF solutions disappear and their h-HF counterparts become complex.

IV.3 Distortion of Cyclobutadiene

Next we consider the cyclobutadiene molecule, which has long been of interest as an archetypal anti-aromatic and highly reactive system.Fratev et al. (1982); Kollmar and Staemmler (1977); Borden et al. (1978); Buenker and Peyerimhoff (1968); Nakamura et al. (1989); Van Voorhis and Head-Gordon (2000); Lyakh et al. (2011); Shen et al. (2008); Demel et al. (2008); Li and Paldus (2009); Deustua et al. (2017); Shen and Piecuch (2012); Musiał et al. (2011); Yang et al. (2013); Rishi et al. (2017); Levchenko and Krylov (2004); Sancho-García et al. (2000); Čársky et al. (1988); Maksić et al. (2000); Eckert-Maksić et al. (2006); Saito et al. (2010); Li and Paldus (2009); Varras and Gritzapis (2018); Voter and Goddard (1986) Correctly describing the ground and excited states of cyclobutadiene has formed the focus of extensive theoretical research.Fratev et al. (1982); Kollmar and Staemmler (1977); Borden et al. (1978); Buenker and Peyerimhoff (1968) In particular methods must balance the multireference nature of the square $\mathcal{D}_{\text{4h}}$ geometry with the single-configurational character of the rectangular $\mathcal{D}_{\text{2h}}$ energetic minimum.

At the $\mathcal{D}_{\text{4h}}$ geometry, the highest occupied molecular orbital (HOMO) is a doubly-degenerate pair of singly occupied orbitals, leading to the $\uppi$-configuration $(\mathrm{a_{2u}})^{2}(\mathrm{e_{g}})^{2}$. This configuration results in a ${}^{1}\mathrm{B_{1g}}$ ground state and a low-lying ${}^{3}\mathrm{A_{2g}}$ first excited state, followed by two higher energy states of symmetry ${}^{1}\mathrm{A_{1g}}$ and ${}^{1}\mathrm{B_{2g}}$.Levchenko and Krylov (2004) As the nuclear coordinates distort towards the rectangular geometry, the molecular symmetry drops from $\mathcal{D}_{\text{4h}}$ to $\mathcal{D}_{\text{2h}}$ and the $\uppi$-configuration becomes $(\mathrm{b_{1u}})^{2}(\mathrm{b_{2g}})^{2}$. The associated descent in symmetry of the ground state (${}^{1}\mathrm{B_{1g}}$ becomes ${}^{1}\mathrm{A_{g}}$) and excited singlet state (${}^{1}\mathrm{B_{2g}}$ becomes ${}^{1}\mathrm{A_{g}}$) leads to a second-order pseudo-Jahn–Teller effect favouring distortion towards the rectangular geometry.Buenker and Peyerimhoff (1968); Nakamura et al. (1989); Borden et al. (1978); Van Voorhis and Head-Gordon (2000)

Significantly, while the ground-state RHF solution provides a suitable reference for the $\mathcal{D}_{\text{2h}}$ geometry, it becomes doubly-degenerate at the $\mathcal{D}_{\text{4h}}$transition state with only one of the degenerate HOMO orbitals being doubly-occupied in each degenerate RHF solution (sketched as RHF 1 in Fig. 4). As a result, single reference methods such as CCSD and CCSD(T) fail to provide even a qualitatively accurate description of the energy surface, with unphysical cusps propagated from the RHF description at the square geometry, as shown in Fig. 6. Removing these cusps requires either the full inclusion of triple excitations (ie. CCSDT),Musiał et al. (2011); Shen and Piecuch (2012); Deustua et al. (2017) or multireference approaches such as multi-configurational SCF,Varras and Gritzapis (2018); Nakamura et al. (1989) generalised valence-bond theory,Voter and Goddard (1986); Čársky et al. (1988) or multireference CC.Lyakh et al. (2011); Shen et al. (2008); Li and Paldus (2009); Sancho-García et al. (2000); Balková and Bartlett (1994)

Starting from the rectangular geometry optimised at the RHF level (see Supporting Information), we replicate the pseudo-Jahn–Teller distortion through the square geometry by simultaneously rotating two opposite corners around the central $\mathcal{C}_{4}$ rotation axis with the distortion angle $\phi_{\text{D}}$, as illustrated in Fig. 3. For each distortion angle $\phi_{\text{D}}$, the \ceC-H and diagonal \ceC-C distances, are held constant. At the square $\mathcal{D}_{\text{4h}}$ transition state geometry we define an active space for SCF metadynamics using the four $\uppi$ orbitals ($\mathrm{a_{2u}}$, $\mathrm{e_{g}}$ and $\mathrm{b_{2u}}$) from the ground state RHF solution and identify twelve low-energy real HF states. After relaxation in the full orbital space, the degeneracies of these states (in order of ascending energy) are 2, 4, 2, 2 and 2. In what follows, we refer to the $n^{\text{th}}$-lowest RHF and UHF states at the square geometry using the notation “RHF $n$” and “UHF $n$” respectively.

Inspecting the $\uppi$ orbitals for each solution — sketched for one state of each degenerate set at the square geometry in Fig. 4 — we observe spatially symmetry-broken orbitals in the UHF 1, UHF 2 and RHF 2 solutions. In contrast, the orbitals of the UHF 3 states preserve spatial symmetry, representing the $(\mathrm{a_{2u}})^{2}(\mathrm{e_{g}})^{2}$ configuration in which both orbitals in the degenerate $\mathrm{e_{g}}$ pair contain one electron each, while the RHF 1 orbitals represent the same valence configuration but with only one of the degenerate $\mathrm{e_{g}}$ orbitals holding both electrons. The degeneracies for each state can be deduced by considering the spatial and spin symmetries of the system. As we move away from the square geometry, the degeneracy of both the RHF 1 and UHF 2 states is broken, splitting into lower (higher) energy states RHF 1a (RHF 1b) and UHF 2a (UHF 2b) with one- and two-fold degeneracies respectively. The single degeneracy of the RHF 1a gtound state leads to the dominant single-reference character at the rectangular geometry. In addition, each of the symmetry-broken solutions coalesces with the symmetry-pure RHF 1a state at a different Coulson–Fischer point, as show in the left panel of Fig. 5. For distortion angles further away from $90\degree$, the h-HF counterparts of the vanishing states continue to exist with complex orbital coefficients (dotted lines in Fig. 5).

Using these HF solutions as a basis for NOCI, we recover continuous and smooth energies that are all variationally lower than their sa-CASSCF $(4,4)$ counterparts (right panel in Fig. 5). Again, NOCI recovers the static correlation required to provide the correct qualitative description of the ground and excited states, while the individually optimised HF basis states capture additional dynamic correlation that quantitively improve the energy relative to sa-CASSCF. Moreover, NOCI yields a smooth relative ground-state autoisomerisation barrier, shown in Fig. 6, with a close agreement to CCSDT. In fact, even a state-specific CASSCF $(4,4)$ ground-state calculation fails to reproduce the relative accuracy of NOCI. We note, however, that the NOCI ground state exhibits an apparent bump at around $\phi_{\text{D}}=99\degree$. This feature results from a complicated “avoided crossing” of the complex h-HF extension of the RHF 2 state with another complex solution in the complex $\lambda$-plane, although a detailed investigation is beyond the scope of the current communication. Crucially, the h-HF states used are consistent in the vicinity of this bump without any coalescence, and thus the NOCI energy remains both smooth and continuous.

Finally, we consider the vertical excitation energies for $\mathcal{D}_{\text{4h}}$ cyclobutadiene. At this geometry, dynamic correlation is known to lower the singlet below the triplet state,Levchenko and Krylov (2004); Kollmar and Staemmler (1977); Nakamura et al. (1989) leading to a violation of Hund’s rules through the dynamic spin-polarisation effect. In Table 1 we compare the excitation energies calculated using NOCI to those computed using CIS, EOM-CCSD, EOM-CCSDT, both state-specific CASSCF $(4,4)$ and sa-CASSCF $(4,4)$, and ex-FCI. In sa-CASSCF $(4,4)$, we simultaneouly optimise the four lowest energy states, while in the state-specific variant we focus on only the lowest energy singlet and triplet states in separate calculations. For the $(^{1}\mathrm{B_{1g}}\rightarrow\ ^{3}\mathrm{A_{2g}})$ transition, NOCI matches the state-specific CASSCF $(4,4)$ result, although both overestimate the excitation energy. In contrast, the uncorrelated CIS leads to an incorrect ordering of the singlet and triplet states, as does the absence of static correlation in EOM-CCSD. For the higher energy singlet-singlet transitions, NOCI gives the closest estimate to the ex-FCI result in comparison to all the methods considered. In contrast, EOM-CCSD significantly overpredicts the excitation energies, while we were unable to converge EOM-CCSDT calculations for these transitions. Overall, NOCI is able to match the performance of state-specific CASSCF $(4,4)$ and outperform EOM-CCSD to provide accurate estimates of multireference excitation energies at a significantly lower cost.

V Concluding Remarks

We have provided a general protocol for using multiple h-HF solutions to construct NOCI wave functions that accurately describe multireference ground and excited energies across all molecular geometries. In particular, we have presented two key algorithmic advances. Firstly, we described an active space SCF metadynamics approach that allows a chemically relevant subset of real HF solutions to be identified. Secondly, to handle any real HF solution that disappears and routinely extend it as a complex h-HF state — and thus ensure smooth NOCI energies — we have presented an approach that exploits the complex topology of multiple h-HF states. We have shown that this combined approach can capture correlation to a similar, if not better, accuracy than CASSCF using an equivalent active space, and can provide accurate multireference excitation energies with remarkably small NOCI expansions. Crucially, this NOCI approach is both systematically improvable and scales very favourably with the system size.

The applications of using fully optimised multiple HF solutions as a basis for NOCI still remain relatively unexplored, and we believe our approach holds great potential for larger systems that are out of reach for conventional multireference techniques. In particular, the inclusion of symmetry-broken solutions that qualitatively represent physical states can provide chemical insight into complicated processes such as electron transfer.Jensen et al. (2018) Moreover, retaining the use of fully optimised HF solutions in the NOCI wave function paves the way for the derivation of molecular gradients and, although NOCI captures mainly static correlation, it also provides a foundation for post-NOCI correlation techniques such as NOCI-MP2.Yost et al. (2013); Yost and Head-Gordon (2016) Ultimately, by combining SCF metadynamics and h-HF theory to create a single systematic and routine approach, we have laid the foundations for using NOCI to describe multireference ground and excited states in general.

Acknowledgements

H.G.A.B. thanks the Cambridge Trust for a studentship and Q-Chem for supporting code development of this work through a summer internship. A.J.W.T. thanks the Royal Society for a University Research Fellowship (UF110161).