Restricted Open-shell cluster Mean-Field theory for Strongly Correlated Systems

In most cluster-based methods, an orbital active space is partitioned into non-overlapping groups, referred to as clusters (subsets of the total available single-fermion states) based on some desired property (e.g., locality, symmetry, etc.). We will index each of these disjoint orbital subsets, clusters, with a Roman index, $I$. Each cluster, $I$, supports a Fock space, for which we will define a “cluster basis”, indexed using Greek letters, $\ket{I_{\alpha}}$. Generally, each cluster state will be written as a linear combination of all possible Slater determinants constructed out of the cluster’s orbitals. However, this is not a strict requirement, and it is possible to use more sophisticated parameterizations of the local cluster states. Assuming each cluster’s Fock space is un-truncated, a basis for the full global Hilbert space can be constructed by forming all possible tensor products of local cluster states. As previously done [23], we further choose to define our cluster states to be eigenvectors of $\hat{N}$ and $\hat{S}_{z}$, which will ensure our quantum states have well-defined local quantum numbers, a feature that will simplify enforcing global symmetries. We will then further index each cluster state with the sector of Fock space it belongs to, $\ket{I_{\alpha}^{{\tt n}_{I}}}$, where $\alpha$ runs over all states in the local Fock sector, ${{\tt n}_{I}}$ of the $I^{th}$ cluster. Each global tensor product state (TPS) serving as a basis vector can be expressed using these local many-body cluster states:

where $\vec{n}=\left({\tt n}_{1},{\tt n}_{2},\dots,{\tt n}_{N}\right)$ is a vector index running over all possible combinations of local Fock sectors (i.e., distributions of electrons among the clusters) such that $\ket{I_{\alpha}^{{\tt n}_{I}}}$ spans the entire Fock space of cluster $I$. These quantum number strings represent lists of eigenvalues of the operators $\hat{N}$ and $\hat{S}_{z}$ for each cluster in the system.[21] The exact wavefunction can then be expressed as a linear combination of such TPS’s:[23]

where $C^{\vec{n}}_{\alpha,\beta,\dots,\gamma}$ is the coefficient tensor. While the formalism presented is general and formally exact, as discussed in the previous section, this representation is only expected to be compact when the interactions within a cluster are stronger than those between clusters, allowing the basis vectors to incorporate a relatively large amount of electron correlation embedded within the local many-body cluster states. As a result, the coefficient tensor, $C^{\vec{n}}_{\alpha,\beta,\dots,\gamma}$ only needs to describe inter-cluster correlation, with all intra-cluster correlation folded into the basis vectors. The resulting wavefunction written in terms of the TPS’s then requires fewer basis vectors than in the traditional Slater determinant basis. By restricting the sum over $\vec{n}$ to only those Fock sector configurations that have the correct total number of alpha and beta electrons, we naturally preserve total particle number, $N$, and spin projection, $S_{z}$, symmetries in our implementation.

Following the formalism defined in the ASD method [17], the standard electronic Hamiltonian in the second quantized form:

can be partitioned into contributions based on the number of distinct clusters involved: one-, two-, three-, and four-cluster terms. These contributions are defined as follows:

In Eq. (II.1), $\hat{H_{I}}$ includes terms where all creation and annihilation operators are within cluster $I$, $\hat{H}_{IJ}$ involves operators from both clusters $I$ and $J$, and so on. As the ab initio Hamiltonian consists solely of two-body interactions, the maximum number of clusters involved in the interactions is limited to four.

The above notation used for describing our tensor product space is general for any complete set of local cluster states. However, the compactness of global states represented in this basis ultimately is determined by the specific states chosen. In principle, we could choose our cluster states, $\ket{I_{\alpha}^{{\tt n}_{I}}}$, to be those which diagonalize the cluster’s local Hamiltonian, i.e., $\hat{H}_{I}\ket{I_{\alpha}^{{\tt n}_{I}}}=E_{\alpha}^{{\tt n}_{I}}\ket{I_{\alpha}^{{\tt n}_{I}}}$. This would then fold all local electron correlations into the cluster states. However, the interactions between clusters ultimately affect the correlation inside of a cluster, thus neglecting all inter-cluster interactions would not create the most efficient representation. Adopting the formalism proposed by Jiménez-Hoyos and Scuseria[1], we choose to define our cluster states by minimizing the energy of a single TPS wavefunction:

