Accurate and interpretable representation of correlated electronic structure via Tensor Product Selected CI

Abstract

The task of computing wavefunctions that are accurate, yet simple enough mathematical objects to use for reasoning, has long been a challenge in quantum chemistry. The difficulty in drawing physical conclusions from a wavefunction is often related to the generally large number of configurations with similar weights. In Tensor Product Selected Configuration Interaction (TPSCI), we use a locally correlated tensor product state basis, which has the effect of concentrating the weight of a state onto a smaller number of physically interpretable degrees of freedom. In this paper, we apply TPSCI to a series of three molecular systems ranging in separability, one of which is the first application of TPSCI to an open-shell bimetallic system. For each of these systems, we obtain accurate solutions to large active spaces, and analyze the resulting wavefunctions through a series of different approaches including (i) direct inspection of the TPS basis coefficients, (ii) construction of Bloch effective Hamiltonians, and (iii) computation of cluster correlation functions.

---

1. Introduction

Computational electronic structure theory has developed into an indisputably powerful tool for understanding the quantum mechanical origins of molecular structure and chemical transformations. Progress over the past several decades (in both hardware and algorithmic improvements) has advanced quantum chemistry to the point where the accuracy can often rival that of experiments, particularly for low-energy molecules near equilibrium. However, as the accuracy of a computation increases, so to does the numerical complexity of the solution, making interpretation more challenging.

The need for achieving both quantitative accuracy and qualitative richness was recognized early on, as computers were first becoming increasingly powerful.^1–3 Gaining access to the underlying driving forces of reactions or structure has proven to be one of the most valuable aspects of quantum chemistry. As such, the ability to extract qualitative insight is perhaps more important than simply arriving at a quantitatively accurate result.

Many approaches to extracting conceptual insight from ab initio calculations involves some sort of “localization”. This is because much of our chemical vocabulary is inherently local (oxidation states, bond order, partial charges, hybridization, etc.). The abundance of local chemical constructs is not an accident; molecular structure is generally highly local. For example, the alcohol group in 1-hexanol behaves very similarly to that in 1-heptanol. As such, understanding the local structure of a functional group in one system extends significant reasoning power to other systems. Consequently, much of the effort spent toward extracting qualitative insight involves the localization of orbitals, such as with natural bond orbital methods (NBOs),^4–7 atomic localized molecular orbital methods (ALMOs),^8–10 localization of the density as in atoms in molecules (AIM),¹¹ or even many-electron states^12–19 (though this list is necessarily far from comprehensive). Localization has also been leveraged extensively for reducing computational complexity. Underlying many of these developments is the fact that the density matrix asymptotically approaches linearly scaling for gapped systems in a localized basis.²⁰

All (most) of the methods discussed above ultimately leverage the fact that a Slater determinant wavefunction (or MP2 or CCSD) is invariant with respect to orbital rotations within the occupied or virtual spaces. Orbitals can be mixed to maximize some localizing objective function, and the resulting wavefunctions can then be analyzed in terms of local or non-local contributions.

In contrast, a tensor product space permits a much more explicit notion of locality, one that naturally exposes the ability to factorize into local quantities (entanglement) and allows clear labeling of the entire Hilbert space in terms of unambiguously local quantities. Recently, we have explored the ability to leverage features of tensor product spaces to decrease the computational cost of large active-space calculations.^21–25 In ref. 23 and 24, we introduced a method called Tensor Product Selected Configuration Interaction (TPSCI), which uses a selected CI algorithm to assemble a basis of tensor products of locally correlated wavefunctions, which can provide accurate approximations to full configuration interaction (FCI). In this paper, we demonstrate the ability of these (still expansive) TPSCI wavefunctions to be meaningfully analyzed and interpreted across a rather wide range of physical systems, including non-bonded chromophores, dichromium spin coupling, and the fully delocalized π system of a graphene nano-flake.

2. Theory

We will start by expressing the electronic Hamiltonian in a basis of active orbitals (p, q, r, s), where h_pq: one-electron integrals; assuming that the chosen active space is large enough to capture the necessary physics. This is generally the most limiting assumption in this paper, and work to include the dynamical correlation arising from external orbitals is currently underway in our lab.²⁶

2.1 Orbital clustering

To make progress toward a compact and interpretable representation, we will assume that the active orbitals can be partitioned into disjoint clusters, or groups of orbitals. We will generally use capital letters, I, to index clusters. This orbital partitioning (or “clustering”) is chosen to maximize the interactions within a cluster while minimizing the interactions between clusters. For example, if one had a bimetallic compound (as we consider later in the paper) then one might define each cluster to include all the orbitals centered on a given metal, such that all local dynamical correlation is included as intra-cluster correlation, and weak spin-coupling is considered as an inter-cluster correlation.

Physically, we will assume that the interactions within a cluster are stronger than the interactions between clusters. This is not a formal requirement, but rather one that affects the convergence of the calculations. Each cluster is effectively a new smaller active space, and thus we can construct correlated many-body wavefunctions, |α_I〉, that are completely localized to each local active space (cluster), I. We will refer to these locally correlated states as cluster states, using them to form an orthonormal basis for the full Fock space on each cluster. Likewise, the list of all tensor products of cluster states forms an orthonormal basis for the global Fock space on the full orbital active space. This allows us to represent an arbitrary wavefunction with s states as a linear combination of cluster-state tensor products:

In this representation, the basis vectors can potentially contain a significant amount of electron correlation folded into the local many-body cluster states. This means that the coefficient tensor, c^s_α,β,…,γ, only needs to describe inter-cluster correlation, with all the intra-cluster correlation being folded into the basis vectors. This is essentially the same basis used in the Active-Space Decomposition (ASD) approach of Shiozaki and coworkers.^27–29

2.1.1 Cluster mean-field (cMF) theory. In order to make the most use out of the representation defined above, it is important that the cluster states are defined carefully, so that they incorporate as much relevant electron correlation as possible. Diagonalizing the Hamiltonian projected onto a single cluster (simply keeping only the terms where all creation/annihilation operators act on orbitals within the cluster), yields a set of correlated many-body cluster states that include an exact description of the intra-cluster correlation. Taking a product of the local FCI ground states provides a reasonable approximation for the global ground state, one that becomes exact in the “clusterable” limit.

However, interesting molecular systems generally have non-trivial interactions between clusters, and so this becomes a rather poor approximation in practice. Fortunately, one can easily obtain a much-improved ground-state estimate by including a mean-field description of the inter-cluster interactions when defining the local cluster Hamiltonian, instead of simply projecting out the inter-cluster terms. An approach, called Cluster Mean-Field (cMF) theory, was introduced by Scuseria and coworkers,^30–33 and used by Gagliardi and coworkers under the name variational localized active-space self-consistent field (vLASSCF).³⁴ In this work, we construct correlated cluster states by diagonalizing the local cMF effective Hamiltonian, Ĥ^cMF_I, for each cluster:

where γ^J_rs = 〈0_J|r̂^†ŝ|0_J〉. The cMF Hamiltonian for cluster I depends on the one-particle reduced density matrix (1RDM) of all the other clusters, requiring the cMF solution to be obtained self-consistently.

Bearing a strong resemblance to traditional Hartree–Fock (HF) theory, the self-consistent solution corresponds to the variational minimization of an unentangled (product state) wavefunction ansatz,

the difference from HF being that in cMF the ansatz only enforces the absence of entanglement between clusters. In further analogy to Hartree–Fock theory, a “generalized Brillouin condition” holds, that rigorously uncouples the (converged) cMF wavefunction (eqn (4)) from tensor product states (TPSs) with a single cluster excited:

2.1.2 Orbital optimization. Once converged, the cMF energy is stationary with respect to the local cluster state wavefunction coefficients (local FCI coefficients). This means that the cMF energy is invariant to intra-cluster orbital rotations (assuming each cluster is solved exactly), but variant with respect to inter-cluster orbital rotations. In order to obtain a further improved product state wavefunction, we can make the cMF energy stationary with respect to all orbital rotations. This is essentially a complete active space self-consistent field (CASSCF) calculation with multiple disjoint active spaces.³⁰ While this clearly has the benefit of providing a lower-energy variational solution, perhaps more importantly, it removes most of the arbitrariness of the orbital clustering. Assuming each cluster has a Hilbert space dimension greater than one, orbital optimization with a single tensor product state wavefunction will naturally tend to localize the orbitals, so as to maximize electron correlation. Consequently, methods that use the cMF wavefunction as a reference state will be well-defined, not dependent on a particular heuristic for orbital localization.^22,31

2.2 Tensor product selected CI (TPSCI)

While cMF provides a qualitatively attractive approximation for the ground state of a clustered molecular system, quantitative accuracy is clearly missing due to the neglect of all inter-cluster entanglement. Analogous to the common approach of using substituted Slater determinants for a basis, we will use substituted tensor product states as a basis for the full Hilbert space, where each TPS is typically taken to be an eigenstate of a cluster's cMF Hamiltonian. This is similar to a CI analogue of the Block-Correlated Coupled Cluster (BCCC) approach of Li.³⁵

