Restricted Open-Shell Hartree–Fock Method for a General Configuration State Function Featuring Arbitrarily Complex Spin-Couplings

In this work, we present a general spin restricted open-shell Hartree–Fock (ROHF) implementation that is able to generate self-consistent field (SCF) wave functions for an arbitrary configuration state function (CSF). These CSFs can contain an arbitrary number of unpaired electrons in arbitrary spin-couplings. The resulting method is named CSF-ROHF. We demonstrate that starting from the ROHF energy expression, for example, the one given by Edwards and Zerner, it is possible to obtain the values of the ROHF vector-coupling coefficients by setting up an open-shell for each group of consecutive parallel-coupled spins dictated by the unique spin-coupling pattern of any given CSF. To achieve this important and nontrivial goal, we employ the machinery of the iterative configuration expansion configuration interaction (ICE-CI) method, which is able to tackle general CI problems on the basis of spin-adapted CSFs. This development allows for the efficient generation of SCF spin-eigenfunctions for systems with complex spin-coupling patterns, such as polymetallic chains and metal clusters, while maintaining SCF scaling with system size (quadratic or less, depending on the specific algorithm and approximations chosen).


INTRODUCTION
−12 The electronic states of these systems are usually characterized by a complex spin alignment of the open-shell electrons, leading to various ferromagnetic or antiferromagnetic-coupled states.The nature of these states is difficult to interpret from the experimental point of view, due to the high density of manyparticle states that may be probed in the various experimental techniques 13 or the strong dependence of the experimental measurements to the actual experimental conditions, often leading to dependencies (e.g., temperature, field, and radiation) that may be difficult to control or interpret.
The theoretical characterization of such systems is also challenging due to the computational complexity associated with the description of the involved ground and excited manyparticle states.The most common practice is to break the spin symmetry of the system and perform mean-field level computations over chosen broken-symmetry Slater determinants on the basis of a variational orbital optimization (the self-consistent field (SCF)-procedure).−25 However, it should be emphasized that BS-DFT gives only a crude representation of antiferromagnetic couplings, as it is unable to correctly describe the electron correlation phenomena in these systems.
In fact, problems of such complexity require elaborate multiconfigurational wave function-based methods.−31 Recent advances of approximate configuration interaction (CI) wave functionbased methods in the framework of the density matrix renormalization group (DMRG), 32 the full configuration interaction quantum Monte Carlo (FCIQMC), 33−35 or 2. THEORY 2.1.General CSF-Based ROHF Formulation.We start by recalling the ROHF energy expression 1 given by Edwards and Zerner 45

E
(2 ) where h pq is the one-electron integral (possibly containing contributions from external point charges or relativistic corrections), J pq and K pq are the Coulomb and exchange integrals, respectively, n I is the occupation number of shell I, C being the closed-shell, and a IJ and b IJ are the so-called vectorcoupling coefficients.In the following text, we use the indices i, j, k, and l for doubly occupied molecular orbitals, t, u, v, and w for singly occupied molecular orbitals (SOMOs), and p, q, r, and s for generic molecular orbitals.
In contrast to methods like CAHF or SAHF that average over configurations and therefore typically feature fractional occupation numbers, the CSF method constrains the occupations to n I = 1.The spin-coupling situation of the unpaired electrons is determined by the number of open-shells and the vector-coupling coefficients.Hence, for setting up the ROHF problem for a given CSF, we need to determine the values of a IJ and b IJ that appropriately describe it.
We can achieve this by constructing the CSFs using one of the many available construction methods for spin-eigenfunctions, 46 for example, the genealogical spin-coupling scheme of Grabenstetter et al. 47 In this scheme, a given CSF is built by sequentially adding the electrons in such a way that the spincoupling information on an electron u in a SOMO is specified with respect to all other SOMOs t < u.
The expectation value of the Born−Oppenheimer Hamiltonian in the second quantization for a given CSF is then given by eq 2 (2 ) −53 From these, we arrive at the following general formulas for the vector-coupling coefficients 1) This way, the determination of the b IJ vector-coupling coefficient reduces to the evaluation of the matrix element ⟨Φ| E u t E t u |Φ⟩ for the CSF in which the ROHF procedure is to be set.We note here that all t,u orbitals within a given open-shell will produce the same value from the ⟨Φ|E u t E t u |Φ⟩ matrix element.The vector-coupling coefficients within a given shell a II and b II are related only by eq 7 and can be defined in any way that satisfies it. 45t is now possible to set up the Fock operator of the openshells of the system.In the implementation of the general ROHF method available in the ORCA program, we employ the Rico-Fernandez 54