with respect to variations in both the cluster states and the orbitals defining the clusters themselves. Analogous to HF theory, this provides a mean-field treatment of all inter-cluster correlations, while treating all intra-cluster correlations explicitly.

In addition to providing the lowest energy reference TPS, this approach additionally offers a more reproducible method for defining orbital clusters rather than relying on arbitrary localization criteria as our approach is based on a variational principle.

The simplest typical model used to describe magnetic behavior is the phenomenological Heisenberg-Dirac-van Vleck (HDvV) Hamiltonian, which can be derived from the Hubbard model at half-filling using quasi-degenerate perturbation theory.¹¹1For cases away from half-filling, the system can be described by different Hamiltonians, such as the $tJ$ or double exchange models. For systems featuring two magnetic centers, $a$ and $b$, the $\hat{H}^{HDvV}$ Hamiltonian is written simply as the product of the spin operators on two centers:

where $\hat{\vec{S}}_{a}=\hat{S}_{a}^{x}\vec{x}+\hat{S}_{a}^{y}\vec{y}+\hat{S}_{a}^{z}\vec{z}$, are the spin operators associated with magnetic center $a$, and $J$ is the strength of coupling between two spin centers. The exchange coupling constant, $J$, dictates the nature of spin alignment, being positive for ferromagnetic (F) and negative for antiferromagnetic (AF) interactions, with its magnitude indicative of interaction strength. These $J$ values that parameterize the HDvV Hamiltonian, completely determine the resulting energy spectrum of low-energy spin states, which for a two-center system is given by the Landé interval rule:

When derived from the Hubbard model with hopping strength $t$ and on-site coulomb repulsion $U$, the second order contribution to the exchange coupling constant is $J=\frac{-4t^{2}}{U}$. This so-called “kinetic exchange” highlights the role that electron delocalization or charge resonance (quantified by the hopping integral, $t$) has toward increasing antiferromagnetic coupling strengths. However, as discussed in Ref. 37, the zeroth order ab-initio Hamiltonian also contains contributions from non-local direct exchange, $K$, in addition to the second order kinetic exchange coming from the Hubbard model. As bare direct exchange favors high-spin states and the second-order kinetic exchange term favors low-spin states, even determining the correct sign for the exchange coupling constant $J$ can be challenging also.

While $J$ is a useful quantity for rationalizing multicenter complexes, ultimately the HDvV Hamiltonian is an approximate, simplified description of the true electronic structure. As such, ab initio calculations are invaluable for not only computing values of $J$ from energy differences between computed states but also for determining the suitability of a phenomenological Hamiltonian for a given complex. However, computing these low-energy states accurately is a well documented challenge for computational chemistry. Many different factors contribute to the spin states energies [38] of multicenter organometallic complexes with open-shell transition metal centers. An important factor is the modulation of magnetic interactions by the ligands. In instances where metal ions are spatially separated by a linear bridging ligand, a “direct” metal–metal interaction is absent. However, in cases of multiply bridged dimers, the interaction strengths can change due to the variable metal-metal distances permitted by the bridging topology. This proximity may result in antiferromagnetic coupling strength deviating from expectations based solely on the chemical nature of the bridges, leading to intricacies in the description of the magnetic interaction. Firstly, orbital mixing between metal and bridging ligand facilitates the delocalization of unpaired electrons and hence thereby increases the magnetic coupling. Dynamical spin and charge polarization can further affect the coupling strengths. Conjugated long bridging units have low-lying $\pi$-$\pi^{*}$ valence excited states, which can lead to a high degree of spin and charge polarizations.