The Hartree–Fock-based Slater determinant basis and the TPS basis are equivalent, in that both span the full space. However, as soon as truncations are made, the two bases span different spaces. One benefit of working in a TPS basis where the local Hamiltonians are diagonal is that, due to the fact that local correlation is folded into the basis states themselves, the low-energy solutions become more heavily concentrated on a smaller number of basis states. Consequently, computational methods that exploit sparsity (e.g., selected CI) might be expected to be more performant in the correlated TPS basis than in a Slater determinant basis. In ref. 23 and 24 we demonstrated that this is often true, and can sometimes be leveraged for computational benefit.

In the Tensor Product Selected CI (TPSCI) method,²³ we use the general CIPSI³⁶ algorithm to discover and exploit, in a bottom-up fashion, the sparsity of the exact wavefunction in a TPS basis. This uses perturbation theory to iteratively discover the non-negligible TPSs that are needed to accurately approximate the exact solution. This is done via the following steps:

(1) Diagonalize Ĥ in the current variational space (this being a list of TPS basis states that are expected to have large amplitude in the exact solution).

(2) Apply the Ĥ to the current variational space eigenvector. This couples the variational space to the external space.

(3) Compute the first-order wavefunction in the external space.

(4) Move the external configurations with large first-order coefficients from the external space to the variational space.

(5) If the variational dimension increases, go back to step 1. If not, exit.

This overall iterative loop is shown in Fig. 1. In principle, all local FCI cluster states would be computed and used to form the basis for state space. While this is tractable for small clusters, for larger clusters, this becomes computationally prohibitive, and high-energy states are generally discarded prior to running TPSCI. We generally use M to refer to the maximum number of cluster states kept in a particular Fock sector of a cluster. For each particle number subspace included, the corresponding lowest M eigenstates are computed. Then the Ŝ⁺ and Ŝ⁻ operators are applied to those cluster states to generate the basis for the higher M_s sectors. More details can be found in ref. 24.

Fig. 1 Schematic illustration of the selected CI algorithm used in TPSCI to build a basis of tensor product states. This is iterated until the dimension of the variational space stops growing.

The computational limitations of conventional (determinant basis) CIPSI are generally determined by the size of the dimension of the variational space. Because TPSCI uses correlated TPS basis vectors, more correlation energy is typically recovered with smaller variational spaces, addressing the most significant bottleneck. This comes at a cost, however, during the matrix-element evaluation. Whereas evaluating matrix elements in the Slater determinant basis is extremely efficient, matrix elements in the TPS basis are significantly more expensive. For a 3-cluster example, consider the following Hamiltonian contribution that contains three operators (p̂^†q̂^†r) on cluster 1, and one operator (s) on cluster 2:

Computing the matrix element of this particular Hamiltonian contribution between two arbitrary TPS basis vectors will require the contraction of an integral sub-block with tensors of local quantities:where χ is a sign determined by the number of electrons in state |α₁〉, and the Γ tensors are the precomputed local operator matrices in the cluster basis, e.g.:Because the matrix elements require a series of tensor contractions, instead of just a single access from an array, the construction of matrix elements becomes the key bottleneck in TPSCI. However, in ref. 23, we compared TPSCI to Heat Bath CI,³⁷ and found that we were able to obtain significantly lower variational energies using TPSCI than with Heat Bath CI. Once the TPSCI variational space has converged, we have also found that it is sometimes beneficial (especially for ground-state problems) to perform a higher-order singular value decomposition (HOSVD) of the resulting wavefunction tensor, to rotate the cluster states into a form that diagonalizes the local cluster reduced density matrices (within subspaces that preserve local particle number and Ŝ^z). More details about the implementation and matrix element construction can be found in ref. 23 and 24.

2.2.1 Bloch effective Hamiltonian. Once an accurate TPSCI wavefunction is obtained, one is often interested in more than just the associated energy. Being able to extract qualitative information to aid in communicating and reasoning about the underlying electronic structure is extremely valuable, and TPS wavefunctions are uniquely interpretable. Since the basis states in TPSCI are essentially diabatic states, it is a very natural extension to use the concept of Bloch effective Hamiltonians to extract quantitative relationships between qualitatively meaningful degrees of freedom^38–47 for analysis.

We start by defining a “model space”, {|ϕ_i〉}, which is taken to be the set of physically meaningful TPSs that qualitatively define the structure or process. To ensure that the model space is actually relevant to the physics computed, the exact low-energy states of the system, |Ψ_s〉, should have relatively large projections onto the model space, i.e.,

where, . Next, we seek a Hermitian effective Hamiltonian which exists only in the model space, but that yields the exact energy spectrum. While this is often done in a bottom-up fashion through approaches like quasi-degenerate perturbation theories or Schrieffer–Wolff transformations, if one already has access to the exact target eigenstates, an effective Hamiltonian (specifically, a Bloch effective Hamiltonian) can be obtained simply in a top-down fashion by direct projection,

Ĥ^Bloch = |[small psi, Greek, tilde]_s〉E_s〈[small psi, Greek, tilde]_s| (10)

where,will always exist when eqn (9) holds. The individual matrix elements of Ĥ^Bloch then contain quantitative relationships between qualitatively meaningful states.

2.2.2 Cluster correlation functions. In addition to the Bloch effective Hamiltonian, which gives us a state-universal description of the interactions between physically intuitive degrees of freedom, we also often want to characterize specific states in terms of physically intuitive variables.

Following the recent work of Luzanov, Krylov, and Casanova,^48,49 the local TPS representation makes it simple to compute cluster correlation functions of various local cluster operators to characterize states in terms of observables. A two-cluster correlation function for operator Ô is given as the covariance between the operator localized onto each individual cluster:

cov(Ô_I,Ô_J) = 〈Ψ_s|Ô_IÔ_J|Ψ_s〉 − 〈Ψ_s|Ô_I|Ψ_s〉〈Ψ_s|Ô_J|Ψ_s〉 (13)

where the covariance of an operator with itself is just the variance, which will be used to measure the local fluctuations in a given cluster. Depending on the system, we will consider correlations between the following cluster operators: local charge (particle number), N̂_I, local spin projection, Ŝ^z_I, local spin Ŝ²_I, and local excitation, Q̂_I = 1̂ − |0_I〉〈0_I|, which indicates that cluster I is excited out of its cMF ground cluster state.

3. Results

In the following sections, we explore the ability to simultaneously obtain quantitative yet interpretable approximations to large active spaces. Because our representation is ideal for separable clusters, to obtain insight into the transferability of the formalism, we explore systems that span a broad spectrum in terms of clusterability, going from a non-bonded tetracene tetramer (Sec. 3 3.1), to an anti-ferromagnetically coupled dichromium complex (Sec. 3 3.2), to a completely delocalized graphene-flake model, hexabenzocoronene (Sec. 3 3.3).

We use PySCF software for performing any necessary geometry optimizations, Hartree–Fock calculations, and integral generation.⁵⁰ All TPSCI calculations are performed using our open-source FermiCG software.⁵¹

3.1 Singlet-fission: tetracene tetramer

Singlet fission is the photophysical process by which a bright singly excited state, |S₁〉, is converted into two lower-energy triplet states, |T₁〉 + |T₁〉, by way of a multiexcitonic intermediate state, |¹TT〉.⁵² Because this process converts a single photon into two excitons (each of which can split into charge carriers), materials that exhibit singlet fission have promising applications in solar cells due to the possibility of overcoming the Shockley–Queisser limit for efficiency.⁵³

A dimer model description of the singlet-fission process requires a total of 8 states: the ground state, |S₀S₀〉, two local bright states, |S₁S₀〉 and |S₀S₁〉, two local triplet states, |S₀T₁〉 and |T₁S₀〉, and the three biexcitonic states arising from the product of two triplet states, |T₁T₁〉. The biexcitonic state is characterized as an entangled pair of local triplets, which can be spin-coupled into either a singlet |¹TT〉, triplet |³TT〉, or quintet |⁵TT〉, represented in the diabatic basis via their Clebsch–Gordan coefficients,

where the first triplet in each state refers to chromophore A and the second triplet to chromophore B. Although |¹TT〉 is the main intermediate, as it is spin-allowed, it has been shown that the triplet and quintet states can play a role in the separation process.⁵⁴

While a dimer model captures the key intermediates, it is too small to describe additional physical effects that occur in bulk systems. For example, the initial bright state is generally understood to delocalize over several monomers, increasing the number of localized biexcitons to which it can couple. Further, it has been seen that singlet fission rates are increased by the involvement of a non-nearest-neighbor biexciton |¹T⋯T〉, which is absent from the dimer model by construction.^55,56 Recent work has further demonstrated the importance of beyond-dimer effects.^57,58

Because methods that rely on single or even double excitations struggle to accurately describe two-electron excitations, multireference methods, such as CASPT2, are often required to capture the wide range of electronic character. While this is suitable for a couple chromophores, active-space methods typically grow exponentially with the number of chromophores, making it difficult to extend to larger systems.

TPSCI (similar to ASD, which preceded it^27–29) is well-suited for treating collections of chromophores because the physical system efficiently maps onto the tensor product basis. In a recent paper,²⁴ we demonstrated that TPSCI was able to provide accurate approximations to complete active space configuration interaction (CASCI) for a large (40e, 40o) active space, which incorporated 10 active orbitals for each of the four tetracene chromophores. In this paper, we explore this further, using the TPS structure to facilitate further analysis and characterization of the resulting wavefunctions.