formulation of the Fock matrix elements of the open-shells
2 ( ) ( ) where A μν I and B μν I are special densities The above relations imply that for a given CSF with an arbitrary spin-coupling pattern, it is possible to obtain the values of the ROHF vector-coupling coefficients for a given CSF by setting an open-shell for each group of consecutive parallel-coupled spins.The procedure outlined above was implemented in a development version of the ORCA computational package and will be part of the next major public release.
2.2.Computational Details.−53 Geometries for the [Fe(SCH 3 ) 4 ] − and [Fe 2 S 2 (SCH 3 ) 4 ] 2− complexes are obtained from refs 55 and 56, respectively.For the [Fe 4 S 4 (SCH 3 ) 4 ] 2− cubane system, the geometry used is taken from ref 57.All other geometries were calculated at the DFT level of theory employing the BP86 functional 58,59 together with Grimme's dispersion correction 60−64 with the triple-ζ Def2-TZVP 65 basis set together with the Def2/J auxiliary basis for the resolution of the identity approxima- a Further examples can be found in Table 2.
The Journal of Physical Chemistry A tion. 66 The Ni chain geometry was fixed to local octahedral geometries, 67 with only the hydrogen atoms being optimized at the same level of theory of the other geometry optimizations.
The coordinates of all the geometries used in this paper can be found in the Supporting Information.CI calculations were performed using the CSF-based version of ICE-CI method. 38,39The triple-ζ quality Def2-TZVP 65 basis set was used for all atoms in the systems containing transition metals, for all ROHF (SAHF, CSF-ROHF, and HS-ROHF), ICE-CI, and CASSCF calculations.For the gadolinium dimer, scalar relativistic effects were considered by using X2C, together with the basis set X2C-TZVPall. 68

DEFINITION OF GUESS STARTING ORBITALS
The methodology described in the theory section above allows for determining the ROHF vector-coupling coefficients on an arbitrary number of open-shells that may consist of a given CSF.A second crucial step is to determine which type of starting molecular orbital may be part of an open-shell.This, in principle, will always involve some freedom of choice from the user's perspective.Nevertheless, we draw some general guidelines with the aim to rendering the entire methodology as closely to the black box as possible.
In our experience, the usage of the atomic valence active space (AVAS) 69 procedure constitutes a systematic way to generate metal centered initial guess orbitals suitable for transition metal dimers, trimers, etc.These orbitals are then used as a guess for a SAHF calculation, where the open-shell is defined as the SOMOs of the system.Since the SAHF problem is already set for the average of spin-coupling situations of the target multiplicity, the resulting orbitals should be better suited as a guess for the ROHF calculation on a specific CSF.The converged SAHF orbitals are localized and used as a guess themselves for the CSF-ROHF.Care must be exercised in order for the localized orbitals to reflect the chemical situation of interest, and we point out that, although recommended, not every chemical problem requires the localization step, as can be seen on the [Co( 1 L N ) 2 ] system presented in this study.
One can also recanonicalize the localized orbitals before using them as a guess for the following CSF-ROHF calculation.In ORCA, this is achieved by using the standalone program orca_blockf, which performs a block diagonalization of the Fock matrix for a given subspace.
An example of the protocol used for the setup of the guess orbitals for the CSF-ROHF method is shown in Figure 2 for the [Ni 2 O(H 2 O) 10 ] 2+ system, where we want to set the ROHF problem for two open-shells, each with the singly occupied e g orbitals from a distinct nickel center.
First, a SAHF calculation for 1 open-shell with 4 orbitals and 4 electrons in an overall singlet state was performed using, as an initial guess, the AVAS orbitals.The SAHF converged orbitals were then localized and recanonicalized before being used as the initial guess for the CSF-ROHF calculation, where two open-shells were set containing the respective nickel The Journal of Physical Chemistry A orbitals.As seen in Figure 2, the CSF-ROHF orbitals differ in comparison to the localized SAHF orbitals, reflecting the relaxation of the former for the given CSF.

