Minimal Active Space for Diradicals Using Multistate Density Functional Theory

This work explores the electronic structure as well as the reactivity of singlet diradicals, making use of multistate density functional theory (MSDFT). In particular, we show that a minimal active space of two electrons in two orbitals is adequate to treat the relative energies of the singlet and triplet adiabatic ground state as well as the first singlet excited state in many cases. This is plausible because dynamic correlation is included in the first place in the optimization of orbitals in each determinant state via block-localized Kohn–Sham density functional theory. In addition, molecular fragment, i.e., block-localized Kohn–Sham orbitals, are optimized separately for each determinant, providing a variational diabatic representation of valence bond-like states, which are subsequently used in nonorthogonal state interactions (NOSIs). The computational procedure and its performance are illustrated on some prototypical diradical species. It is shown that NOSI calculations in MSDFT can be used to model bond dissociation and hydrogen-atom transfer reactions, employing a minimal number of configuration state functions as the basis states. For p- and s-types of diradicals, the closed-shell diradicals are found to be more reactive than the open-shell ones due to a larger diabatic coupling with the final product state. Such a diabatic representation may be useful to define reaction coordinates for electron transfer, proton transfer and coupled electron and proton transfer reactions in condensed-phase simulations.


Introduction
"It is hard to imagine complex chemistry without diradicals", said Hoffmann and coworkers [1]. Diradicals are reactive species and have been extensively studied both theoretically and experimentally [1][2][3][4][5]. The spin coupling between two unpaired electrons results in a singlet state and a low-energy triplet state [6][7][8][9][10]. Both states, having an M S = 0, are of multiconfigurational character that cannot be adequately treated by Kohn-Sham density functional theory (KS-DFT) [6][7][8]. However, the energy-degenerate M S = ±1 spin multiplets of the triplet state can be represented by a single determinant using KS-DFT [6,7]. Furthermore, electron correlation plays an important role between the two closed-shell configurations in a two-electron and two-orbital active space, which ultimately determines the relative energies of the singlet and triplet ground states [9,10]. Consequently, the need to use a multiconfigurational method to treat open-shell diradicals poses a significant challenge to KS-DFT, and often a broken-symmetry approach is used to estimate the singlet-triplet-energy gap [11][12][13]. On the other hand, the difficulties of KS-DFT can be easily overcome using multiconfiguration self-consistent-field (MCSCF) methods in wave function theory (WFT), such as the complete-active-space self-consistent-field (CASSCF) approach, but, in this case, it is necessary to use a large active space, followed by including