In this subsection, we use TPSCI to go beyond the minimal dimer example for singlet fission and explore dressed Hamiltonains and local cluster operator correlations for a tetracene tetramer (shown in Fig. 2(a)) using a large (40e, 40o) active space.

Fig. 2 Bare and Bloch effective Hamiltonians for a tetracene tetramer with diagonal entries subtracted to show off-diagonal couplings in meV. Active space of (40e, 40o) with system and clusters labeled. (a) Tetracene tetramer with the associated clusters labeled; (b) Hamiltonian matrix for the model space (diabatic basis); (c) Bloch effective Hamiltonian obtained by projecting TPSCI wavefunctions onto the model space. The non-singlet states have been removed for clarity.

3.1.1 Active-space selection and clustering. In order to construct an active space that contains the relevant orbitals for describing both the local S₁ and T₁ states, we used the configuration interaction singles-natural orbitals (CIS-NO)⁵⁹ approach to select our active space in the 6-31G*⁶⁰ basis. We first performed a configuration interaction singles (CIS) calculation for the first singlet and triplet on each chromophore and averaged these states into the one-particle reduced density matrix (1RDM). By diagonalizing the 1RDM, we obtained a set of natural orbitals from which we extracted the 40 most correlated orbitals (i.e., those that have the most fractional occupancy) as our active space. We then localized these 40 orbitals using the Pipek–Mezey method⁶¹ and then grouped the orbitals into four (10e, 10o) clusters on each chromophore for an overall active space of 40 orbitals and 40 electrons (40e, 40o). We are currently developing approximate solvers (such as restricted active space configuration interaction (RAS-CI)) for obtaining the local cluster states. This will allow us to consider clusters larger than the relatively small ten-orbital clusters used here.

3.1.2 Bloch effective Hamiltonian. To analyze our TPSCI results, we start by computing a Bloch effective Hamiltonian by projecting the TPSCI eigenvectors onto our diabatic basis (i.e., model space). The diabatic basis for a tetracene tetramer includes the biexciton diabatic states for each pair in addition to the singly excited states on each chromophore. In the tetramer, there are six possible dimer configurations and each generates three M_s = 0 spin components, |T⁺T⁻〉, |T⁰T⁰〉, and |T⁻T⁺〉, which leads to a total of 18 diabatic biexciton states. We also observe in the tetracene monomer that both T₁ and T₂ are lower in energy than the first singlet excited state, S₁; thus all three of these states must be represented for each chromophore in our model space. In total, our model space includes 31 diabatic states. However, to simplify the picture, we focus on the singlet model space, where the biexcitonic states have been mixed using their Clebsch–Gordon coefficients to form proper |¹TT〉 diabatic states. This reduces our model space from 31 states to 11 states.

In Fig. 2(b) and (c), we plot the bare Hamiltonian and effective Hamiltonian in the model space as described above, with the columns arranged to correspond with the cluster labels in Fig. 2(a). The diagonal energies are subtracted to better reveal the off-diagonal activity in meV. In both plots, the Hamiltonian is blocked by state type: the ground state, four singlet excited states, and six singlet biexcitons, where the non-singlet states have been omitted for clarity. As expected, the singlet excitons couple very strongly to each other, both in the bare and effective Hamiltonians, which ultimately gives rise to bright-state delocalization. There is only negligible coupling between the S₁ and biexcitons in the bare plot, but after including external space correlations, we see a significant growth in the strength of the effective coupling. These are listed explicitly in Table 1. Whereas the |¹TT〉 states on the 3 herringbone dimers develop significant S₁ coupling after the inclusion of the external space, the planar dimers remain uncoupled from the bright spectrum. In addition to strengthening the coupling between the bright states and the biexcitonic states, the inclusion of higher-energy states also induces couplings between the biexcitons themselves.

Table 1 TPSCI Effective Hamiltonian to show coupling strengths between S₁ and |¹TT〉 in meV

Ĥ^eff		|^1TT〉
		1,2	1,4	2,3	1,3	2,4	3,4
〈S1|	1	−6.43	19.46	0.33	0.80	1.31	0.01
	2	3.66	2.47	10.16	1.04	0.50	0.00
	3	0.37	−0.82	14.36	−0.64	−0.08	−0.07
	4	−0.23	−6.62	0.42	0.10	0.37	0.08

---

3.1.3 Correlation analysis. As mentioned in Section 2.2.2, correlation functions of various local cluster operators can be used to characterize the adiabatic TPSCI wavefunction in terms of physically meaningful relationships between clusters. In Table 2, we list the expectation values of the cluster excitation operators, which measure the amount of excited-state character in each state, and by summing over clusters, the total excitation rank of each excited state. This allows us to unambiguously identify the biexcitons, which we use to label the states accordingly in Fig. 3, where we compute the inter-cluster cumulants for cluster particle numbers (N̂_I), cluster spin projections (Ŝ^z_I), and cluster excitations (Q̂_I) for each of the six singlet biexcitons ¹(TT) and four bright states.

Table 2 Tetracene tetramer local excitation strengths. Columns correspond to the expectation values of the local cluster excitation operator, Q̂_I = 1̂ − |0_I〉〈0_I|. A value of 0 indicates that the cluster is in its local ground state. A value of 1 indicates that the cluster is always in a local excited state. All the local excitation ranks of the states are summed

State	〈Q̂1〉	〈Q̂2〉	〈Q̂3〉	〈Q̂4〉	Excitation rank
6	0.90	0.89	0.20	0.04	2.01
7	0.63	0.54	0.81	0.03	2.02
10	0.95	0.10	0.01	0.95	2.01
12	0.45	0.54	1.00	0.02	2.01
19	0.05	0.11	0.91	0.95	2.02
21	0.05	0.86	0.11	0.99	2.01
24	0.63	0.27	0.13	0.11	1.14
25	0.24	0.57	0.08	0.16	1.05
28	0.18	0.08	0.81	0.02	1.09
31	0.06	0.20	0.01	0.80	1.07

---

Fig. 3 Tetracene tetramer correlation functions. (top row) Particle number, N̂_I. (middle row) Spin projection, Ŝ^z_I. (bottom row) Cluster excitation, Q̂_I. The first 6 plots from the left are for the |¹TT〉 states. The last 4 plots from the left are for the |S₁〉 states. Each matrix column/row corresponds to cluster 1 to 4, as labeled in Fig. 2. The color scale for each correlation function is shown on the far right.

Looking first at the particle number correlations, N̂_I, we see that overall, the “dark” |TT〉 states are relatively quiet compared to the charge correlations present in the bright states. This is entirely expected based on the physical characteristics of the bright vs. dark states. Bright states have relatively large amounts of charge-transfer character mixed in. The presence of charge transfer makes the local particle number less well defined, which increases a cluster's charge variance, and similarly increases the anti-correlation between two clusters' charge states (when the donor is cationic, the acceptor has a high probability of being anionic). Although weaker than the lowest S₁ state, we see clear signatures of charge correlation present in a couple of the biexcitons (state 6 and 10). This is consistent with an analysis of the wavefunction itself. If we compute the amount of charge transfer present in each state, we find that out of all of the biexcitons, states 6 and 10 have the highest percentage of CT character, 4.6% and 4.5%, respectively (see the ESI† for all state CT compositions). Charge correlations between clusters entangled into a biexciton also indicate significant superexchange, which stabilizes the low-spin biexciton.^62,63

Considering next the Ŝ^z_I correlations, we see the opposite trend, where the biexciton states have significantly more pronounced correlations, and the bright states are featureless (as would be expected from the lack of local triplet character). This is also consistent with the nature of the different sets of states. A |¹TT〉 biexciton is characterized as two entangled triplet states coupled into a singlet state. Because the total M_s is zero, when the first monomer is in the M_s = 1 microstate, the entangled partner is very likely to be found in the M_s = −1 microstate. This entanglement leads to a very strong Ŝ^z_I covariance.

Using the spin correlation as a way to label the biexcitons,⁴⁸ we can identify that state 6 is a (1,2) biexciton (meaning that it primarily exists on chromophores 1 and 2), states 7 and 12 are superpositions of (1,3) and (2,3) biexcitons, and state 21 is a (2,4) biexciton. Looking at state 19, we see what resembles a (3,4) biexciton, although the overall magnitude is much smaller. In order to understand this, we can look at the Ŝ² expectation value of the state (shown in the ESI†). We find that in this case, the (what we labeled to be) |¹TT〉 state has significant spin contamination of about 1.1. This is a consequence of the fact that the |¹TT〉 and |³TT〉 states are approximately degenerate, meaning that any arbitrary mixing of the two states is also an eigenstate. This mixing of the two spin states essentially creates a “broken-symmetry” state, where one chromophore is M_s = 1 and the other is M_s = −1. By locking the local spin vectors, the local Ŝ^z_I fluctuations are diminished, and hence the ability to have significant covariance with any other cluster. This could be corrected by tightening our TPSCI convergence, or by rediagonalizing Ŝ² in nearly-degenerate subspaces.