SCALING AND COMPARISON WITH CASSCF TIMINGS
In the next step, we compare, in this section, the performance of the aforementioned CSF-ROHF procedure.Following the process described in Section 3, for each chain size, guess orbitals were obtained by first converging a SAHF calculation where the open-shell consists of the degenerate set of e g orbitals of the Ni centers to low multiplicity (S = 0 for even numbers of Ni centers and S = 1 for odd numbers).The converged orbitals were then localized using the Pipek− Mezey 70 localization scheme, which were then used as initial guess orbitals for all CSF-ROHF and CASSCF calculations.The CSF-ROHF calculations were set up following the respective branching diagram CSF for the alternating spincoupling (Figure 3).The CASSCF calculations were performed for the lowest root of a minimal active space consisting of the d orbitals of the nickel centers with appropriate multiplicity.The HS-ROHF calculations were set using quasi-restricted orbitals (QROs), 71 obtained from a UHF calculation on the same multiplicity as an initial guess.The usage of QROs usually ensures quick convergence of the HS-ROHF.
The computed average times per SCF iteration and total number of SCF iterations needed to achieve convergence are presented in Figure 4  As shown in Figure 4a, the CASSCF calculation presents the expected steep scaling with the increasing number of Ni(II) centers and becomes prohibitively expensive beyond 4 centers [corresponding to active spaces larger than (32,20)].By contrast, HS-ROHF and CSF-ROHF scale proportionally to N basis 2 if no precautions are taken for the calculation of the Coulomb term.We point out that in CSF-ROHF, a number of Fock matrices proportional to the number of open-shells needs to be processed, while in HS-ROHF, there are always only 2 Fock matrices in the procedure.The presented scaling for the CSF-ROHF method allows its application on problems where CASCF calculations are unfeasible.
It is also worth mentioning the number of SCF iterations needed to achieve convergence.Overall, the CSF-ROHF calculations converge in a reasonable number of iterations, ranging between 40 and 135 iterations, showing that, at least for the systems studied in this work, a rather smooth convergence behavior was observed.While this convergence behavior is, of course, closely related to the specific algorithm used to solve the SCF equations, these results at least indicate that the CSF-ROHF method is not expected to create any convergence issues that are worse than those of the parent HS-ROHF method.

