Spin-density calculation via the graphical unitary group approach

In this work we discuss the calculation of the spin-density matrix from fundamental spin principles as implemented in the COLUMBUS Program System employing the graphical unitary group approach (GUGA). First, a general equation for the spin-density matrix is derived in terms of the one-and two-particle reduced density matrices, quantities that are spin-independent and readily available within the GUGA formalism. Next, the evaluation of this equation using the Shavitt loop values is discussed. Finally, the spatially resolved counterpart of the spin-density matrix, the spin distribution, is calculated for the phenalenyl radical and structures produced by heteroatoms with mono- and di-substitutions. The physical meaning of the spin-density along with its computational description using various methods is discussed putting special emphasis on negative contributions to the spin-density and their quantification via a spin-promotion index. GRAPHICAL ABSTRACT


Introduction
Modern multiconfigurational electronic structure methods applied to molecules rely typically on configuration state functions (CSFs) as a basis for the expansion of the many-body wave function that describes the spatial and spin distributions of the system. The graphical unitary group approach (GUGA) [1][2][3][4][5] can be employed to construct a CSF expansion space that makes no reference to its determinantal representation. In this formalism, each CSF may be represented as a walk from the tail to the head of a Shavitt graph, or as a sequence of nodes within the graph, or as a step vector containing the occupation and spin coupling of each spatial orbital.
In the absence of external fields, the nonrelativistic molecular Hamiltonian operator commutes with the spin on M, and then to compute the property expectation value in that representation. This approach becomes inefficient for large expansions, since the determinantal basis is larger and increases faster than the CSF basis with molecular size. The best approach is to compute the spin-dependent properties directly within the GUGA formalism, making the effort to compute spin-dependent properties similar to that for spinindependent properties.
Spin-density distributions obtained from density functional theory calculations can be sensitive to the choice of functional [7,[17][18][19]. In this scenario, an ab initio approach to obtain reliable spin-densities is essential to achieve unambiguous results [8]. Furthermore, accurate ab initio spin-density distributions could be used for developing density functional theory (DFT), especially in the spin-DFT formalism proposed by von Barth and Hedin [20]. In this framework, the spin-density is employed as a fundamental property, and highly accurate results could be employed as benchmarks to construct reliable exchange-correlation functionals that approximate the exact molecular electronic structure.
The main goal of this article is to describe how to obtain the spin-density matrix within the GUGA formalism and to present results acquired from the new implementation in the Columbus Program System [5,[21][22][23]. The current implementation focuses on multiconfiguration self-consistent field (MCSCF) calculations, but it may be directly extended also to multireference configuration interaction (MRCI) and multireference average quadratic coupled cluster (MR-AQCC) calculations.
An interesting molecule to explore the recently added spin-density function in Columbus is the phenalenyl (PLY) radical, C 13 H 9 . PLY is an open-shell non-Kékule molecule composed of three fused benzene rings with a triangular topology containing 13 π-electrons. This molecule is studied widely due to its fascinating electronic structure properties [24,25], its role in soot formation [26,27], as well as potential applications in energy storage [28] and organic electronics [29,30]. Structural changes to PLY by incorporating heteroatoms into the structure drastically affects the electronic properties due to the change of the number of π -electrons and to the delocalisation of the unpaired electron [31]. Furthermore, modified phenalenyl provides the interesting possibility of inverted singlet-triplet gaps (i.e. S 1 lying below T 1 ) [32].

Spin-density matrices
In this section, the spin algebra will be reviewed and the equation to obtain the spin-density matrix from the one-particle and two-particle reduced density matrices (RDM) will be derived. For more details about this subject, the reader is referred to Ref. [33].
A wave function |ψ; S, M is an eigenfunction of the operatorŜ 2 with quantum number S ≥ 0, and of the operatorŜ z that determines the projection on the z-axis with quantum number Atomic units, in which = 1, are used herein. The quantum numbers S and M assume integer values for even numbers of electrons N and half-integer values for odd N.
It is usual to define the raising and lowering operators (Ŝ + =Ŝ x + iŜ y andŜ − =Ŝ x − iŜ y ) to further develop the projections of the spin along the axes. A sameorbital-for-different-spins (SODS) spin-orbital basis of dimension 2n is defined as {ϕ k (r) : k = 1 . . . n} ⊗ {α, β} where ϕ k (r) is an orthonormal molecular orbital (MO), and α and β are the two single-electron spin functions |1/2, ±1/2 . A member of this set is denoted ϕ kσ ≡ ϕ k (r)σ with σ ∈ {α, β}. Employing the second quantisation, creation and annihilation operators are denoted in this spin-orbital basis asâ It is convenient to collect the four SODS single-excitation operatorsâ † pαâqα ,â † pαâqβ ,â † pβâ qα , andâ † pβâ qβ into spintensor form [34,35]: with one singlet operator and the three components of a triplet operator. TheÊ pq singlet operator is a generator of the unitary group and is discussed below. With these factors and phases, the spin-tensors satisfy the standard relations The '±' convention in this equation and hereafter is that either the sequence of top signs or the sequence of bottom signs may be taken. This spin-tensor notation is useful because it allows the normal spincoupling relations to be used for operator-operator, wavefunction-wavefunction, and operator-wavefunction products. In particular, matrix elements of the form ψ ; S 3 , M 3 |T S 2 ,M 2 pq |ψ; S 1 , M 1 can be nonzero only when M 3 = M 1 + M 2 and when the three spin values S 1 , S 2 , and S 3 satisfy the triangle inequality, e.g. |S 1 − S 2 | ≤ S 3 ≤ S 1 + S 2 . In the spin-orbital basis, the spin operators can be written aŝ N σ is the number operator for σ type electrons, and it commutes with any operator that conserves the number of σ electrons. The total number operatorN =N α + N β commutes with any operator that conserves the total number of electrons.
The application of these spin operators to a wave function results in TheŜ 2 operator may be written aŝ The commutation relations, follow from Equations (6) and (7). It is useful to link the spin operators to the generators of the unitary group [2,4]. This generator may be defined with no reference to spin [36], but here it is convenient to use the formulation in Equation (5) based on spinorbitals, also follow from Equations (4), (6), and (7). The molecular Hamiltonian can be written in terms of these generators aŝ (12) in which h pq = p|h|q and v pqrs = p(1)r(2)|v 12 |q(1) s(2) are one-and two-electron integrals over the spatial orbitals. From Equation (12), it is possible to recognise the two-electron operator e pqrs =Ê pqÊrs − δ qrÊps (13) that will be used to relate the spin-density matrix to the spin-independent properties promptly obtained in the GUGA formalism. Consider the identitieŝ Performing a sum over the k index results in Recognising the operatorsN σ ,Ŝ + andŜ − [Equation (7)], Replacing theN α andN β operators by 1 2N ±Ŝ z , this equation can be written as Taking the expectation value of Equation (17), using the definitions for the charge-density matrix and for the spin-density matrix and recognising that the last two terms in Equation (17) vanish for M = S, results in the final expression for the spin-density matrix calculation for M = S in terms of 1-and 2-RDM elements, (20) No assumptions have been made regarding the wave function form in Equations (18)- (20), making this equation completely general and applicable to ground and excited states and to arbitrary full-CI and restricted expansions, including MCSCF and MRCI. Note that other density matrix normalisation conventions are sometimes used that include factors of 1/ √ 2 or 1/2.
For other members of the multiplet with M = S, the operator identities, lead to the Wigner-Eckart relation [37][38][39], Note that ψ; 0, 0|T 1,0 pq |ψ; 0, 0 = 0 because the spin values (0,1,0) violate the triangle inequality. It follows from these relations that the spin-density matrix is zero for any M = 0 spin eigenfunction, including in particular the spin-density for a singlet wave function. The members within a multiplet satisfy the relation D (1,0 The charge-density matrix in the MO basis is hermitian (symmetric for real matrices), its eigenvalues satisfy 0 ≤ μ k ≤ 2 ∀ k, and Tr(D (0,0) ) = N. It follows from Equation (11) that the charge density matrix does not depend on M, i.e. each member of a multiplet has the same charge-density matrix. Given the charge-density matrix for the wave function |ψ; S, M , the spatial charge distribution at a point r is given by The spatial charge distribution, like the charge-density matrix, is independent of the quantum number M within a wave function multiplet. Integration over all space of the charge distribution, gives exactly the total number of electrons, with no approximation due to the finite MO basis. The spin-density matrix in the MO basis is also hermitian (symmetric for real matrices), its eigenvalues satisfy −1 ≤ μ k ≤ 1 ∀ k, and Tr(D (1,0;M) ) = 2M. Given the spin-density matrix for the wave function |ψ; S, M , the spatial spin-distribution at a point r is given by Equation (25) shows that the spin-distribution for a wave function that is a spin eigenfunction with M = 0 vanishes everywhere in space. It is shown in the Appendix that this is not true for an M = 0 spin-contaminated wave function. The members within a multiplet satisfy the relation Integration over all space of the spin-distribution, gives exactly the value 2M, with no approximation due to the finite MO basis.
Taking the appropriate linear combinations of Equations (23) and (25) gives the spin-dependent spatial charge-distributions, The spin-dependent charge density matrices satisfy Thus in the GUGA formulation, the spin-dependent occupations, charge density matrices, the spin-dependent spatial charge distributions, and the spatial spin-distributions may all be computed from the spin-independent 1-and 2-RDM elements, quantities that are already available within the MCSCF computational procedure.

Natural spin-density orbitals
The spin-density matrix D (1,0;S) is hermitian (symmetric) and may be diagonalised by a unitary (orthogonal) transformation matrix U, where the eigenvalues μ (S) p are real and are assumed hereafter to be in increasing order. If the orbitals are transformed as ϕ = ϕU, then the spatial spin-distribution takes the simple form (29) in which the spatial orbital factors are nonnegative at every point in space, regardless of the orbital signs or nodal structures. These orbitals may be called the natural spin-density orbitals (NSDO), or sometimes other designations [40]. The positive μ (S) p values will therefore result in positive contributions to the spatial spin-distribution, and negative μ (S) p values will result in negative contributions. In principle, this equation holds for the inactive and virtual orbital contributions, but μ (S) p = 0 for these orbital subsets since they have equal α and β occupations, so nonzero contributions to the spin-distribution can arise only from the active orbital subspace. For CASSCF expansions, direct-product expansions, and some other common wave function expansion forms, the NSDO transformation is redundant and may be applied to resolve the active orbitals, thereby simplifying the representation of the orbitals and CSF coefficients in the wave function. The correct treatment of such redundant transformations is essential during the orbital optimisation process and also during the computation of analytic energy gradients and nonadiabatic coupling between electronic states [41]. Once the NSDOs have been determined for M = S, these same orbitals also apply to the other M = S members of the multiplet. The eigenvalues for the other members of the multiplet are given by p , which follows from Equation (22). The relation Tr(μ (M) ) = 2M holds for the eigenvalues of the NSDO basis for each multiplet member.
The interpretation of these orbitals follows from the min/max condition of the hermitian eigenvalue equation, Equation (28). Specifically, the highest eigenpair, indexed by k, results from the maximisation of the expectation value, E(u k ) ≡ u T k D (1,0;S) u k /u T k u k with respect to the elements of the vector u k . This is typically the orbital that results in the most strongly dominant α spin-density. The maximisation of E(u k−1 ) subject to the constraint that u T k−1 u k = 0 gives the orbital with the second most strongly dominant α spin-density. The lowest eigenpair results from the minimisation of the expectation value, E(u 1 ). This is typically the orbital that results in the most strongly dominant β spin-density. The second lowest eigenpair is equivalent to the minimisation of the expectation value E(u 2 ) subject to the constraint that u T 2 u 1 = 0. This is the orbital that results in the second most strongly dominant β spin-density. Continuing in this manner, each successive interior eigenpair may be regarded as the result of both a minimisation problem, subject to orthogonalisation constraints to the previously determined lower eigenpairs, and a maximisation problem, subject to the orthogonalisation constraints to the previously determined higher eigenpairs. Thus the NSDOs are those that exhibit the most extreme spin-density values. Although there are some special cases, e.g. involving wave functions that have diagonal charge and spin density matrices, for a general correlated wave function the charge-and spin-density matrices do not typically commute. This means that the natural charge-density orbitals that diagonalise D (0,0) , and thereby exhibit the most extreme charge density values, are distinct from the NSDOs that diagonalise D (1,0;S) .
A restricted Hartree-Fock (RHF) wave function is a special case of the CASSCF expansion. It is a single determinant wave function CAS(N act ,N act ) with maximal S = M = N act /2 spin and where each active orbital has occupation η k = 1. The inactive and virtual orbital subspaces are invariant as usual, but in this case the active orbital invariance results from the fact that the single CSF transforms to itself for any active orbital transformation U. The spin-density matrix in the active space is diagonal, D (1,0;S) qp = +δ pq , and therefore the spin-distribution is always given by the simple, NSDO, form of Equation (29) with positive eigenvalues. The RHF spin-distribution is therefore nonnegative at every point in space. The other members of the multiplet, generated for example by application of theŜ − operator, are generally multideterminantal wave functions, but the spin-distributions satisfy nevertheless the scaling relation, Equation (25). This means that the spin-distribution of each member of the RHF wave function multiplet is nonnegative everywhere in space for the M > 0 members, identically zero everywhere in space for the M = 0 member (for even N), and nonpositive everywhere in space for the M < 0 members. It follows that a wave function that has both positive and negative regions of spin-distribution cannot be of RHF form, or a member of an RHF multiplet, or equivalent through orbital transformation to these wave functions. Such wave functions must consist of other types of CSFs, they must include electron correlation through the mixing of CSFs, or, as discussed in the Appendix, they must be spin-contaminated.
Further insight is given by drawing an analogy to excited-state difference density matrices, widely used to characterise excited states via the analysis of attachment-detachment densities [42,43]. The analogy is apt since D (1,0;M) does indeed correspond to the difference between the D [α;M] and D [β;M] density matrices. Although the NSDO factors in Equation (29) are all positive, at a given point r in space there can be cancellations between the orbitals with positive and negative eigenvalues. These eigenpairs may be collected into two subsets: The inactive and virtual orbitals, along with any active orbitals with μ (M) p = 0 need not be included in either subset. Thus μ (+;M) ∪ μ (−;M) ⊆ {p : p = 1 . . . n}, and either or both subsets may be empty. The spindistribution in Equation (29) is partitioned analogously  (r) promotion distributions themselves provide a visual perspective of the essential contributions to the spatial spin-distribution.

Implementation in the GUGA framework
The spin-density matrix D (1,0;S) is computed from the spin-independent 1-and 2-RDM elements using Equation (20). Within GUGA, these RDM elements are constructed from individual contributions of Shavitt loops (see Figure 2 in Ref. [44], Figure 8 in Ref. [3], or Figure 4 in Ref. [45]). The charge-density matrix D (0,0) is the normal 1-RDM that is already computed and is available within the MCSCF optimisation procedure, where |ψ = m c m |m is assumed to be normalised. m |Ê pq |m is a one-electron coupling coefficient. The D The first term is computed from Shavitt loop type 14a, and the last terms are computed from loop type 11b. Consider next the off-diagonal D (1,0;S) qp spin-density matrix elements for p < q. There are five possible ranges for the index k, namely k < p < q, k = p < q, p < k < q, p < k = q, and p < q < k. These result in the contributions, The terms in Equation (37) correspond respectively to Shavitt loop types 4b, 12a, 6b, 12b, and 8b. The general computational procedure is the same as for the 1-RDM except that each contribution is accumulated both into the spin-density matrix and into the symmetrised 2-RDM, Equation (35). The p > q matrix elements are determined implicitly from the matrix index symmetry.

Computational details
The geometry of the pristine phenalenyl was first optimised using the TPSS [46] density functional and the def2-tzvp [47] basis set using the unrestricted Kohn-Sham (UKS) method. The resulting geometry presented the D 3h symmetry and the subsequent calculations were performed using the C 2v point group (subgroup of the D 3h point group). The computed wave functions have A 1 symmetry in D 3h , which correlates to A 2 symmetry in C 2v . Next, heteroatom substitutions were applied to selected positions substituting the CH group by N and NH and the structures were reoptimised. The final geometries presented C s symmetry for the N atom monodoped structure and C 2v symmetry for the N and NH didoped structures, and these point groups were employed in the following calculations. The sets of cartesian coordinates are available within the supplemental material (SM). At these optimised geometries, the spin-density distributions were obtained from single point calculations employing the unrestricted approach for the TPSS density functional [46] and the def2-tzvpp basis set [47]. This functional was employed because it was observed for iron-containing complexes that nonhybrid functionals (like TPSS) predict spin-density distributions closer to CASSCF than hybrid functionals [7,8]. This issue will be discussed in more detail below. All TPSS calculations were performed with M = +1/2. These calculations were carried out with the Orca 5.0.1 software [48].
The MCSCF calculations employing the 6-311G * * basis set [49] were performed using the same optimised geometries attained by the TPSS/def2-tzvp approach. Seven active orbitals were selected in relation to the pristine phenalenyl to make a CAS(7,7) wave function which is augmented with 3 RAS orbitals and 3 AUX orbitals. This results in the RAS(3)/CAS(7,7)/AUX(3) MCSCF wave function with S = 1/2 that describes the 13 π -orbitals and 13 π-electrons. This set of active orbitals is maintained for the doped structures, in which one or two CH groups are substituted. The resulting number of electrons obtained are as follows: (i) each N atom contributes one electron to the π -system (pyridinic doping) keeping the same number of electrons as for PLY, and (ii) each NH substitution contributes two electrons to the π -system (graphitic doping) adding one additional electron to the active space relative to PLY. Restricted Hartree-Fock (RHF) calculations are also shown. These calculations have a single active electron in a single active π-orbital, equivalent to CAS(1,1), with all other occupied orbitals inactive. These calculations were carried out with the Columbus package [23].

Results and discussion
The new implementation of the spin-density in Columbus was first applied to calculations on the PLY radical. PLY possesses 13 π-electrons meaning that it contains at least one unpaired electron. For the wave functions described below, there is typically a single MO with occupation η k ≈ 1 which is called the singly occupied molecular orbital (SOMO). As a consequence, the system is classified as a non-Kekulé benzenoid molecule, i.e. there is no valence bond (VB) structure where all π-electrons are paired. Indeed, there are seven possible resonance structures. One of these options is displayed in Figure 1, in which the unpaired electron is located in position 2, and the possible other positions for the open radical are indicated in blue. Notably, if the unpaired electron is localised at an edge-carbon, a resonance structure with a Clar sextet is present, whereas this is not the case when the radical is at the centre. Hence, one expects enhanced radical character and, therefore, spin-density at the edge positions. PLY is an alternant hydrocarbon, which means the C atoms may be partitioned into two subsets, starred and unstarred, with no direct chemical bonds between any two atoms within one of the subsets. UHF spindistributions for such systems using model Hamiltonians have been studied by Pauncz. [54] The starred C atoms correspond to the blue radical sites in Figure 1 of the VB structures. All doping substitutions were chosen to keep an odd number of electrons since singlet systems present a zero spin-density matrix [Equation (22)]. Figure 1 displays the carbon numbering scheme for the PLY structure. Doping was always performed in positions 2 and 12 since these (and the symmetry equivalent positions) are the atoms where the singly occupied molecular orbital is located in the resonance structures. The spin-distributions obtained with the MCSCF, UHF, and RHF methods as well as the PBE50, PBE0, PBE and TPSS functionals are displayed in Figure 2 for the pristine PLY. The wave functions were computed with A 2 symmetry within the C 2v point group, but they display the full A 1 symmetry within the D 3h point group. Colour is used to distinguish between the regions of positive and negative spin-distributions. Blue is used to designate the α-rich positive regions, and red is used to designate the β-rich negative regions. All spin-distributions are computed with M = +1/2, thus the positive regions overall dominate the negative regions. The formal radical resides within the α-rich region, and additional spinpolarisation is reflected through the further positive and negative spin-distribution regions.
First considering the spin distribution for the pristine PLY (Figure 2), all methods present alternance of positive and negative spin densities on neighbouring carbons residing mostly in the π -system, a similar pattern as also seen for triangulene structures with larger number of rings [55]. The exception to this is the RHF result in Figure 2(c) which shows only positive spindistribution values. As discussed in Section 2.2 this is the expected result for an RHF wave function with M > 0. For the other methods, there are 7 α-rich and 6 β-rich sites, consistent with the starred and unstarred C atoms in the alternant hydrocarbon. The VB framework provides some insight into the spin distributions. Here, spin population on a position can be seen as the result of a weighted combination of the resonance structures presenting open radicals on that site. Enhanced positive spin (blue) is found for the seven positions that can hold the free radical in the VB picture as indicated in blue in Figure 1. As discussed above, the radical on the central atom does not allow the formation of a Clar sextet. Thus, this resonance structure should possess a higher energy and present a lower contribution to the spin distribution. This difference is clearly represented in the MCSCF  Figure 2(b)] displays a strikingly different spin-distribution, which is, first, significantly more pronounced than the MCSCF one and, second, shows a diminished difference between outer and inner atoms (0.784e vs. 0.724e). RHF [ Figure 2(c)] presents the opposite limiting case with no spin-density in the centre. Proceeding to UKS, we find that reducing the amount of HFX in the functional along the series PBE50 (d), PBE0 (e), PBE (f) reduces the overall spin-density on the atoms and restores the difference between the outer and central carbon atoms. Finally, TPSS (g), which is a local meta-GGA based on PBE, provides a similar appearance to PBE. In agreement with Refs. [7,8], we find that using little or no HFX tends to improve the results. Specifically, the Mulliken population analysis for the MCSCF and TPSS calculations differs by at most 0.0267e (Table S1), highlighting that TPSS is a suitable choice for the description of this molecule. Figure 2 illustrates the challenges present in describing spin-densities using UHF/UKS approaches. Dramatic differences and qualitatively diverging results are observed between different functionals, and, were the MCSCF reference not available, there would be no simple way to assess the quality of these descriptions. From a more formal viewpoint it is worth noting that our GUGA-based MCSCF strictly produces spineigenfunctions. The negative spin-distributions seen in Figure 2(a) arise via electron correlation. By contrast, negative spin-distributions from a UHF/UKS single determinantal treatment can only arise via spincontamination, i.e. the α and β spin-orbitals have different shapes and, as a consequence, the resulting determinant is no longer an eigenfunction ofŜ 2 , (see the Appendix for further details). Following Equation (1), the expectation value Ŝ 2 for a doublet state should be 0.75, and the deviations from that value can be used to quantify the spin-contamination. Increasing the HFX, we find that spin-contamination increases sharply and, specifically, Ŝ 2 values of 0.766 (PBE), 0.824 (PBE0), 0.973 (PBE50), and 2.01 (UHF) are obtained. Further discussions of spin-contamination in UKS theory are quite subtle [12,56,57] and outside the scope of this work, but two points are particularly relevant. First, the KS determinant is an auxiliary construction whose purpose is to produce values for charge and spin densities whereas its Ŝ 2 expectation value does not have a physical meaning.
There is no reason to expect that Ŝ 2 would reflect the true spin of the system. Second, UHF, as a wave function theory, is required to produce a valid Ŝ 2 expectation value, [Equations (A6) and (A15)]. In this sense, the spincontamination in UHF is more troublesome than in UKS. Therefore, UHF as such, and by extension its admixture to UKS within a hybrid functional, is more problematic ) for the pristine phenalenyl computed for the 2 A 2 ( 2 A 1 ) ground state obtained from the MCSCF calculation. The isovalue is ±0.001e · Å −3 for ρ [−;1/2] ± (r) and ±0.05 Å −3/2 for the NSDOs. Blue and red represent positive and negative values, respectively. than using the original 'pure' UKS theory. Hence, in line with the observed results, a pure functional seems more suitable than a hybrid functional for the description of spin-densities.
The RHF results in Figure 2(c) are consistent with an uncorrelated (single determinant in this case) spineigenfunction with spin-promotion index (1/2) = 0e in Equation (33). The RHF spin-distribution arises only from the delocalisation of the SOMO, with no contributions either from electron correlation or from spincontamination. The RHF SOMO has a 2 (a 1 ) symmetry in the C 2v (D 3h ) point group. The central carbon atom lies at the intersection of three vertical nodal planes of the SOMO and therefore acquires no significant spindistribution from this delocalisation. In contrast, this central carbon has significant positive spin-density due to correlation in the MCSCF wave function in Figure 2(a).
The NSDOs of the MCSCF wave function, the spinpromotion distributions, and the spin-promotion numbers for PLY are shown in Figure 3. The natural chargedensity orbitals obtained from the MCSCF calculation are displayed in Table S2. The spin-promotion numbers p The positive spin-promotion distribution is concentrated on the six alternating edge carbons, as suggested from the VB analysis, and also consistent with the delocalisation of the a 1 SOMO. A smaller positive spin-promotion distribution occurs on the central carbon, a consequence of the valence correlation. This density does not arise from the SOMO, which has three intersecting vertical nodal planes centred on that atom. The smaller negative spinpromotion distribution is spread among the other three edge carbons and the remaining three carbons, consistent with their unstarred C-atom designations. Due to these spatial separations, there is relatively little cancellation of positive and negative spin-promotion distribution to form the total spin-distribution ρ [−;1/2] (r) shown in Figure 2(a). The square of the orbital amplitude times the μ (1/2) p eigenvalue for each orbital gives the contribution of that orbital to the spin-promotion distribution. Thus all of the orbitals on the left in Figure 3 contribute only positive values and all of the orbitals on the right contribute only negative values to the respective spinpromotion distribution and to the total spin-distribution.
The spin distributions for the N-doped structures are displayed in Figures 4 and 5. These substitutions are isoeletronic to PLY and lead to a pyridinic form, not affecting the number of electrons in the π -system of the molecule. These substitutions break the D 3h symmetry of the molecule, reducing it to C 2v for the C 11 N 2 H 7 di-substitution and to C s for the C 12 NH 8 mono-substitution.
Concerning the spin distribution for the doping by two N atoms (Figure 4), there is no appreciable qualitative change in comparison to the pristine case. The MCSCF spin-promotion index is (1/2) = 0.552e, which is comparable to the PLY value. For both the MCSCF and TPSS methods, based on the Mulliken analysis, the spin population decreases in position 2 of the molecule in comparison to PLY, from 0.25e to 0.21e for the MCSCF  method and more drastically from 0.28e to 0.18e for the TPSS functional, as seen in Table S3. The natural chargedensity orbitals and NSDOs calculated with the MCSCF method are displayed in Tables S4 and S5.
For the spin-density of the single N-doped structure ( Figure 5), a comparison of the spin distribution obtained from the RAS(3)/CAS(7,7)/AUX(3) calculation with the pristine counterpart ( Figure 2) shows that the MCSCF method predicts a slightly larger spin distribution on the nitrogen and central carbon atoms than the TPSS functional (Table S6). The former has values equal to 0.22e (N) and 0.07e (C) while the latter has values equal to 0.20e (N) and 0.02e (C). The MCSCF spin-promotion index is (1/2) = 0.560e, which is comparable to the PLY value. The natural charge-density orbitals and NSDOs attained from the RAS(3)/CAS(7,7)/AUX(3) wavefunction are displayed in Tables S7 and S8. Figure 6 shows the NH di-substitution structures. This molecule presents a C 2v symmetry and adds two more electrons on the π-system, resulting in nine electrons in the complete active space, RAS(3)/CAS(9,7)/AUX(3). The SOMO for this wave function has a 2 symmetry, and the wave function itself has A 2 symmetry. The MCSCF spin-promotion index is (1/2) = 0.160e, which is much less than the PLY value and the previous N-substitution values. Within the 13-orbital active space, the lowest five NSDOs have negative μ (1/2) p eigenvalues and the remaining eight NSDOs have positive eigenvalues. A spin-distribution closer to an ROHF wave function is observed in Figure 6. According to results obtained from highly correlated multireference methods, this graphitic doping is expected to have the potential of enhancing the biradicaloid character for larger polyaromatic hydrocarbons [58][59][60][61][62]. Moreover, among the doping possibilities considered here, this type of substitution is the only one that changes qualitatively the α and β alternance in neighbouring atoms ( Figure 6). For both methods considered, 10 of the 12 edge atoms present α spin while only the carbons in position 1, 7 and 13 (central) present β spin. These carbon atoms lie in the vertical nodal plane of the a 2 SOMO, and thereby receive no significant spin-distribution from the delocalisation of this orbital. For the TPSS functional, a small β spindistribution is also observed between the carbons (i.e. within the chemical bonds) presenting σ character. The RAS(3)/CAS(9,7)/AUX(3) wave function cannot validate these contributions since there are no σ orbitals in the active space. This feature will be examined in more detail in future work using more accurate wave functions. The spin population for this structure is listed in Table S9, and the MCSCF and the DFT values differ by no more than 0.024e. The natural charge-density orbitals and the NSDOs obtained from the RAS(3)/CAS(9,7)/AUX(3) method are displayed in Tables S10 and S11.

Conclusion
In this work, the theoretical background for the spindensity matrix calculation is presented. The matrix equation is written in terms of one-and two-particle charge density matrix elements that do not depend on the spin component quantum number and are available within the GUGA formalism. This approach avoids the projection of the wave function in a basis of spindependent Slater determinants, and it uses information that is already computed and is available within the MCSCF optimisation procedure. This approach is also applicable to MRCI and MR-AQCC methods.
The properties of the spin-density matrix and its spatial spin-distribution are discussed. Among several insights, the spin-distribution for a wave function that is a spin-eigenfunction with M = 0 vanishes everywhere in space, while this is not true for spin-contaminated M = 0 wave functions. In the latter case, it is only the integral over all space that is zero. Moreover, M = 0 multiconfigurational wave functions that are spin-eigenfunctions are able to produce spin-distributions with both positive and negative regions due to electron correlation, in contrast to unrestricted wave functions in which this feature arises due to spin-contamination. The natural spin-density orbitals provide a convenient basis for the discussion of the spin-distributions. Spin-promotion distributions and spin-promotion numbers provide insight into the spin-polarisation of the MCSCF wave function.
The implementation of the spin-density matrix calculation for the MCSCF method in Columbus was applied to the phenalenyl benzenoid radical and three doped structures, substituting the CH group by the N-atom and NH group. The pristine phenalenyl structure presents alternating αand β-rich regions on neighbouring carbons, as do the N-atom mono-and di-substituted radicals. These radicals all display comparable spinpolarisation. In contrast, the NH di-substituted radical does not show alternating regions, and it displays much smaller spin-polarisation.
When compared to DFT functionals, it is observed that UKS functionals with more Hartree-Fock exchange produce exaggerated spin-distributions, and by comparison with the MCSCF spin-distribution, a pure functional like TPSS is more reliable than hybrid functionals.
We believe this work provides solid basis on how to obtain spin properties from the multiconfigurational GUGA formalism and provides trustworthy tools to explore and shed light on the spin properties of molecular systems. and they are not necessarily related to each other throughŜ ± or through simple spin-tensor relations such as Equation (6).
Another important difference between spin-eigenfunctions and spin-contaminated wave functions is that typically D [−;χ ;0] = 0, for M = 0, in contrast to Equation (22). In principle, a spin-contaminated wave function may be expanded in a spineigenfunction basis, with S min = |M| and S max = min 1 2 N, n − 1 2 N . This expansion is discussed in more detail in Ref. [64]. In general, the spin-eigenfunction basis functions |ψ; S, M are multideterminantal expansions in the DODS spin-orbital basis, even when the wave function |ψ; * , M is a single UHF determinant. In terms of this expansion, the UHFŜ 2 expectation value Equation ( The S = |M| term in this expansion is typically the desired wave function, and any other nonzero terms, each of which is positive so there are no cancellations within the summation, are regarded as spin-contamination. For a DODS wave function expansion, the spin-density operator in the ϕ [α] orbital basis may be written This operator, along with its counterpart operator in the ϕ [β] basis, is consistent with the AO spin-density matrix D x S x S+1 ψ; S, 0|D 1,0 qp +D 1,0 pq |ψ; S + 1, 0 . (A20) These spin-density matrix elements are composed entirely of spin-contamination contributions from off-diagonal adjacent terms in the spin-eigenfunction expansion. The nonzero eigenvalues of this matrix must be both positive and negative in order to satisfy Tr(D (1,0;M) ) = 2M = 0, and thus a nonzero spatial spin-distribution at some point r in space can be either positive or negative, with integration over all space giving zero. The possibility of a nonzero spatial spin-distribution is also apparent from Equation (A12). This is in contrast to the M = 0 spindistribution for a wave function that is a spin-eigenfunction, which vanishes at every point in space, Equation (25).