Looking more closely at state 12, we see that clusters 1, 2, and 3 all have significant Ŝ^z_I fluctuations, but only the 1,3 and 2,3 pairs are spin correlated. Considering the cluster excitation Q̂_I covariance plot, we see that while cluster three has zero fluctuations (it is consistently excited), clusters 1 and 2 are strongly correlated. This means that cluster 3 is always excited, but when cluster 1 is excited, cluster 2 is most likely in the ground state, and vice versa. This suggests a situation where a triplet state on cluster 3 forms a biexciton with a triplet delocalized between clusters 1 and 2. This is consistent with an analysis of the individual TPS state coefficients, were we find that 98% of the wavefunction is characterized as a superposition of the 1,3 and 2,3 biexcitons, ψ₁₂ ≈ 0.67¹(T₁S₀T₁S₀) + 0.73t¹(S₀T₁T₁S₀).

3.2 Cr₂-complex effective Hamiltonian

Multi-center organometallic complexes often exhibit interesting physics, such as single-molecule magnetism,^64,65 or valuable catalytic capabilities, such as water oxidation,⁶⁶ or nitrogen fixation.⁶⁷ While having a computational method that could efficiently compute the low-energy structure of organometallic compounds would be highly valuable, several physical features of these systems make this challenging. When multiple weakly interacting metal centers possess unpaired electrons, the resulting low-energy states are highly multiconfigurational, making conventional approaches like perturbation theory or coupled cluster theory inappropriate, as they require a qualitatively correct single determinant wavefunction as a reference.

Because the product of multiple high-spin centers leads to a large number of spin states, it is not always possible to predict, a priori, the spin multiplicity of a multi-center organometallic complex. Consider a simple biradical Hubbard model, represented in a local minimal orbital basis, as shown in Fig. 4(a). Here, the low-energy configurations are both open-shell broken-symmetry states, and proper eigenstates must be superpositions of these two configurations, one being a singlet, and the other the M_s = 0 component of the triplet. Because the ionic (or charge transfer) configurations are both singlets, they can only mix with the singlet combination of neutral configurations, which is ultimately the origin of the antiferromagnetic coupling in biradicals.

Fig. 4 Cartoon illustration of the electronic structure of an biradical approximately mapped onto an isotropic Heisenberg model at the large U limit, both restricted to the relevant M_s = 0 subspace. (a) The 4 possible Slater determinants. Assuming a localized orbital basis, the bottom two correspond to neutral excitations and the top two are ionic or charge-transfer configurations. (b) The 2 possible spin configurations after being mapped onto a Heisenberg model with quasi-degenerate perturbation theory (QDPT).

This picture can often be simplified significantly. For systems where the electron repulsion is sufficiently large, such that hopping of an electron from one center to another incurs a high energy barrier (i.e., U ≫ t), the influence of charge fluctuations can be approximately downfolded using quasi-degenerate perturbation theory into an effective spin Hamiltonian, called the Heisenberg–Dirac–van Vleck Hamiltonian:⁶⁸

Ĥ^HDvV = −2JŜ₁Ŝ₂, (17)

where . For this model, the low-energy spectrum is completely determined by the value of J.

Starting from the ab initio Hamiltonian instead of the Hubbard model, one finds that the zeroth-order term also contains the non-local direct exchange integral. Since bare exchange stabilizes the high-spin states, and the second-order super-exchange term stabilizes the low-spin states, even getting the sign correct for the exchange coupling constant can be difficult, as the value of J is determined by the subtle interactions coupling the metals with each other and with the ligands. However, once known, the relative ordering of the spin states can be directly written down in terms of J. While this is indeed an approximate description of the low-energy electronic structure, it is profoundly useful, as the prediction of J is also often what connects experiment to theory, where J is commonly fitted to experimental magnetic susceptibility measurements.

The most common approach to computing J from ab initio quantum chemistry is to use DFT, where one of the degenerate broken-symmetry configurations is optimized, followed by a spin projection formula, originally proposed by Noodleman,⁶⁹ and then improved by Yamaguchi.^70,71 While this approach has been widely used due to its conceptual simplicity and computational efficiency, there are downsides. First, the formalism doesn't actually ever compute the proper low-spin wavefunction, and so only the energy is generally able to be extracted. Second, the results become highly functional-dependent. While all systems demonstrate some density-functional dependence, spin-coupled complexes are intrinsically more sensitive because the percentage of exact exchange directly affects the relative energies of the high-spin and broken-symmetry states.^72,73 Finally, DFT doesn't offer a path toward systematic improvements, making it difficult to compare results.

Because of these reasons, multireference methods like CASSCF and CASPT2 are often used to model exchange-coupled systems. However, the associated computational cost limits the active-space size, making it difficult to converge results to the point where quantitative comparison to experiment is possible. Often one finds that dynamical correlation (involving interactions with non-magnetic orbitals) has a significant impact on the value of J, generally strengthening the antiferromagnetic interactions. As such, density matrix renormalization group algorithm (DMRG) has emerged as the standard benchmark method for computing exchange coupling constants in organometallic compounds,^74–83 although if only the value of J is needed, efficient approaches that combine spin-flip methods have also been useful.^44,45,47,84