VALIDATION
The validation of the presented CSF-ROHF methodology was performed in a study set that consisted of a series of monometallic and polymetallic complexes belonging to the transition metal and lanthanide families (Figure 5).In particular, the monomer [Fe(SCH 3   = , following Hund's rule, where each of the 3d orbital of the iron center is singly occupied.The CSF representing this spin-coupling situation is [+ + + + +], and its branching diagram representation can be seen in Figure S1. Keeping the occupation number of the five 3d orbitals, it is also possible to construct CSFs for the non-Hund states with S = , representing the aforementioned non-Hund states of [Fe(SCH 3 ) 4 ] − .Each of these CSFs has a unique branching diagram, which allows for the determination of the b IJ vector-coupling coefficients for each case, as shown in Table 2.
The number of open-shells for each CSF is defined by the number of sets of parallel spin-coupling orbitals.As an example, the CSF [+ + − + +] has three groups of parallel spins: the first two +, the following -, and the last two + +, leading to a ROHF problem with 3 open-shells.With the values of the ROHF vector-coupling coefficients determined, the calculation of the total electronic energy for each of the CSF-ROHF cases can be carried out.The total energies are presented in Table 3, as well as the relative energy with respect to the converged HS state.
It should be emphasized that for the case of the unique sextet CSF, both CSF-ROHF and CASSCF (5,5) solve the same problem, since from this CSF within the defined active space, no excitations can be performed.Hence, the CASSCF wave function will also consist of only one sextet CSF and the orbital optimization step reduces to the ROHF problem.This is true for any HS situation where the active space is defined to include only the singly occupied orbitals.

Probing Ferromagnetic and Antiferromagnetic
Couplings.We now turn to the case of the iron dimer [Fe 2 S 2 (SCH 3 ) 2 ] 2− , which consists of two HS Fe(III) centers bridged by sulfides.Assuming two local S = systems, we constructed the ROHF problem for the ferromagnetic coupling (S = 5) and antiferromagnetic coupling (S = 0) situations.For S = 5, only one CSF can be constructed, where the 10 unpaired electrons are all coupled parallel to each other.However, there are 42 CSFs that can be constructed for 10 electrons in 10 orbitals coupled to a singlet state; these can be referred to as the "neutral" CSFs since both Fe centers have the same charge.Among these 42, we selected the CSF with [+ + + + + − − − −] coupling as the illustrative example of the antiferromagnetic coupling between the two Fe centers.
The CSF-ROHF problem was then set up in the following way.For the HS S = 5 state, we have 1 open-shell with 10 orbitals and 10 electrons.For the LS S = 0 state, two openshells are employed, each containing 5 electrons and the 3d orbitals of the Fe(III) centers.For comparison, CASSCF calculations were also performed on both spin states using the same set of guess orbitals.The energy differences between both spin states are shown in Table 4.
Inspection of Table 4 shows that, in a first approximation, there is a disagreement between the CSF-ROHF and CASSCF energy differences between the two spin states, with the CSF-ROHF energy difference favoring the S = 5 state.The reasoning behind this discrepancy follows from the mechanism of antiferromagnetic coupling between two metal centers.As discussed in more detail elsewhere, 56,72−74 in order to properly describe the coupling between the two Fe centers, not only the "neutral" CSFs are needed but also the "ionic" CSFs, in which one electron is transferred from one iron center to the other.Since the ROHF calculation was performed exclusively for one of the "neutral" CSFs, it fails to capture the antiferromagnetic coupling in its entirety.Hence, in a subsequent step, the missing "ionic" CSFs can be included on the basis of CAS-ICE calculations with the ROHF-optimized orbitals (Figure 6).This lowers the energy of the S = 0 state by 2084 cm −1 , due to mixing with the "ionic" CSFs, providing consistent energy splittings between CASSCF and CSF-ROHF/CAS-ICE.We emphasize that even though there is an energy lowering resulting from the mixing of other CSFs, the wave function  5).For comparison, CAS-ICE calculations were also performed using the HS ROHF solution, where the obtained orbitals were just localized before the CI.In this case, the resulting wave function is no longer dominated by the [+ + + + + − − − − −] CSF, showing the inadequacy of the HS solution as a basis for subsequent correlated calculations of the lower multiplicities (consistent with ref 35).
Likewise to the monomer case, the mixing of CSFs discussed above is not present in the HS state, which is uniquely defined by a single CSF.This renders the CASSCF calculation to be exactly the ROHF problem, leading to the same final electronic energy obtained by both methods (Table S1).
The same procedure used for the iron dimer is not restricted only to transition metals, but it is generally applicable to any situation where the interacting open-shells can be defined.As an example, we used a model Gd(III) dimer [Gd 2 Cl 11 ] 5− where the gadolinium centers have f 7 electronic configuration in a local S = spin state.Once again, by employing CSF-ROHF, CASSCF, and CAS-ICE calculations, both ferromagnetic (S = 7) and antiferromagnetic (S = 0) coupling situations may be probed (Table 4).In contrast to [Fe 2 S 2 (SCH 3 ) 2 ] 2− , in the case of [Gd 2 Cl 11 ] 5− , the energy difference between the ferromagnetic and antiferromagnetic coupling situations is negligible with all methods applied.This is due to the lack of covalency between metal and ligands, causing the two gadolinium centers to be essentially isolated from each other, as is also reflected in the ROHF canonical MOs obtained for both spin-coupling situations (see Supporting Information).
Going up in complexity, we address the model iron−sulfur cluster [Fe 4 S 4 (SCH 3 ) 4 ] 2− (Figure 5).It is generally accepted 57 that the spin situation in this cluster involves mixed valence Fe(II)−Fe(III) pairs, also referred as Fe 2.5+ −Fe 2.5+ , that couple ferromagnetically within themselves and antiferromagnetically with the pair on the opposite face of the cube, resulting in an overall singlet state.
With the accepted spin situation in mind, the open-shells of the system are to be specified using a set of 18 orbitals and 18 electrons resulting from 2 Fe(II) and 2 Fe(III) atoms.This gives rise to the question of how to define the open-shells for this system.Although we can use the CSF represented on the branching diagram of Figure 7, there is still ambiguity on which orbitals to include on the two open-shells.Since there are 4 iron centers that can be paired and the optimized structure of this model is not a perfect cube, there are 3 distinct possible iron pairings to consider.
Following the procedure described in Section 3, by first performing a SAHF calculation with just one open-shell containing all iron atoms and an overall singlet state, we obtain converged orbitals that, after localization with the Pipek− Mezey scheme, 70 are already paired in a distinct situation.These orbitals can be used as the guess orbitals on the CSF-ROHF problem for S = 0 CSF, where the two open-shells are defined according to the pairing of SAHF orbitals.For comparison, we also performed the conventional HS ROHF  The Journal of Physical Chemistry A for the S = 9 state.The ROHF energy differences between the two spin-coupling situations are listed in Table 4. Again, the correct description of the antiferromagnetic coupling can be obtained only after the inclusion of the missing "ionic" CSFs in a subsequent CI calculation for the S = 0 problem.In line with the [Fe 2 S 2 (SCH 3 ) 2 ] 2− dimer results, after CI, the wave function is still dominated by the CSF [+ + + + + + + + + − − − − − − − − −] (Table 5), this time for both CI calculation employing MOs obtained from the two spincoupling ROHF problems.
As has been already shown, 35,57 probing the actual magnetic structure of the ground state of these systems requires expansion of the active space to include all relevant metal and ligand-based orbitals participating in metal−ligand covalent interactions.In the context of approximate CI, further CSF mixing with the "neutral" CSF will be introduced, which, in principle, would require further orbital relaxation than the CSF-ROHF provides before robust results could be achieved.This is beyond the scope of the current work.

Probing Complex Spin-Coupling Situations.
In the previous section, we encountered the problem of defining the open-shells in spin-coupling situations where this choice was somewhat straightforward.However, this is not always the case, particularly in multimetallic systems with more than two metal centers.A characteristic case refers to the copper trimer [Cu 3 (OH) 3 (en) 3 ] 3+ .In this system, there are 3 Cu(II) centers each with local S 1 2 = .In this scenario, the three centers can couple to an overall quartet state or to a doublet.
For the quartet state, there is only one possible CSF that can be constructed, where all unpaired spins are coupled parallel to each other.In fact, for this case, the open-shell definition is unique since there is only one open-shell with 3 electrons and 3 orbitals.As can be seen in Table 6, the three possibilities result in similar single point energies, in accordance with the local ∼ C 3 symmetry of the Cu centers of the complex (Figure 8  The CSF-ROHF optimized orbitals are shown in Figure 9, and the calculated ROHF energies are presented in Table 7.One can readily see that all ROHF solutions for the CSFs with 3 orbitals and 3 electrons are lower in energy than the closedshell doublet.We also performed CAS-ICE calculations using different ROHF solutions in order to obtain a more correct description of the spin-coupling situation of the system.

The Journal of Physical Chemistry A
As can be seen in Table 8, after the CI calculation using either the ROHF orbitals obtained for the [+ − +] or [+ + −] CSFs, the wave function is over 70% composed of the [+ + −] one.This, together with the energies ca.847 cm −1 lower than the quartet state reflects the diradical character of the ligands, in accordance with the accepted electronic structure. 75The difference in the percentage composition for these two CI calculations reflects the two CSF-ROHF solutions for the [+ − +] and [+ + −] CSFs.

CONCLUSIONS
In this study, we presented a method to calculate the ROHF vector-coupling coefficients for a given CSF, which allows for the calculation of proper spin-eigenfunctions of arbitrary spincoupling situations.The method is generally applicable to systems where one defines an open-shell consisting of only singly occupied orbitals, meaning that the occupation number n μ for the open-shell part of the restricted open-shell Fock operator is always 1.This allows the formulation of a general spin ROHF methodology starting from a given CSF that can be used to provide a SCF solution for a given spin-coupling situation.
CSF-ROHF makes use of an efficient infrastructure to compute the vector-coupling vectors, by computing the needed two-body ⟨Φ|E u t E t u |Φ⟩ matrix elements directly in the CSF basis of the problem, with no expansion in Slater determinants required.It scales close to quadratically with system size and shows satisfactory convergence performance.
Ongoing efforts are directed toward the development of dynamic correlation methods that start from the CSF-ROHF reference states.Next to many-body perturbation and coupled cluster approaches, obvious choices are methods related to approximate full-CI schemes such as DMRG, 32 FCI-QMC, 33−35 or the CIPSI 36,37 /ICE 38,39 family of selecting CI methods.We foresee exciting applications of such methods to challenging chemical problems in (bio)chemistry, catalysis, and material sciences.

The Journal of Physical Chemistry A
expression for a IJ and b IJ , the relations, eqs 4 D μνJ is the density matrix element for shell J α IJ and β IJ are related to a IJ and b IJ by

2 =
As a numerical example, let us consider the cases of CSFs consisting of (a) two unpaired electrons coupling to an opensinglet [+ −] with S = 0, and (b) three unpaired electrons coupling to open-doublets (S 1 ), where two CSFs are possible: [+ − +] and [+ + −], as shown in Figure 1.As stated above, the open-shells are defined as each set of consecutive parallel-coupled spin.Therefore, the open-singlet [+ −] CSF contains two open-shells and the open-doublets [+ − +] and [+ + −] contain three and two open-shells, respectively.The values obtained for the open-singlet and the two opendoublet cases are presented in

Figure 2 .
Figure 2. Example sequence of molecular orbital transformations used in the protocol for spin-couplings between metal centers.From left to right: the AVAS orbitals; the SAHF orbitals obtained using the AVAS guess; the localized orbitals obtained from the canonical SAHF ones; the recanonicalized orbitals obtained by using the orca_blockf utility; and the CSF-ROHF-optimized orbitals obtained from the localized and recanonicalized SAHF orbitals.
For this purpose, we perform a series of calculations on [Ni(H 2 O)6] 2+ as well as the [Ni(H 2 O) 5 ] n O n−1 , n = 2−10 chain, where the number of nickel atoms is gradually increased from 2 to 10 Ni atoms (from 303 to 2535 basis functions).Each nickel added is coupled antiferromagnetically with the previous one, resulting in a series of CSFs following the spin-coupling scheme, as shown in Figure 3.
as a function of the number of Ni centers in [Ni(H 2 O)6] 2+ and [Ni(H 2 O) 5 ] n O n−1 , n = 2−10.All comparisons were performed in an Intel cluster using 8 cores and 4 GB of RAM per core.

Figure 3 .
Figure 3. Model Ni chain used for timings and the branching diagram representation of the alternating spin-coupling between the Ni centers.Dashed line represents a continuation of the diagram as more Ni centers are added.

Figure 4 . 5 . 1 .
Figure 4. a) Average SCF iteration time for CSF-ROHF, HS-ROHF, and CASSCF calculations on the Ni chain with an increasing number of Ni atoms.(b) Number of SCF cycles needed to achieve convergence of the respective methods.

Figure 5 .
Figure 5. Structures of all systems studied in the present paper.

The Journal of Physical Chemistry A
Being able to directly evaluate the ⟨Φ|E u t E t u |Φ⟩ element, the value of the vector-coupling coefficients can be obtained as follows: By comparing the terms involving only the open shells of eq 2 with the open-shell terms of the ROHF energy expression 1 for n I = 1, one arrives on eq 3 https://doi.org/10.1021/acs.jpca.4c00688J. Phys.Chem.A 2024, 128, 5041−5053 5042

Table 1 .
We note that the values of the [+ −] and [+ − +] CSFs are in agreement with the ones determined by Edwards and Zerner, which are presented on lines 6 and 8 of Table 1 of ref 45, together with the b II values of the given cases.

Table 2
also shows the respective numbers of open-shells for each CSF studied.

Table 2 .
Number of Open-Shells and Determined Values for b IJ

Table 7 .
Calculated Energies for the Three CSFs of Complex [Co( 1 L N ) 2 ] and the Energy Difference in Respect to the One Open-Shell Doublet State

Table 8 .
CAS-ICE Wavefunction Composition UsingDifferent CSF-ROHF Orbitals (Figure9) in Terms of the Possible CSFs for 3 Electrons in 3 Orbitals a aA number 2 represents a doubly occupied orbital and a number 0 represents a virtual orbital.