Because of all the myriad of contributions affecting these low-energy spin states, weakly interacting metal centers with unpaired electrons exhibit high degrees of multiconfigurational complexity. As mentioned in the introductory section, conventional approaches such as perturbation theory or coupled cluster theory are unsuitable because they rely on a qualitatively accurate single Slater determinant wavefunction as a reference. Given that the single-determinant representation only provides direct access to the high-spin (HS) state, $J$ values in DFT are typically derived using broken-symmetry (BS) and high-spin (HS) states [39, 40]. Although this approach is routinely capable of providing qualitative accuracy, the dependency of results on the chosen DFT functionals needs careful consideration, with quantitative accuracy often depending on the specific functional employed, making it impossible to systematically improve the results.

In this article, we apply RO-cMF theory to obtain clustered representations for multi-radical systems such as these kinds of transition metal complexes, with an eye toward providing a compact reference state for post-cMF methods, such as TPSCI [23], TPS-CEPA [22], etc. We have computed the exchange coupling constants for $[$Fe₂OCl${}_{6}]^{2-}$, $[$L₂Cr(III)${}_{2}(\mu-$OH$)_{3}]^{3+}$, L = N,N^′,N^′′- trimethyl-1,4,7-triazacyclononane, and phenalenyl-dimers using both RO-cMF with State Mixing, as well as with a PT2 corrected RO-cMF-PT2.

The exchange coupling within $[$Fe₂OCl${}_{6}]^{2-}$ (see Fig. 4(a)) has been investigated with various theoretical methods, including unrestricted HF (UHF),[41] density functional theory (DFT),[42, 43] internally contracted MRCI (IC-MRCI) [44], and DMRG [45, 46].

In a previous study, Morokuma et al.[45] investigated the impact of both basis sets and active-space selection on the calculated $J$ values for dinuclear complexes using density matrix renormalization group algorithm (DMRG), which has become a standard benchmark method for computing exchange coupling constants in transition metal complexes [47, 48, 49, 50, 51, 52, 53, 54, 55, 45]. However, if one needs only the value of $J$ or just the highest couple spin states, spin-flip methods can be highly effective as well [56, 57, 58, 59]. More recently [37], we have illustrated that TPSCI can also be used to compute low-energy states of organometallic compounds.[37]

One of the challenges in using active space methods is the need for the user to choose which orbitals to include. We have attempted to partially automate this process for this work. We start by optimizing the ROHF wavefunction for the high-spin ($S=5$) in 6-31G* basis. We have constructed a 52 orbital active space that consists primarily of the $3d$ orbitals of each Fe center, the $2p$ and $3p$ bridging oxygen orbitals, and the $3p$ and $4p$ orbitals of the six chlorine atoms. This is done by separately projecting the occupied, open-shell, and virtual ROHF orbitals onto the associated orthogonalized atomic orbital functions, choosing the largest overlapping singular vector from each of the ROHF subspaces. The orbitals associated with each atomic center are grouped into separate clusters (nine in total).

Using the described clustering, we get two sextet clusters and seven singlet clusters. As described above, to obtain the density matrices used in RO-cMF, we simply average over the 1RDMs for all degenerate $M_{s}$ values of each cluster’s ground state. For this system, this corresponds to six $M_{s}$ values for the two Fe clusters ($\frac{5}{2}$, $\frac{3}{2}$, $\frac{1}{2}$, $-\frac{1}{2}$, $-\frac{3}{2}$, and $-\frac{5}{2}$) and only a single $M_{s}$ value for the singlet ligand clusters. After optimizing the RO-cMF wavefunction, the resulting orbitals are shown in Fig. (3).