Theoretical Background
The matrix element of the Hamiltonian matrix functional H[D(r)] of the multistate density D(r) in the subspace spanned by N states, R N ∈H, is given by [18,20,23] where the terms on the right-hand side of the equation are, respectively, the kinetic energy, Coulomb (Hartree) and exchange energy, the external potential energy, and the exchangecorrelation energy for the interactions between states A and B. If A = B, Equation (1) is equivalent to KS-DFT for the density D AA (r) = ρ A (r) represented by the determinant Φ A . However, for A = B, each term is a functional of the transition density D(r). Both state and transition densities are related to the one-particle density matrix P AB (below); it is important to emphasize that the matrix density D(r) is not to be confused with P AB , which is not a function of r.
D AB (r) = m ∑ µν χ µ (r) > (P AB ) µν < χ ν (r)| = χ T (r)P AB χ(r), where χ(r) is a column vector of m atomic orbital basis functions, and the one-particle density matrix P AA (A = B) and transition density matrix P AB (A = B) for determinants Φ A and Φ B are related to the coefficient matrices C A and C B of occupied orbitals by with R being the overlap matrix of the basis functions [18,20]. For the molecular systems considered in this work, the external potential in Equation (1) is simply the nucleus-electron Coulomb energy, and the first three terms can be expressed together by where h and G are standard matrices of one-electron integrals (kinetic and nuclear attraction integrals) and two-electron Coulomb-exchange integrals. For the last term of Equation (1), Kohn-Sham exchange-correlation functional can be used to approximate the diagonal element of the multistate matrix functional, E xc AA [D AA (r)] = E KS xc [D AA (r)], but the offdiagonal elements are the transition density correlation functional (TDF), which does not exist in KS-DFT [20]. In special situations, such as the present singlet-triplet-state energy difference involving spin-pairing interactions between two unpaired electrons, the TDF energy can be obtained by enforcing the spin-multiplet degeneracy condition of the triplet states [24][25][26]29].
In Equation (5), ] are correlation energies computed using KS-DFT correlation functional for the triplet state and spin-contaminated configuration. The dynamic correlation energy contribution to the electronic coupling, E xc AB [D AB , Φ ↑↓ A , Φ ↓↑ B ] between two spin-mixed determinants treated by KS-DFT ensures that the resulting triplet state with M S = 0 is degenerate with the M S = +1 component [25,26]. In turn, the pure singlet state energy is also determined. We note that the M S = +1 component of a triplet state can be adequately treated by KS-DFT with a single Slater determinant. Since the KS exchangecorrelation functional leads to the exact ground-state energy, the TDF energy for diabatic electronic coupling in Equation (5) is exact, even though the specific functional form of TDF is unknown.
For non-spin-coupled determinants, we use the KS-DFT energy-weighted correlation energy to approximate the correlation energy of the TDF: where H NO AB =< Φ A H Φ B > is the nonorthogonal matrix element of the determinants Φ A and Φ B , and the energies in the denominator are Hartree-Fock energies using KS orbitals.
The adiabatic eigenstate energies are obtained by diagonalizing the Hamiltonian matrix function H[D(r)]: where I denotes the Ith adiabatic eigenstate, and A and B specify determinant functions {Φ A } representing the matrix density D(r). Equation (7) shows that E MS I [D(r)] is an implicit functional of the eigenstate density ρ a AI a BI D AB (r), which is not directly used as an input in a Kohn-Sham density functional approximation; a KS functional is optimized corresponding to the residual kinetic energy from a single determinant.

Computational Details
Molecular geometries were optimized using the Gaussian09 software with the M06-HF functional (guess = mix option was used in open-shell calculations) [38]. For comparison, we also performed CASSCF and CASPT2 calculations using the MolPro2012 package [39]. Dunning's correlation consistent valence triple-zeta basis set (cc-pVTZ) was used unless specifically noted in the text [40].
In MSDFT calculations, the diagonal matrix elements of the Hamiltonian matrix functional is determined either by block-localized KS-DFT, in which effective valence bond-like configurations are represented by fragment block-localized orbitals [18], or by targeted optimization using a delta-SCF approach with a pre-selection transformation of the target orbitals [41]. For the off-diagonal elements, the transition density correlation functional (TDF) was obtained by applying the spin multiplet degeneracy constraint for spin coupling interactions (Equation (5)), or by a determinant-coupling weighted correlation energy for other situations (Equation (6)). The latter scales similarly to that with an overlap weight [20,23]. All determinant states were optimized using unrestricted KS-DFT, constrained either by fragment block-localization [23] or by a non-aufbau configuration with specific target orbital occupations [41], i.e., delta-SCF optimization. Additional computational details are provided in the text and in Supporting Information. MSDFT calculations were performed using a program interfaced with a locally modified GAMESS program for integral evaluations [42].