In this subsection, we present the first application of TPSCI to a transition-metal compound, a tris-hydroxy-bridged Cr(III) dimer, [L₂Cr(III)₂(μ-OH)₃]⁺³ (L = N,N′,N′′-trimethyl-1,4,7-triazacyclononane (Fig. 5(c)), which has a J value that was experimentally fitted to a value of −66 cm⁻¹.⁸⁵ Recently, Pantazis studied this system using DMRG to solve the low-energy states in a large orbital active space of up to (30e, 22o),⁷⁴ calculating an exchange coupling constant of −23.9 cm⁻¹.

Fig. 5 Convergence and extrapolation of Cr₂ low-energy spectra for the def2-SVP basis and a (32e, 38o) active space. (a) Plot of the TPSCI variational energy as a function of the computed PT2 correction. The solid line is the linear fit. Units in milliHartree. M = 100. (b) Clustering of the 38-orbital active space. (c) Molecular structure.

3.2.1 Active-space selection. For our calculations, we extend the size of the active space by including the orbitals that overlap most strongly with the 3d and 4d orbitals on each Cr center, as well as the 2p and 3p oxygen orbitals on each bridging OH⁻¹ ligand, leading to an overall orbital active space of 38 orbitals. Our active space was obtained by first optimizing the restricted open-shell Hartree Fock (ROHF) wavefunction for the high-spin heptet state. We then define a set of atomic orbitals for which we would like to span as closely as possible without destroying the ROHF reference. For this system, we take as our projection space (μ_A) the Cr 3d and 4d atomic orbitals, and the bridging O 2p and 3p orbitals, leading to a total of 38 orbitals. We then separately project the doubly occupied (i), singly occupied (s), and virtual (a) orbital spaces onto this atomic orbitals (AO) subspace, providing matrices C_μA,i,C_μA,s,C_μA,a. We then perform separate SVD's on each projected subspace, and keep the largest singular vectors from each orbital space. As such, we start with 38 atomic orbitals and end up with 38 molecular orbitals. Following this automated procedure produced an active space consisting of 13 doubly occupied, 6 singly occupied (the full ROHF open-shell space), and 19 virtual orbitals, leading to a (32e, 38o) active space. While this is not the only way to yield an active space, it was convenient for our purposes, as the resulting orbital active spaces are already localized to our target system. In the future, more extensive tests will be performed for automating the construction of localized active spaces. The active orbitals are shown in the ESI.†

3.2.2 Clustering. The 38 active orbitals described above were then organized into 5 clusters, as depicted in Fig. 5(b). Here, each Cr atom cluster defined a local (7e, 10o) active space, and each oxygen a (6e, 6o) local active space. For each cluster, a local many-body basis was constructed from the M lowest-energy states for each sector of Fock space, which contained up to N_I ± δ_elec number of electrons, where we set δ_elec = 3 for these calculations. For example, each Cr cluster has a basis of up to M states obtained by diagonalizing the CMF Hamiltonian for each of the following active spaces: (4e,10o), (5e,10o), (6e,10o), (7e,10o), (8e,10o), (9e,10o), and (10e,10o). Similarly, each bridging OH⁻¹ ligand had 7 different active spaces centered at (6e, 6o).⁸⁶

As is commonly done with selected CI approaches, significantly improved approximations to the energy can obtained by performing a series of selected CI calculations with varying thresholds and extrapolating to the zero-error limit, which is taken as an approximation to full CI. While more sophisticated extrapolation schemes have been proposed,⁸⁷ we use the common approach of assuming a linear relationship between the variational energies and the PT2 correction.⁸⁸ In Fig. 5(a), we plot the variational TPSCI energy of the 4 lowest eigenstates as a function of the PT2 correction to each state. By extrapolating this linear relationship to zero PT2 correction, we can obtain an estimate of the exact eigenvalues of the Hilbert space defined by M. In Fig. 5(a), we show the extrapolation for the M = 100 calculations.

In order to compute J, we can use the Landé interval rule derived from the energy spectrum of a two-site Heisenberg model, J = (E(S − 1) − E(S))/2S. After computing the lowest energy S = 0, 1, 2, and 3 spin states, we could use any of the gaps to compute J. If the ab initio system was to be perfectly described by the Heisenberg model, then the computed J value would be independent of the particular energy gaps we were to choose. However, the Heisenberg model is rarely exact, and so we can partially quantify how approximate the model is by comparing J values computed with different energy gaps. We list the various J values in Table 3, using either the best variational TPSCI energies, the TPSCI + PT2 corrected energies, or the extrapolated energies. Results for both M = 100 and M = 200 are included.⁸⁹

Table 3 Exchange coupling constants (cm⁻¹) for the Cr₂ compound with the def2-SVP basis and (32e, 38o) active space. J refers to H = −2JŜ₁·Ŝ₂. J(S0,S1) denotes which spin states are used to compute J via the Landé rule. “TPSCI” refers to the best variational energy obtained, using ε_CIPSI = 2e − 4. “TPSCI + PT2” is the best variational energy plus the state-specific PT2 correction. “Extrapolated” uses differences between the extrapolated energies. M is the maximum number of cluster states computed for each cluster Fock sector. The dimension of the sparse variational TPSCI subspace is 97 357 for M = 100 and 127 493 for M = 200

	J (S0,S1)	J (S1,S2)	J (S2,S3)
TPSCI
M = 100	−25.4	−25.7	−26.6
M = 200	−26.6	−27.0	−27.7
TPSCI + PT2
M = 100	−26.7	−27.3	−28.3
M = 200	−28.3	−28.9	−29.9
Extrapolated
M = 100	−28.0	−28.7	−30.0
M = 200	−29.3	−30.3	−31.3

---

Inspecting first the effect of energy extrapolation, we find that the extrapolated J values are larger than the TPSCI + PT2 values by only around 1 cm⁻¹, and that doubling the size of M from 100 to 200 only increases the J value by another 1 cm⁻¹, despite the fact that this also increases the dimension of the total accessible Hilbert space significantly from 2.7 × 10¹⁴ to 6.1 × 10¹⁵. Compared to the recent (30e, 22o) DMRG-SCF calculations, which yielded a J value of −23.9 cm⁻¹,⁷⁴ our computed values are slightly larger, in good agreement with the reported CASSCF(6e,10o)-NEVPT2 value of −31.8 cm⁻¹.⁷⁴

3.2.3 Bloch effective Hamiltonian. In order to further analyze the results listed in Table 3, we compute a Bloch effective Hamiltonian, providing access to the individual effective (dressed) interactions between the various spin microstates that lead to the low-energy spectrum. A qualitative description of this complex assigns each Cr center a well-defined oxidation state (III) and spin state . As such, the low-energy spectrum is expected to be dominated by the 16 spin configurations that form a basis for the S = 3, 2, 1, and 0 states, providing a clear definition for our model space. However, because the Hamiltonian preserves spin, we can restrict our focus to only the global M_s = 0 subspace, and take our model space to be the corresponding 4-dimensional subspace. In Fig. 6, we plot both the bare Hamiltonian (Fig. 6(a)) and the Bloch-effective Hamiltonian (Fig. 6(d)) in the model space. The corresponding low-energy spectra are provided in Fig. 6(b) and (c), respectively.

Fig. 6 Bare and Bloch effective Hamiltonians for the Cr₂(III) complex. (32e, 38o) active space. Units of J in cm⁻¹. (a) Hamiltonian matrix in the basis of M_s = 0 tensor products of local spin states, where all the bridging ligands are in singlet states, and the Cr centers are in (b) Energy spectrum of bare Hamiltonian in a local basis. (c) Energy spectrum of extrapolated TPSCI results. (e) Bloch effective Hamiltonian obtained by projecting TPSCI wavefunctions onto a local basis.

There are two main features of the effective Hamiltonian that emerge from the implicit inclusion of the external space correlation: the spin-coupling interactions have their signs flipped, and their magnitudes are increased. The result of this is that the system changes from ferromagnetic to antiferromagnetic coupling and the gaps between the spin states are approximately doubled.

This qualitative result is consistent with a recent study (ref. 90) where J values for a similar dichromium(III) complex were computed using the vLASSCF-SI method. Because the vLASSCF method works in a similar TPS basis, they were able to evaluate the impact of explicitly including some charge-transfer configurations. They, too, found that this effected a qualitative change in the sign of J, switching from ferromagnetic to anti-ferromagnetic.

In order to further analyze this Cr₂(III) system, we can compute cumulants of local observables, as mentioned above.^48,49 In Table 4, we list the expectation values, the variances, and the covariances of a few operators local to the Cr₂(III) centers, including the local particle number, Ŝ^z_Cr, and Ŝ²_Cr. We also include the global Ŝ² because our basis is not spin-adapted, so there is potential for spin-contamination, although our calculations are converged tightly enough to reduce spin-contamination to the 4th decimal place.

Table 4 Local expectation values for the Cr₂(III) complex. Variance is computed as: var(Ô_Cr) = 〈Ô²_Cr〉 − 〈Ô_Cr〉². Covariance is computed as: cov(Ô_Cr) = 〈Ô_CrAÔ_CrB〉 − 〈Ô_CrA〉〈Ô_CrB〉

Root	1	2	3	4
〈Ŝ^2〉	0.000	2.000	6.000	12.000
〈N̂Cr〉	7.010	7.010	7.010	7.010
var(N̂Cr)	0.019	0.018	0.016	0.014
cov(N̂Cr,N̂Cr)	−0.004	−0.004	−0.002	0.000
〈Ŝ^zCr〉	0.019	0.013	0.003	0.003
var(Ŝ^zCr)	1.248	2.042	1.246	0.453
cov(Ŝ^zCr,Ŝ^zCr)	−1.245	−2.039	−1.242	−0.449
〈Ŝ^2Cr〉	3.741	3.742	3.745	3.750
var(Ŝ^2Cr)	0.064	0.061	0.056	0.047
cov(Ŝ^2Cr,Ŝ^2Cr)	0.017	0.014	0.008	−0.001

---

Considering first the local particle number values, we find that the average number of electrons is quite consistent across the different spin states, staying just barely above 7 electrons (which corresponds to a Cr(III) oxidation state, with 2 doubly occupied ligand orbitals). However, the fluctuations in the oxidation state noticeably depend on the spin state, increasing as the global spin is decreased. This is easily understood as a consequence of the super-exchange mechanism, whereby coupling to electron-transfer states stabilizes the low-spin states relative to the high-spin states. Again, this is consistent with what is seen in ref. 90, and the more general treatment of formal magnetic interactions from Malrieu and coworkers.⁴³ For the global singlet state, the statistical correlation between the oxidation-state fluctuations on the two Cr centers is only 21.1%, indicating that the majority of the oxidation-state fluctuations are due to electron exchanges with the bridging ligands.

Because both Cr centers are re-coupled into global eigenvectors of Ŝ², the local Ŝ^z_Cr is no longer a good quantum number, and becomes maximally uncertain. In fact, we can further test how closely the system follows Heisenberg–Dirac–von Vleck physics by comparing the analytic values of the local S^z_Cr variance using the Clebsch–Gordon coefficients for a product of two spins. For the singlet, triplet, quintet, and heptet states, the analytic var(Ŝ^z_Cr) values are −1.25, −2.05, −1.25, and −0.45, respectively. Our computed correlations are only slightly different from these analytical values: −1.245, − 2.039, − 1.242, and − 0.449, further indicating good consistency with the Heisenberg model.

Inspecting the local Ŝ²_Cr values, we see a complementary picture to that provided by the particle-number fluctuations. As the global spin is decreased, the local S² values also tend to decrease, while the variance increases. This is consistent with the superexchange mechanism stabilizing global low-spin states by coupling to charge-transfer states. When an electron transfers from one Cr to the another, the number of unpaired electrons decreases. As a result, the local S²_Cr values decrease, and develop a positive covariance between the centers.

3.3 Conjugation in 2D

In the earlier sections, the tetracene tetramer (Sec. 3 3.1) served as an example of a completely non-bonded system, which is clearly quite easily clusterable. In Sec. 3 3.2, we demonstrated that the concepts of oxidation state and local spin allowed us to treat the Cr₂(III) complex in a clustered representation. In contrast, conjugated π-systems are characterized primarily by the highly delocalized nature of their electronic structure. In this section, we investigate the ability to compute and analyze the full π active space for a large delocalized π system.

Hexabenzocoronene (C₄₂H₁₈) is a polycyclic aromatic hydrocarbon (shown in Fig. 7(a)) where an additional benzene ring is fused to the outside of a central coronene ring. The π electrons are delocalized across the entire molecule, which contributes to its unique electronic properties. As we have seen previously,²³ the delocalized hexabenzocoronene system serves as a nice edge case for evaluating the ability of TPSCI to provide both accurate and insightful results for systems that are not obviously clusterable.

Fig. 7 Hexabenzocoronene. (a) Molecular structure and cluster indices. Active space (42e, 42o) includes all π orbitals. (b) Extrapolation of the energy of the singlet ground state. Units of milli-Hartree. (c) Charge covariance matrix, cov(N̂_I,N̂_J).

3.3.1 Active-space selection and clustering. For these results, we have considered the full π-system orbital active space (42 orbitals consisting of the 2p_z orbitals on each carbon) using the cc-pVDZ basis set.⁹¹ Viewing hexabenzocoronene as a collection of seven benzene rings, we partition the 42 orbitals into 7 clusters of 6 orbitals. A depiction of this clustering is shown in Fig. 7(a). The geometry is optimized at the B3LYP/cc-pVDZ level of theory. The cMF optimization (including orbital rotations between clusters) is performed using our open-source Julia package ClusterMeanField.jl.⁹² For each cluster, a local many-body basis was constructed using the embedded Schmidt truncation (EST) approach, where we define the cluster basis as the singular vectors of the FCI ground state on the cluster plus an orbital bath. We discarded Schmidt vectors with singular values smaller than a threshold value of 1 × 10⁻⁴. Detailed analysis of the EST approach has been carried out in our recent TPSCI paper.²³

3.3.2 Convergence of the TPSCI ground state. In Fig. 7(b), we plot the extrapolation of the TPSCI variational energy as a function of the PT2 correction. We use ε_FOIS = 1 × 10⁻⁶ (threshold on the external space couplings), and tightest ε_CIPSI = 1.5 × 10⁻⁴ through the bootstrapping HOSVD approach for this calculation.²⁴ Here we see that our ground-state TPSCI + PT2 energy is only about 5 mH away from the extrapolated result. Further, this was with a variational dimension of only around 114k. The extrapolated total energy of the molecule is −1601.6971 ± 0.0001 au. It's important to emphasize that the uncertainty here is due to the linear fit and does not imply a variational guarantee. If we knew that the relationship was indeed linear over the full range, then our energy would be correct to 0.1 mH.

We note that the performance of TPSCI on such delocalized π systems depends significantly on the topology of the molecule. For the complex considered here, there is a well defined Clar's structure that suggests a unique clustering. We expect this to be key to achieving accurate solutions. In contrast, our recent work²³ revealed that π systems without a well-defined clustering into a Clar's structure (such as coronene) are significantly slower to converge. We plan to explore this topic in the future for a more extensive set of polyaromatic hydrocarbons.

3.3.3 Correlations in between clusters. In Fig. 7(c), the charge covariances between the clusters are depicted as a heatmap, with each row/column corresponding to a given cluster labeled by the number on the diagonal. Looking first at the diagonal of the matrix (the charge variances), we see that the outside clusters (1–6) all have equivalent charge fluctuations, while the central benzene unit has significantly larger fluctuations in the ground state. Because the variance quantifies the uncertainty in the number of electrons in a given cluster, each of the outer and inner clusters has a number of electrons of 6.0 ± 0.5 and 6.0 ± 0.7, respectively (using a 3σ confidence interval). Considering next the off-diagonal matrix elements, we see that all nearest-neighbor cluster couplings are approximately the same. Assuming that two-body correlations dominate, this means that the central carbon should have a variance that is about twice that of the outer clusters, just based on the fact that it has twice as many nearest neighbors, which is consistent with the observed results. Because small differences are difficult to see in the heatmap, we have listed the unique covariance quantities in Table 5. Very similar results exist for the Ŝ^z_I correlations, as can be seen in Fig. S2 in the ESI.†

Table 5 Unique cluster correlations in hexabenzocoronene. N̂_I is the number operator for cluster I. Ŝ^z_I is the spin magnetization operator for cluster I. Q̂_I is the projector onto the orthogonal complement of the cMF ground state for cluster I. Cluster pair indices correspond to the labeling in Fig. 7(a) with the description of the interaction type in parenthesises

Cluster pair	cov(N̂i,N̂j)	cov(Ŝ^zi,Ŝ^zj)	cov(Q̂i,Q̂j)
1,1 (outer)	0.02551	0.00925	0.03307
7,7 (inner)	0.05729	0.02103	0.06867
1,2 (nearest neighbor)	−0.00759	−0.00276	0.01079
1,3 (meta)	−0.00018	−0.00005	−0.00010
1,4 (para)	−0.00040	−0.000134	−0.00006
1,7 (outer-inner)	−0.00956	−0.00350	0.01342

---

By considering the wavefunction directly, we notice that about 89.6% of the wavefunction is attributable to TPSs that have 6α and 6β in each cluster, whereas 9.1% is due to charge-transfer configurations, and 1.3% is due to local spin-flip configurations.

We also note that in Table 5, we see stronger particle number and spin correlations between benzenes connected in the para position than between those with meta connections, despite being further in distance, indicating a slight directing effect of the central benzene. However, the opposite is seen with the cluster excitation, Q̂_I, correlations.

Although the particle number covariance between each pair of neighboring clusters is negative (indicating charge transfer), the Ŝ²_I correlations are positive between neighboring clusters. This is consistent with the description of charge correlation. When an electron from a cluster hops into another cluster, then a doublet state will be formed in both of the clusters. One cluster will be one electron deficient, giving rise to a cationic doublet state, and the extra electron forms an anionic doublet state in the neighboring cluster. Consequently, when one cluster is in a doublet state, its neighbors have a higher probability of also being found in a doublet state, making the correlation positive. The entanglement between clusters leads each to acquire a non-zero average Ŝ²_I value. The outer clusters have total spin of 0.034 ± 0.55, whereas the central cluster as a local S²_I value of 0.078 ± 0.84, where uncertainty is given as 3σ.

4. Conclusions

In this paper, we have explored the ability of TPSCI to provide accurate yet interpretable approximations to FCI on relatively large orbital active spaces. Because TPSCI works in a basis consisting of products of local FCI states, the more separable a system is, the easier it should be to simulate. As such, in this paper, we consider three example systems, which range in the degree of separability: (i) a completely non-bonded tetramer of tetracene molecules, (ii) a more strongly interacting dichromium organometallic complex that, while bonded, is still characterized with local quantities like oxidation state, and (iii) a completely delocalized π system of hexabenzocoronene.

For the dichromium example, this was the first TPSCI calculation applied to open-shell biradical systems. We found that TPSCI was able to compute exchange coupling constants that are larger in magnitude (presumably more accurate) than recent DMRG calculations.

For each of these systems, we characterized the resulting wavefunctions using quantities that are easily accessible from the TPS basis. By leveraging the natural diabatic character of the TPS basis, we are able to easily construct Bloch effective Hamiltonians, which provide quantitative relationships between physically relevant degrees of freedom. This provided access to quantities such as the effective coupling between the bright states and the multiexcitonic states, 〈S₁|Ĥ^eff|¹TT〉, which implicitly includes the downfolded effects from charge-transfer couplings, which are substantially enhanced compared to the direct coupling.

We additionally used correlation functions of quantities like particle number, spin, and cluster excitation to provide a more detailed analysis of the various variational TPSCI eigenstates, and ultimately compare the results to inspection of the wavefunction itself, which is particularly interpretable due to the diabatic nature of the basis. This work helps lay out approaches for extracting more insight from TPSCI wavefunctions (and other TPS-based wavefunctions) in the future.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Science Foundation (Award No. 1752612). NMB also acknowledges the support provided by the Graduate School Doctoral Assistantship from the Department of Chemistry at Virginia Tech.

References

1. K. Ruedenberg, M. W. Schmidt, M. M. Gilbert and S. T. Elbert, Are atoms intrinsic to molecular electronic wavefunctions? I. The FORS model, Chem. Phys., 1982, 71, 41–49.
2. K. Ruedenberg, M. W. Schmidt, M. M. Gilbert and S. T. Elbert, Are atoms intrinsic to molecular electronic wavefunctions? III. Analysis of FORS configurations, Chem. Phys., 1982, 71, 65–78.
3. K. Ruedenberg, M. W. Schmidt and M. M. Gilbert, Are atoms sic to molecular electronic wavefunctions? II. Analysis of fors orbitals, Chem. Phys., 1982, 71, 51–64.
4. J. P. Foster and F. Weinhold, Natural hybrid orbitals, J. Am. Chem. Soc., 1980, 102, 7211–7218.
5. E. D. Glendening, C. R. Landis and F. Weinhold, Natural bond orbital methods, Wiley Interdiscip. Rev. Comput. Mol. Sci., 2012, 2, 1–42.
6. E. D. Glendening and A. Streitwieser, Natural energy decomposition analysis: An energy partitioning procedure for molecular interactions with application to weak hydrogen bonding, strong ionic, and moderate donor–acceptor interactions, J. Chem. Phys., 1994, 100, 2900–2909.
7. A. E. Reed, R. B. Weinstock and F. Weinhold, Natural population analysis, J. Chem. Phys., 1985, 83, 735–746.
8. R. Z. Khaliullin, M. Head-Gordon and A. T. Bell, An efficient self-consistent field method for large systems of weakly interacting components, J. Chem. Phys., 2006, 124, 204105.
9. R. Z. Khaliullin, E. A. Cobar, R. C. Lochan, A. T. Bell and M. Head-Gordon, Unravelling the Origin of Intermolecular Interactions Using Absolutely Localized Molecular Orbitals, J. Phys. Chem. A, 2007, 111, 8753–8765.
10. Y. Mao, M. Loipersberger, P. R. Horn, A. Das, O. Demerdash, D. S. Levine, S. Prasad Veccham, T. Head-Gordon and M. Head-Gordon, From Intermolecular Interaction Energies and Observable Shifts to Component Contributions and Back Again: A Tale of Variational Energy Decomposition Analysis, Annu. Rev. Phys. Chem., 2021, 72, 641–666,  DOI:10.1146/annurev-physchem-090419-115149.
11. R. F. W. Bader, A quantum theory of molecular structure and its applications, Chem. Rev., 1991, 91, 893–928.
12. R. J. Cave and M. D. Newton, Calculation of electronic coupling matrix elements for ground and excited state electron transfer reactions: Comparison of the generalized Mulliken–Hush and block diagonalization methods, J. Chem. Phys., 1997, 106, 9213–9226.
13. R. J. Cave and M. D. Newton, Generalization of the Mulliken-Hush treatment for the calculation of electron transfer matrix elements, Chem. Phys. Lett., 1996, 249, 15–19.
14. H.-C. Chen, Z.-Q. You and C.-P. Hsu, The mediated excitation energy transfer: Effects of bridge polarizability, J. Chem. Phys., 2008, 129, 084708.
15. S. Fatehi, E. Alguire and J. E. Subotnik, Derivative couplings and analytic gradients for diabatic states, with an implementation for Boys-localized configuration-interaction singles, J. Chem. Phys., 2013, 139, 124112.
16. C.-P. Hsu, Z.-Q. You and H.-C. Chen, Characterization of the Short-Range Couplings in Excitation Energy Transfer, J. Phys. Chem. C, 2008, 112, 1204–1212.
17. K. Ruedenberg and G. J. Atchity, A quantum chemical determination of diabatic states, J. Chem. Phys., 1993, 99, 3799–3803.
18. J. E. Subotnik, S. Yeganeh, R. J. Cave and M. A. Ratner, Constructing diabatic states from adiabatic states: Extending generalized Mulliken–Hush to multiple charge centers with Boys localization, J. Chem. Phys., 2008, 129, 244101.
19. J. Zheng, Y. K. Kang, M. J. Therien and D. N. Beratan, Generalized Mulliken-Hush analysis of electronic coupling interactions in compressed pi-stacked porphyrin-bridge-quinone systems, J. Am. Chem. Soc., 2005, 127, 11303–11310.
20. R. Baer and M. Head-Gordon, Sparsity of the Density Matrix in Kohn-Sham Density Functional Theory and an Assessment of Linear System-Size Scaling Methods, Phys. Rev. Lett., 1997, 79, 3962–3965.
21. V. Abraham and N. J. Mayhall, Cluster many-body expansion: A many-body expansion of the electron correlation energy about a cluster mean field reference, J. Chem. Phys., 2021, 155, 054101.
22. V. Abraham and N. J. Mayhall, Coupled Electron Pair-Type Approximations for Tensor Product State Wave Functions, J. Chem. Theory Comput., 2022, 18, 4856–4864.
23. V. Abraham and N. J. Mayhall, Selected Configuration Interaction in a Basis of Cluster State Tensor Products, J. Chem. Theory Comput., 2020, 16, 6098–6113.
24. N. M. Braunscheidel, V. Abraham and N. J. Mayhall, Generalization of the Tensor Product Selected CI Method for Molecular Excited States, J. Phys. Chem. A, 2023, 127, 8179–8193.
25. N. J. Mayhall, Using Higher-Order Singular Value Decomposition To Define Weakly Coupled and Strongly Correlated Clusters: The n -Body Tucker Approximation, J. Chem. Theory Comput., 2017, 13, 4818–4828.
26. N. P. Bauman, E. J. Bylaska, S. Krishnamoorthy, G. H. Low, N. Wiebe, C. E. Granade, M. Roetteler, M. Troyer and K. Kowalski, Downfolding of many-body Hamiltonians using active-space models: Extension of the sub-system embedding sub-algebras approach to unitary coupled cluster formalisms, J. Chem. Phys., 2019, 151, 014107.
27. S. M. Parker, T. Seideman, M. A. Ratner and T. Shiozaki, Communication: Active-space decomposition for molecular dimers, J. Chem. Phys., 2013, 139, 021108.
28. S. M. Parker and T. Shiozaki, Quasi-diabatic States from Active Space Decomposition, J. Chem. Theory Comput., 2014, 10, 3738–3744.
29. S. M. Parker, T. Seideman, M. A. Ratner and T. Shiozaki, Model Hamiltonian Analysis of Singlet Fission from First Principles, J. Phys. Chem. C, 2014, 118, 12700–12705.
30. C. A. Jiménez-Hoyos and G. E. Scuseria, Cluster-based mean-field and perturbative description of strongly correlated fermion systems: Application to the one- and two-dimensional Hubbard model, Phys. Rev. B, 2015, 92, 085101.
31. A. Papastathopoulos-Katsaros, C. A. Jiménez-Hoyos, T. M. Henderson and G. E. Scuseria, Coupled Cluster and Perturbation Theories Based on a Cluster Mean-Field Reference Applied to Strongly Correlated Spin Systems, J. Chem. Theory Comput., 2022, 18, 4293–4303.
32. A. Papastathopoulos-Katsaros, T. M. Henderson and G. E. Scuseria, Linear combinations of cluster mean-field states applied to spin systems, arXiv, 2024, arXiv:2402.06760 [cond-mat, physics:physics] version: 1 DOI:10.48550/arXiv.2402.06760.
33. A. Papastathopoulos-Katsaros, T. M. Henderson and G. E. Scuseria, Symmetry-projected cluster mean-field theory applied to spin systems, J. Chem. Phys., 2023, 159, 084107.
34. M. R. Hermes and L. Gagliardi, The Variational Localized Active Space Self-Consistent Field Method, arXiv, 2020, arXiv:2003.02995 [physics] DOI:10.48550/arXiv.2003.02995.
35. S. Li, Block-correlated coupled cluster theory: The general formulation and its application to the antiferromagnetic Heisenberg model, J. Chem. Phys., 2004, 120, 5017.
36. B. Huron, J. P. Malrieu and P. Rancurel, Iterative perturbation calculations of ground and excited state energies from multiconfigurational zeroth-order wavefunctions, J. Chem. Phys., 1973, 58, 5745.
37. A. A. Holmes, N. M. Tubman and C. J. Umrigar, Heat-Bath Configuration Interaction: An Efficient Selected Configuration Interaction Algorithm Inspired by Heat-Bath Sampling, J. Chem. Theory Comput., 2016, 12, 3674–3680.
38. C. Bloch, Sur la théorie des perturbations des états liés, Nucl. Phys., 1958, 6, 329–347.
39. J. des Cloizeaux and J. J. Pearson, Spin-Wave Spectrum of the Antiferromagnetic Linear Chain, Phys. Rev., 1962, 128, 2131–2135.
40. C. J. Calzado, J. Cabrero, J. P. Malrieu and R. Caballol, Analysis of the magnetic coupling in binuclear complexes. II. Derivation of valence effective Hamiltonians from ab initio CI and DFT calculations, J. Chem. Phys., 2002, 116, 3985–4000.
41. R. Maurice, R. Bastardis, C. d. Graaf, N. Suaud, T. Mallah and N. Guihéry, Universal Theoretical Approach to Extract Anisotropic Spin Hamiltonians, J. Chem. Theory Comput., 2009, 5, 2977–2984.
42. A. Monari, D. Maynau and J.-P. Malrieu, Determination of spin Hamiltonians from projected single reference configuration interaction calculations. I. Spin 1/2 systems, J. Chem. Phys., 2010, 133, 044106.
43. J. P. Malrieu, R. Caballol, C. J. Calzado, C. de Graaf and N. Guihéry, Magnetic Interactions in Molecules and Highly Correlated Materials: Physical Content, Analytical Derivation, and Rigorous Extraction of Magnetic Hamiltonians, Chem. Rev., 2014, 114, 429–492.
44. N. J. Mayhall, P. R. Horn, E. J. Sundstrom and M. Head-Gordon, Spin-flip non-orthogonal configuration interaction: a variational and almost black-box method for describing strongly correlated molecules, Phys. Chem. Chem. Phys., 2014, 16, 22694–22705.
45. N. J. Mayhall and M. Head-Gordon, Computational Quantum Chemistry for Multiple-Site Heisenberg Spin Couplings Made Simple: Still Only One Spin–Flip Required, J. Phys. Chem. Lett., 2015, 6, 1982–1988.
46. N. J. Mayhall, From Model Hamiltonians to ab Initio Hamiltonians and Back Again: Using Single Excitation Quantum Chemistry Methods To Find Multiexciton States in Singlet Fission Materials, J. Chem. Theory Comput., 2016, 12, 4263–4273.
47. P. Pokhilko and A. I. Krylov, Effective Hamiltonians derived from equation-of-motion coupled-cluster wave functions: Theory and application to the Hubbard and Heisenberg Hamiltonians, J. Chem. Phys., 2020, 152, 094108.
48. D. Casanova and A. I. Krylov, Quantifying local exciton, charge resonance, and multiexciton character in correlated wave functions of multichromophoric systems, J. Chem. Phys., 2016, 144, 014102.
49. A. V. Luzanov, D. Casanova, X. Feng and A. I. Krylov, Quantifying charge resonance and multiexciton character in coupled chromophores by charge and spin cumulant analysis, J. Chem. Phys., 2015, 142, 224104.
50. Q. Sun, et al., Recent developments in the PySCF program package, J. Chem. Phys., 2020, 153, 024109.
51. N. J. Mayhall, V. Abraham, N. M. Braunscheidel and A. Bachhar, FermiCG, 2024, https://github.com/mayhallgroup/FermiCG, accessed 03-03-2024.
52. M. B. Smith and J. Michl, Recent Advances in Singlet Fission, Annu. Rev. Phys. Chem., 2013, 64, 361–386.
53. M. C. Hanna and A. J. Nozik, Solar conversion efficiency of photovoltaic and photoelectrolysis cells with carrier multiplication absorbers, J. Appl. Phys., 2006, 100, 074510.
54. J. C. Johnson, Open questions on the photophysics of ultrafast singlet fission, Commun. Chem., 2021, 4, 1–3.
55. R. D. Pensack, A. J. Tilley, S. R. Parkin, T. S. Lee, M. M. Payne, D. Gao, A. A. Jahnke, D. G. Oblinsky, P.-F. Li, J. E. Anthony, D. S. Seferos and G. D. Scholes, Exciton delocalization drives rapid singlet fission in nanoparticles of acene derivatives, J. Am. Chem. Soc., 2015, 137, 6790–6803.
56. G. D. Scholes, Correlated Pair States Formed by Singlet Fission and Exciton-Exciton Annihilation, J. Phys. Chem. A, 2015, 119, 12699–12705.
57. T. C. Berkelbach, M. S. Hybertsen and D. R. Reichman, Microscopic theory of singlet exciton fission. III. Crystalline pentacene, J. Chem. Phys., 2014, 141, 074705.
58. T. C. Berkelbach, in Advances in Chemical Physics, Vol 162, ed. S. A., Rice, A. R., Dinner, 2017, vol. 162, pp. 1–38.
59. Y. Shu, E. G. Hohenstein and B. G. Levine, Configuration interaction singles natural orbitals: an orbital basis for an efficient and size intensive multireference description of electronic excited states, J. Chem. Phys., 2015, 142, 024102.
60. R. Ditchfield, W. J. Hehre and J. A. Pople, Self-Consistent Molecular-Orbital Methods. IX. An Extended Gaussian-Type Basis for Molecular-Orbital Studies of Organic Molecules, J. Chem. Phys., 1971, 54, 724–728.
61. S. Lehtola and H. Jónsson, Pipek–Mezey Orbital Localization Using Various Partial Charge Estimates, J. Chem. Theory Comput., 2014, 10, 642–649.
62. T. C. Berkelbach, M. S. Hybertsen and D. R. Reichman, Microscopic theory of singlet exciton fission. II. Application to pentacene dimers and the role of superexchange, J. Chem. Phys., 2013, 138, 114103.
63. V. Abraham and N. J. Mayhall, Simple Rule To Predict Boundedness of Multiexciton States in Covalently Linked Singlet-Fission Dimers, J. Phys. Chem. Lett., 2017, 8, 5472–5478.
64. A. Caneschi, D. Gatteschi, R. Sessoli, A. L. Barra, L. C. Brunel and M. Guillot, Alternating current susceptibility, high field magnetization, and millimeter band EPR evidence for a ground S = 10 state in [Mn12O12(Ch3COO)16(H2O)4].2CH3COOH.4H2O, J. Am. Chem. Soc., 1991, 113, 5873–5874.
65. S. Demir, M. I. Gonzalez, L. E. Darago, W. J. Evans and J. R. Long, Giant coercivity and high magnetic blocking temperatures for N2 3- radical-bridged dilanthanide complexes upon ligand dissociation, Nat. Commun., 2017, 8, 2144.
66. M. Gil-Sepulcre and A. Llobet, Molecular water oxidation catalysts based on first-row transition metal complexes, Nat. Catal., 2022, 5, 79–82.
67. Y. Tanabe and Y. Nishibayashi, Recent advances in catalytic nitrogen fixation using transition metal–dinitrogen complexes under mild reaction conditions, Coord. Chem. Rev., 2022, 472, 214783.
68. This is true for the Hubbard model at half-filling. Situations with away-from half-filling can be mapped onto different Hamiltonians, such as the tJ model.
69. L. Noodleman, Valence bond description of antiferromagnetic coupling in transition metal dimers, J. Chem. Phys., 1981, 74, 5737–5743.
70. K. Yamaguchi, F. Jensen, A. Dorigo and K. N. Houk, A spin correction procedure for unrestricted Hartree–Fock and Møller-Plesset wavefunctions for singlet diradicals and polyradicals, Chem. Phys. Lett., 1988, 149, 537–542.
71. S. Yamanaka, T. Kawakami, H. Nagao and K. Yamaguchi, Effective exchange integrals for open-shell species by density functional methods, Chem. Phys. Lett., 1994, 231, 25–33.
72. M. Reiher, O. Salomon and B. Artur Hess, Reparameterization of hybrid functionals based on energy differences of states of different multiplicity, Theor. Chem. Acc., 2001, 107, 48–55.
73. M. Reiher, Theoretical Study of the Fe(phen)2(NCS)2 Spin-Crossover Complex with Reparametrized Density Functionals, Inorg. Chem., 2002, 41, 6928–6935.
74. D. A. Pantazis, Meeting the Challenge of Magnetic Coupling in a Triply-Bridged Chromium Dimer: Complementary Broken-Symmetry Density Functional Theory and Multireference Density Matrix Renormalization Group Perspectives, J. Chem. Theory Comput., 2019, 15, 938–948.
75. S. Szalay, M. Pfeffer, V. Murg, G. Barcza, F. Verstraete, R. Schneider and r. Legeza, Tensor product methods and entanglement optimization for ab initio quantum chemistry, Int. J. Quantum Chem., 2015, 115, 1342–1391.
76. A. Baiardi and M. Reiher, The density matrix renormalization group in chemistry and molecular physics: Recent developments and new challenges, J. Chem. Phys., 2020, 152, 040903.
77. U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Ann. Phys., 2011, 326, 96–192.
78. S. Wouters and D. Van Neck, The density matrix renormalization group for ab initio quantum chemistry, Eur. Phys. J. D, 2014, 68, 272.
79. J. J. Eriksen, et al., The Ground State Electronic Energy of Benzene, J. Phys. Chem. Lett., 2020, 11, 8922–8929.
80. R. Olivares-Amaya, W. Hu, N. Nakatani, S. Sharma, J. Yang and G. K.-L. Chan, The ab-initio density matrix renormalization group in practice, J. Chem. Phys., 2015, 142, 034102.
81. S. Sharma, K. Sivalingam, F. Neese and G. K.-L. Chan, Low-energy spectrum of iron-sulfur clusters directly from many-particle quantum mechanics, Nat. Chem., 2014, 6, 927–933.
82. S. Sharma, G. Jeanmairet and A. Alavi, Quasi-degenerate perturbation theory using matrix product states, J. Chem. Phys., 2016, 144, 034103.
83. T. V. Harris, Y. Kurashige, T. Yanai and K. Morokuma, Ab initio density matrix renormalization group study of magnetic coupling in dinuclear iron and chromium complexes, J. Chem. Phys., 2014, 140, 054303.
84. N. J. Mayhall and M. Head-Gordon, Computational quantum chemistry for single Heisenberg spin couplings made simple: Just one spin flip required, J. Chem. Phys., 2014, 141, 134111.
85. A. Niemann, U. Bossek, K. Wieghardt, C. Butzlaff, A. X. Trautwein and B. Nuber, A New Structure–Magnetism Relationship for Face-Sharing Transition-Metal Complexes with d3–d3 Electronic Configuration, Angew. Chem., Int. Ed. Engl., 1992, 31, 311–313.
86. For each active space, we solve for the lowest M states in the lowest M_s-subspace, and then generate the higher M_s vectors by directly applying Ŝ±. This is done to, not only reduce the cost, but also to ensure that every TPS in our basis can be properly spin-adapted to form proper eigenstates of Ŝ², as first noted in ref. 24.
87. H. G. A. Burton and P.-F. Loos, Rationale for the Extrapolation Procedure in Selected Configuration Interaction, arXiv, 2024, arXiv:2312.12530 [cond-mat, physics:nucl-th, physics:physics] DOI:10.48550/arXiv.2312.12530.
88. A. A. Holmes, C. J. Umrigar and S. Sharma, Excited states using semistochastic heat-bath configuration interaction, J. Chem. Phys., 2017, 147, 164111.
89. We have also explored directly extrapolating the J values themselves, and found the results to be essentially the same.
90. R. Pandharkar, M. R. Hermes, C. J. Cramer and L. Gagliardi, Localized Active Space-State Interaction: a Multireference Method for Chemical Insight, J. Chem. Theory Comput., 2022, 18, 6557–6566.
91. T. H. Dunning, Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen, J. Chem. Phys., 1989, 90, 1007–1023.
92. N. J. Mayhall, V. Abraham, N. M. Braunscheidel and A. Bachhar, ClusterMeanField.jl, 2023, https://github.com/nmayhall-vt/ClusterMeanField.jl, accessed 11-24-2023.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4fd00049h

---