By diagonalizing the Hamiltonian in the relevant ($M_{s}=0$) sub-block (see the state mixing approach described in Fig. 2(b)), the spin-state energies are computed and shown in Fig. 4(a). By neglecting all ionic terms, this zeroth-order approximation incorrectly predicts a ferromagnetic alignment, with a positive $J$ value ($J=26.5$ cm^-1) derived using the Landé interval rule. When we apply perturbation theory to the RO-cMF states (Fig. 4(b), which naturally includes the inter-cluster hopping terms, the spin state ordering is reversed to the correct antiferromagnetic alignment, yielding a $J$ value of $-11.38$ cm^-1.²²2We note that both of these numbers correspond to the $E(S=0)-E(S=1)$ energy gap. However, we found a negligible deviation from the Landé intervals, and so all gaps would give nearly identical exchange coupling constants, with the largest differences on the order of $0.1$ cm^-1. While this is still far from the experimentally derived value of -117 cm^-1 [61], even the approximate inclusion of the ionic terms is able to provide qualitatively correct $J$ values. More accurate results can be obtained by treating these ionic terms non-perturbatively, as done in TPSCI [37].

Examination of chromium(III) dimers, as a paradigmatic class of antiferromagnetically coupled systems, has revealed ambiguities regarding the optimal description of their magnetic properties. Initially, ligand-mediated superexchange[62] was proposed as the predominant mechanism, but empirical observations revealed a correlation between the antiferromagnetic coupling strength and the metal-metal distance in octahedral Cr₂(III) complexes that share faces. This observation prompted the proposition of through-space interaction as the principal mechanism governing antiferromagnetic coupling in these complexes.[63, 47]

As we mentioned in the previous section, DFT has been instrumental in analyzing exchange-coupled systems. However, a noteworthy challenge emerges in the context of Cr₂(III) complexes, where DFT reportedly falls short in providing a qualitative description of the antiferromagnetic coupling. While BS-DFT demonstrated success in electronically similar high-valent manganese complexes, its reported failure for Cr₂(III) complexes highlights the challenges in distinguishing between different coupling mechanisms, such as direct exchange versus ligand-mediated superexchange [64, 65]. A recent study using BS-DFT, CASSCF, and DMRG by Pantazis et al.[47] has concluded that the dominance of direct through-space interaction is corroborated, yet the additional role of superexchange introduces an additional contribution to the overall magnetic behavior.

Tris-$\mu-$hydroxo Cr₂(III) (Fig. 5) is a face-sharing $d^{3}-d^{3}$ complex. For this complex, we decided to use five total clusters: one for each of the two Cr(III) centers, and one for each of the three bridging OH^-1 ligands, yielding two quartet clusters and three singlet clusters. Analogous to the Fe complex above, each Cr(III) cluster is a quartet, requiring an average of the 1RDMs over all the degenerate $M_{s}$ microstates {$\frac{3}{2}$, $\frac{1}{2}$, $-\frac{1}{2}$, $-\frac{3}{2}$}. Using the same projection based approach, used for the Fe(III) complex, we defined each of the Cr(III) clusters by projecting onto both $3d$ and $4d$ orbitals, and each bridging OH cluster by projecting onto $2p$ and $3p$ oxygen orbitals, resulting in a total active space of (32e, 38o).

We note that during the projection procedure, the Cr(III) clusters ended up adding 2 doubly occupied orbitals, so each Cr(III) cluster is defined as a local (7e, 10o) active space instead of the expected (3e, 10o). The RO-cMF orbitals are shown in Fig. (6). As shown in Fig. 5, the state mixing of RO-cMF states yields a ferromagnetic $J$ value of $14.7$ cm^-1, while RO-cMF-PT2 corrects this change this into an antiferromagnetic $J$ value of $-10.2$ cm^-1. This is nearly identical to the results obtained by Gagliardi and coworkers in Ref.  33, using a non-perturbative vLASSCF-SI approach, on a simplified model of this same complex. Again, while RO-cMF-PT2 is not sufficient for quantitative accuracy, the qualitative description is insightful, and can be made quantitatively accurate by more sophisticated post-cMF treatments like TPSCI, which in Ref. 37 found that TPSCI increased the magnitude of $J$ to $-31.3$ cm^-1, in very close agreement with CASSCF-NEVPT2 results of $-31.8$ cm^-1 [47].

Organic radicals are molecules or molecular fragments containing one or more unpaired electrons. These unpaired electrons make organic radicals highly reactive species, influencing their chemical behavior and potential applications in organic electronics, spintronics, and molecular magnetism. For example, they have been explored as components in organic light-emitting diodes (OLEDs),[66] organic field-effect transistors (OFETs),[67] and organic photovoltaics (OPVs).[68]

Phenalenyl is a polyaromatic hydrocarbon $\pi$-radical characterized by its odd-alternant structure and relative stability. Because the unpaired electron is able to delocalize throughout the conjugated $\pi$ system, the molecule exhibits a unique stabilization relative to other organic radicals. Despite this stability, phenalenyl radical and its substituted derivatives readily dimerize, resulting in a unique (2e, 12c) $\pi$-$\pi$ stacking bonding interaction of two radical units, termed “pancake bonding” [69]. This yields a singlet biradicaloid state where the two unpaired electrons from the multicenter SOMO orbitals of the radical monomers spin-recouple. This results in closer distances and stronger binding compared to standard dispersion bound complexes. The stabilization of the $\pi$ stacking configuration through pancake bonding renders it energetically competitive with $\sigma$-bonding. In $\pi$ stacking, the hexagonal arrangement of SOMO allows for both eclipse and staggered stacking. However, shorter $\pi$-$\pi$ bonding distance favors the staggered stacking more than eclipse which also gets destabilized due to smaller atom-atom repulsions.

In this section, we use RO-cMF for computing the effective exchange coupling constants for both phenyl 2,5,8-substituted and t-butyl (${}^{t}Bu$)-substituted phenalenyl radical dimers. In these two examples, we include all 26 $\pi$ electrons in the active space, which is then partitioned into two clusters, one for each phenalenyl radical. For the phenyl substituted dimer, we also include the $\pi$ electrons present in phenyl R-groups, as the radical can further delocalize to some extent onto the R-groups. Consequently, the active space for the phenyl substituted system is (62e, 62o), broken up into six (6e, 6o) clusters and two (13e, 13o) clusters. In contrast, because the unpaired electron is expected to be localized to the central phenalenyl units on the t-butyl substituted system, we have only considered monomer radical units as clusters that make (26e, 26o) active space, partitioned into two (13e, 13o) clusters.

Following the same procedure as with the organometallic complexes, we obtain a zeroth-order manifold of spin-states by diagonalizing the TPS’s with non-zero occupations in the RO-cMF density matrix (see Fig. 2(b)), providing spin-pure, yet physically deficient singlet and triplet states, from which an effective exchange coupling constant can be extracted. Based on our discussion above, since the RO-cMF State Mixing Hamiltonian is constructed from only neutral states, we expected the resulting $J$ value to be positive. While this is true when the clusters are single Slater determinants and the only interaction is exchange, the correlation present in the local FCI states seems to provide some additional stabilization to the low-spin states, and the resulting $J$ values are correctly predicted to be negative, albeit small negative values ($-3.8$ $cm^{-1}$ and $-9.9$ $cm^{-1}$ for t-butyl and phenyl substituted phenalenyl dimers, respectively). Perturbative correction stabilizes the low-spin triplet state more than the high-spin singlet state which gives a $J$ value of $-240.1~{}\text{cm}^{-1}$ and $-484~{}\text{cm}^{-1}$ for t-butyl and phenyl substituted phenalenyl dimers, respectively. The spin-states are shown in Fig. (7). The perturbative correction incorporates inter-cluster hopping terms (kinetic exchange) which correctly reveals that radical units in $\pi$-dimers are strongly antiferromagnetically coupled, so much so that they form a weak covalent bond. The difference in exchange coupling constant value arises because of the decreased distance between two carbon atoms in phenyl substituted phenalenyl dimer (3.13 angstrom) than in t-butyl substituted one (3.5 angstrom). Phenyl groups help to delocalize the unpaired electron in phenyl substituted units, and this extended conjugation further seems to help the two radicals come closer to each other that increases the overlap, thus making the $J$ value more negative. However, while the PT2 interactions are significantly stronger than the bare neutral-state interactions, comparison with TPSCI reveals that the PT2 numbers are still almost half as large as they should be (TPSCI yields $-403.2$ cm^-1 for t-Bu substituted phenalenayl dimer and $-782.3$ cm^-1 for phenyl substituted complex).