Potential Energy Curves of Hydrogen Molecule
We first illustrated the computational procedure of the MSDFT-NOSI method with a minimal active space (MAS) by computing the potential energy curves for the lowest 4 states in the M S = 0 manifold (three singlet and one triplet) of a hydrogen molecule (H 2 ) [1,43,44]. The minimal representation of these four states includes two doubly occupied closedshell (CS) configurations, and a pair of singly occupied open-shell (OS) configurations by placing two electrons in the two lowest-energy orbitals. Thus, the state densities of the subspace R 4 for the lowest four states of H 2 are represented by the following configuration state functions: whereÂ is the antisymmetry operator, φ  (8) and (9), χ H and χ L are, respectively, the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). Note that although the same symbols χ H and χ L are used to denote Kohn-Sham orbitals, they are generally different in different determinants and are nonorthogonal.
Because of spin symmetry, the triplet state Θ T OS (Equation (11)) is orthogonal to the singlet states (Equations (8)- (10) where Equations (12)- (15) are the four lowest adiabatic states of H 2 , and a y x are the coefficients of CSFs. Note that all adiabatic states are of a multiconfigurational character, including the closed-shell states, although Θ 20 CS makes the predominant contribution to the ground state near bonding distances. Figure 1 illustrates that the ground-state potential energy curve of H 2 along the bond-distance coordinate can be roughly divided into three regions, corresponding to closed-shell (I), diradicaloid (II) and diradical (III) [1]. In region I, at the bonding distance, the energy of the Θ 20 CS state is below the T 1 state (Ψ T 1u ), contributing dominantly to the singlet ground state Ψ S 0 1u ( Figure 2a). In region II, the Ψ S 0 1u singlet ground state is characterized as a diradicaloid [1], which has an increasing amount of multireference character with a large and rapidly increasing diradical index in Figure 2a [9]. Finally, in region III, the Θ 20 CS = φ 20 CS configuration (Equation (8)) has a higher energy than that of the triplet state ( Figure 1), but Ψ S 0 1u is still below the triplet state as a diradical species. In this region, the singlet and triplet ground states are nearly degenerate (degenerate at infinite separation). Here, the singlet state can be equivalently viewed as a multiconfiguration combination of two CS configurations (Equation (12)) and two localized open-shell diradicals as in valence bond (VB) representation. Here, the triplet state results from the combination of two spin-mixed OS determinants (Equation (15)). Nevertheless, the geometry at which the energies of the closed-shell CSF Θ 20 CS and the triplet-state Ψ T 1u switch order may be used as the transition point from the closed-shell diradicaloid bonding character to a singlet diradical species.     MSDFT-NOSI calculations along with the FCI results (dotted curves). The agreement between MSDFT and FCI results is good, although the relative energy between the singlet and triplet states from MSDFT-NOSI is somewhat greater than that of FCI, suggesting that electronic coupling between the two CS configurations may be overestimated. Significantly, the correct trend of the potential energy curve for H-H dissociation is obtained by MSDFT-NOSI through state interaction between the incorrect dissociation behavior of the   (11)), and the lowest singlet Ψ S 0 1u and triplet Ψ T 1u adiabatic states for H 2 from MSDFT-NOSI calculations along with the FCI results (dotted curves). The agreement between MSDFT and FCI results is good, although the relative energy between the singlet and triplet states from MSDFT-NOSI is somewhat greater than that of FCI, suggesting that electronic coupling between the two CS configurations may be overestimated. Significantly, the correct trend of the potential energy curve for H-H dissociation is obtained by MSDFT-NOSI through state interaction between the incorrect dissociation behavior of the spin-mixed determinants (Θ 20 CS and Θ 02 CS ). The area in blue between 1.45 Å and 2.2 Å in Figure 2b corresponds to region II in Figure 2a, which extends to longer bond distances than the latter (CASSCF) calculations without including the dynamic correlation. Figure 2 highlights the gradual transition in bonding character from a closed-shell configuration to a diradical state. The crossing point between the Θ 20 CS diabatic state and the triplet state Ψ T 1u is located at 1.98 Å from MSDFT, consistent with the earlier studies [1,46].

Singlet-Triplet-Energy Gap
The singlet-triplet-energy splitting, ∆E ST = E S − E T , is an important property of diradicals [47][48][49][50]. Here, we summarize the results for three sets of compounds, including benzyne isomers, cyclobutadiene and polyacenes. Table 1 lists the computed singlet-triplet-energy gaps for the three benzyne isomers (meta, ortho and para) using DFT and WFT approaches [51] along with experimental values [52]. We found that methods that included both static and dynamic correlations clearly outperformed those that lacked either dynamic (CASSCF) or static correlation (UDFT) in comparison with experimental data. MSDFT employing the M06-HF functional and SF-CCSD, both of which included dynamic and strong correlations, had the lowest mean unsigned errors (MUEs), ranging from 0.5 to 0.7 kcal/mol. CASPT2 also yielded reasonable results with a somewhat larger MUE (1.3 kcal/mol) than MSDFT and SF-CCSD approaches. CASSCF systematically underestimated ∆E ST for this series of compounds, having the largest unsigned errors (6.8 kcal/mol) among multiconfigurational methods. Interestingly, the use of unrestricted orbitals (UDFT) produced results of a similar quality as that of CASSCF calculations. It is interesting to note that the energies determined for o-benzyne were the same for RDFT and UDFT, strongly suggesting that a closed-shell configuration was the singlet ground state of this molecule [52]. On the other hand, pand m-benzyne showed significant spin-contamination, evidenced by the computed S 2 values (Table S5) and the large difference between ∆E ST from UDFT (underestimate) and RDFT (overestimate). Clearly, static correlation plays an important role to adequately describe the singlet diradicals.

Cyclobutadiene
The change in ∆E ST versus the geometry interconversion between the two rectangular structures were examined (Figure 3). The transition-state structure at the squared structure at D 4h symmetry on the singlet-state potential energy surface was a minimum of the triplet state [53]. At this geometry, the lowest singlet state corresponds to a combination of two doubly occupied (closed-shell) configurations [54,55], which can undergo Jahn-Teller distortion, leading to two rectangular geometries with D 2h symmetries [56]. The potential energy profiles determined with the CASSCF, multi-reference configuration interaction (MRCI) and MSDFT methods are shown in Figure 3 for the lowest singlet and triplet states. In Figure 3, the reaction coordinates for the interconversion between the two rectangular structures is interpolated to match the difference between the two rectangular sides (bond lengths). At the square geometry (D 4h ), the interpolation reaction coordinate is 0, and the two rectangular structures correspond to a unitless value of ±5. We used a minimum active space of two electrons in two orbitals (Equations (8)-(11)) both for MSDFT and CASSCF, whereas an active space of 4 electrons in 4 orbitals was also employed in the multireference MRCI calculations. However, we emphasized that a common set of orbitals were used both in CASSCF and in MRCI calculations, but nonorthogonal orbitals were adopted in the present MSDFT-NOSI. Figure 3 shows that, at the optimal geometries of the rectangular minima, all methods produce the correct order in energy where the singlet state is below the triplet state. Interestingly, triplet-state energies are essentially the same for all methods, but there is significant variation in the energy of the singlet state among the different methods. In particular, the singlet state is above the triple at the D 4h geometry by CASSCF (2,2). The correct ∆E ST difference is obtained from MRCI and MSDFT, both using a minimal active space of two electrons in two orbitals. This is confirmed further by the comparison with MRCI calculations using a larger active space of four electrons in four orbitals [57]. Interestingly, the correct ∆E ST difference at the square geometry can also be achieved by expanding the active space to four electrons in four orbitals, even without perturbation correction ( Figure 3). Clearly, the amount of dynamic correlation introduced by expanding the size of the active space is sufficient to recover the correct energy difference. The effect introduced by using a larger active space has been called spin-dynamic correlation in the literature [58][59][60]. Of all the methods examined, MSDFT-NOSI, employing M06-HF/cc-pVTZ, yields a somewhat large ∆E ST value at 23.2 kcal/mol. SF-EOM-CCSD results are found at 16 kcal/mol [57]. In closing this discussion, we note that Jahn-Teller effects are found in a variety of systems, for example, recent computational studies of gold nanoclusters that include 25 gold atoms in different oxidation states [61], a classical example of the compressed and elongated conformers resulting from ionization of benzene [20,25,62], and coronene radical cation [63]. interaction (MRCI) and MSDFT methods are shown in Figure 3 for the lowest sin triplet states. In Figure 3, the reaction coordinates for the interconversion between rectangular structures is interpolated to match the difference between the two rec sides (bond lengths). At the square geometry (D4h), the interpolation reaction co is 0, and the two rectangular structures correspond to a unitless value of ±5. W minimum active space of two electrons in two orbitals (Equations (8)-(11)) MSDFT and CASSCF, whereas an active space of 4 electrons in 4 orbitals was ployed in the multireference MRCI calculations. However, we emphasized that a c set of orbitals were used both in CASSCF and in MRCI calculations, but nonort orbitals were adopted in the present MSDFT-NOSI. Figure 3. Computed potential energy curves for the ground-state singlet S0 and triplet T1 cyclobutadiene along the interpolated reaction coordinate between the two rectangular (D etries at a unitless value of ±5 through the transition state geometry (D4h). The cc-pVTZ bas used in all calculations. The energy of triplet state at the square geometry was set as the r since it has the least multiconfigurational character. The M06-HF functional was used in calculations. Figure 3 shows that, at the optimal geometries of the rectangular minima, all m produce the correct order in energy where the singlet state is below the triplet sta estingly, triplet-state energies are essentially the same for all methods, but there icant variation in the energy of the singlet state among the different methods. In pa the singlet state is above the triple at the D4h geometry by CASSCF (2,2). The corre difference is obtained from MRCI and MSDFT, both using a minimal active spac electrons in two orbitals. This is confirmed further by the comparison with MRCI tions using a larger active space of four electrons in four orbitals [57]. Interestin correct ST E  difference at the square geometry can also be achieved by expan Figure 3. Computed potential energy curves for the ground-state singlet S 0 and triplet T 1 states of cyclobutadiene along the interpolated reaction coordinate between the two rectangular (D 2h ) geometries at a unitless value of ±5 through the transition state geometry (D 4h ). The cc-pVTZ basis set was used in all calculations. The energy of triplet state at the square geometry was set as the reference, since it has the least multiconfigurational character. The M06-HF functional was used in MSDFT calculations.

Polyacenes
The singlet-triplet-energy gaps in the series of polyacenes have been extensively studied experimentally and computationally [17,[64][65][66][67][68]. Since there is a large number of calculations in the literature, we chose to use the work by Yang and coworkers for a comparison with the present results [67]. Shown in Figure 4 are the computed ∆E ST in the range of n = 3 to 15 rings, using both restricted and unrestricted KS-DFT (RDFT and UDFT) and MSDFT, employing the M06-HF functional and 6-31G(d) basis set [69], along with the computational results obtained by Yang and coworkers [67]. RDFT results incorrectly show an ∆E ST energy inversion between singlet and triplet states at n = 8, but the values both from UDFT and from MSDFT calculations indicate that the energies of the singlet state are uniformly below those of the triplet state in this series of compounds. This finding is consistent with a previous study that employed the density matrix renormalization group (DMRG) method, where the ground state remains to be singlet as the chain length increases with a finite singlet-triplet gap in the infinite chain limit (2-12 kcal/mol) [1,17]. The computed ∆E ST from UDFT exhibits a maximum at n = 9, in contrast to MSDFT results and DMRG studies. Interestingly, the particle-particle random phase approximation (pp-RPA) results seem to follow the UDFT trend [67], although at a much slower rate ( Figure 4). Nevertheless, the difference in ∆E ST propagation between UDFT and RDFT calculations may indicate that there is a transition from a predominantly closed-shell diradicaloid configuration to a diradical state as the number of fused rings increases beyond 8. Figure 4 illustrates that MSDFT-NOSI with a minimal active space of [2,2] nonorthogonal configurations can adequately balance static and dynamic correlations through this transition.

Hydrogen-Atom Transfer Reactions
The hydrogen-atom transfer (HAT) reactions of SiH4 by the singlet diradicals of cyclobutadiene and para-benzyne may be formally considered as a concerted electron-proton transfer (CEPT) process [1], and we previously introduced a diabatic state representation of the different natures in electronic structure between HAT and CEPT using MSDFT [30]. Since these two reactions are known to be hydrogen atom abstraction, we restricted our discussion to the HAT diabatic states. Here where the molecular fragment blocks of the diradical (DR) and the substrate (SiH4) in the reactant state and those of the HAT product (DRH) and the SiH3 free radicals (Scheme 1)

Hydrogen-Atom Transfer Reactions
The hydrogen-atom transfer (HAT) reactions of SiH 4 by the singlet diradicals of cyclobutadiene and para-benzyne may be formally considered as a concerted electron-proton transfer (CEPT) process [1], and we previously introduced a diabatic state representation of the different natures in electronic structure between HAT and CEPT using MSDFT [30]. Since these two reactions are known to be hydrogen atom abstraction, we restricted our discussion to the HAT diabatic states.
Here, we used a minimal active space of three spin-adapted diabatic states, two for the reactant state corresponding to the closed-shell configuration Θ R CS and the open-shell state Θ R OS of the diradical species, and one for the product state Θ P HAT consisting of two separate radical species that are spin coupled ( Figure 5). Each of the diabatic states are expressed in terms of two fragment block-localized determinants where the molecular fragment blocks of the diradical (DR) and the substrate (SiH 4 ) in the reactant state and those of the HAT product (DRH) and the SiH 3 free radicals (Scheme 1) are grouped in parentheses, the superscript core denotes a product of doubly occupied orbitals that do not change occupation during the reaction, the superscripts α and β specify the electron spin, the orbitals χ H and χ L are the highest-occupied and lowest-unoccupied block-localized KS (BLKS) orbitals of the diradical fragment, and γ hyd represents the 1 s orbital of the hydrogen atom (it is a BLKS orbital of Si-H bond in the reactant state) [30]. represents the 1 s orbital of the hydrogen atom (it is a BLKS orbital of Si-H bond in the reactant state) [30].

The Hamiltonian matrix functional [ ( )]
H D r is expressed in terms of the BLKS de terminant functions as a 6 × 6 matrix. Its elements are determined as follows, keeping in mind that the determinant states are individually optimized to yield the variational dia batic states (VDS) that are best valence-bond-like representations of these Lewis structure [70][71][72][73].


The diagonal elements of [ ( )] H D r are directly determined as the BLKS-DFT ener gies of the corresponding determinants. The two reactant diabatic states can be sep arately obtained as the lower root of a 2 × 2 NOSI diagonalization of the two states in Equations (16) and (17) for illustration in Figure 6, although it is not needed to deter mine the potential energy surface of the adiabatic ground state. The four determi nants in Equations (16) and (17) involve block-local excitations, for which the optimi zation has been detailed in reference [41]. For the product state, only ground-state BLKS optimization is sufficient;  The spin-coupling matrix element between the spin-pair determinants in Equation (18) is evaluated using Equations (1) and (5), which requires a separate BLKS calcu lation of the triplet state with MS = +1;


To examine the variations of state interactions as the HAT occurs, we also computed the effective diabatic coupling values between the reactant (Equations (16) and (17)) and product (Equation (18) , where R = 1 and 2, P = 3, and g  is the adiabatic ground-state energy. The two processes in Scheme 1 correspond to hydrogen-atom abstractions by a πtype diradical (cyclobutadiene) and a σ-type diradical (p-benzyne), and the former has been previously described [1], providing a reference of data for comparison with the present MAS in MSDFT-NOSI calculations.
Shown in Figure 6  . In addition, the effective reactant-product diabatic coupling terms are displayed. Evidently, the open-shell diabatic state (Equation (17)) has the highest energy in both reactions (Scheme 1), and its effective coupling (V23) with the product state (dotted green curves) is weak with little variations along the reaction coordinate ( Figure 6). On the other hand, interactions between the CS reactant state and product state (V13) is strong, reaching its maxima near the transition states for the two HAT reactions. The substantial difference in the effective diabatic couplings between Θ and Θ diabatic states indicates that the reactivities of both the π and σ types of biradicals in the singlet states are dominantly determined by the closed-shell states. This result is in good accord with previous analyses by Hoffman and coworkers who found that both singlet diradicals are best characterized by two closed-shell configurations [1].

•
The spin-coupling matrix element between the spin-pair determinants in Equation (18) is evaluated using Equations (1) and (5) (6)); • To examine the variations of state interactions as the HAT occurs, we also computed the effective diabatic coupling values between the reactant (Equations (16) and (17)) and product (Equation (18)) states, denoted as V 13 and V 23 according to V RP = H RP − ε g S RP , where R = 1 and 2, P = 3, and ε g is the adiabatic ground-state energy.
The two processes in Scheme 1 correspond to hydrogen-atom abstractions by a π-type diradical (cyclobutadiene) and a σ-type diradical (p-benzyne), and the former has been previously described [1], providing a reference of data for comparison with the present MAS in MSDFT-NOSI calculations.
Shown in Figure 6 are the reaction energy profiles of the adiabatic ground state (S ad 0 ) and the diabatic states for the two reactant states, Θ R CS and Θ R OS , along with the HAT state Θ P HAT . In addition, the effective reactant-product diabatic coupling terms are displayed. Evidently, the open-shell diabatic state (Equation (17)) has the highest energy in both reactions (Scheme 1), and its effective coupling (V 23 ) with the product state (dotted green curves) is weak with little variations along the reaction coordinate ( Figure 6). On the other hand, interactions between the CS reactant state and product state (V 13 ) is strong, reaching its maxima near the transition states for the two HAT reactions. The substantial difference in the effective diabatic couplings between Θ R CS and Θ R OS diabatic states indicates that the reactivities of both the π and σ types of biradicals in the singlet states are dominantly determined by the closed-shell states. This result is in good accord with previous analyses by Hoffman and coworkers who found that both singlet diradicals are best characterized by two closed-shell configurations [1].
For comparison, the computed energy barrier for the hydrogen-atom abstraction of SiH 4 by cyclobutadiene is 12 kcal/mol MSDFT-NOSI employing the M06-HF functional and the cc-pVTZ basis set ( Figure 6), which is in reasonable accord with a value of 15 kcal/mol using CCSD(T)/cc-pVTZ//B3LYP/cc-pVTZ [1]. Employing the same structures reported by Hoffmann and coworkers [1], we found that the difference in the reaction energy between the two methods is greater (2 vs. −2 kcal/mol) [1]. There is no reported computational study of the HAT between SiH 4 and p-benzyne. We optimized the reaction pathway using M06-HF/cc-pVTZ and obtained an energy barrier of 24 kcal/mol, significantly greater than the reaction by cyclobutadiene, and an energy of reaction of −3 kcal/mol. For comparison, they are, respectively, 28 kcal/mol and −1 kcal/mol from MSDFT-NOSI calculations. Interestingly, the reaction involving the CS diabatic state (Equation (16)) is formally an electron transfer from the σ Si-H bond to the lowest unoccupied BLKS orbital of the diradical reactants. Although both reactions occur dominantly via the charge transfer pathway, it is remarkable to note that through-bond interactions between the two σ-frontier orbitals of p-benzyne lead to significant energy separations between the in-phase and out-phase combinations, resulting in a CS diradical character. In the case of cyclobutadiene, on the other hand, the breaking of orbital degeneracy is due to the Jahn-Teller distortion of the molecular geometry. For comparison, the nature of orbital interactions in the OS diabatic state (Equation (17)) follows a spin-exchange mechanism between the SiH 4 and the diradical configurations ( Figure 5).

Conclusions
In this work, we described the use of a minimal active space of two electrons in two orbitals to treat diradical species using multistate density functional theory (MSDFT). Because dynamic correlation is included in the first place in each basis state function, it is possible to yield quantitative results for these systems that would otherwise require a much larger active space plus correction for dynamic correlation in wave function theory. The optimization of the Hamiltonian matrix density in terms of its total subspace ensemble energy with respect to the multistate matrix density D(r) can be performed variationally. In the present study, we adopted a nonorthogonal state interaction (NOSI) approach to approximate the self-consistent-field optimization in which both orbitals and state coefficients were simultaneously changed. We note here that since the Slater determinants in the MAS were used to represent the matrix density D(r) of interacting states of the subspace, the procedure was a state interaction rather than a configuration interaction. In this article, we first illustrated the computational procedure and its performance on the well-understood case of hydrogen molecule dissociation, and then we extended the same MAS (in terms of constrained Kohn-Sham determinant states) to other prototypical diradical cases, including singlet-triplet-energy splitting of benzyne isomers, Jahn-Teller structural tautomerization of cyclobutadiene and polyacenes up to 15 fused rings. The method was also used to investigate the competition between the charge transfer and spin exchange mechanism in the hydrogen abstraction reaction of SiH 4 by the π-type diradical cyclobutadiene and σ-type diradical p-benzyne. In comparison with accurate results, we found that MAS-MSDFT can be used to adequately model the energies and reactivities of the diradicals examined in this work, and we anticipate that this trend can be extended to other diradical species. The overall computational cost was slightly greater than that needed for N-sperate KS-DFT calculations, with N being the number of determinants in the MAS, plus the time needed to evaluate the nonorthogonal matrix element.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data that supports the findings of this study are available within the article and its Supplementary Materials. Additional information is available from the corresponding author upon reasonable request.