General Correlated Geminal Ansatz for Electronic Structure Calculations: Exploiting Pfaffians in Place of Determinants

We propose here a single Pfaffian correlated variational ansatz that dramatically improves the accuracy with respect to the single determinant one, while remaining at a similar computational cost. A much larger correlation energy is indeed determined by the most general two electron pairing function, including both singlet and triplet channels, combined with a many-body Jastrow factor, including all possible spin–spin, spin–density, and density–density terms. The main technical ingredient to exploit this accuracy is the use of the Pfaffian for antisymmetrizing a highly correlated pairing function, thus recovering the Fermi statistics for electrons with an affordable computational cost. Moreover, the application of the diffusion Monte Carlo, within the fixed node approximation, allows us to obtain very accurate binding energies for the first preliminary calculations reported in this study: C2, N2, and O2 and the benzene molecule. This is promising and remarkable, considering that they represent extremely difficult molecules even for computationally demanding multideterminant approaches, and opens therefore the way for realistic and accurate electronic simulations with an algorithm scaling at most as the fourth power of the number of electrons.


INTRODUCTION
The accurate determination of the many-electron wave function (WF) has always been a challenging task starting from the early stage of quantum mechanics. 1 So far, several attempts have been made toward this direction ranging from CCSD(T) 2 to tensor network 3 and density matrix renormalization group (DMRG), 4 up to the very recent breakthrough with the use of machine learning methodologies. 5 All these schemes pay the price of being computationally demanding, with a computational complexity ranging from a large degree polynomial of the number of electrons to exponential complexity.
Quantum Monte Carlo (QMC) techniques for electronic structure calculations have proven to be very successful in describing the electronic correlation encoded in a many-body WF. 6−8 In particular, the variational Monte Carlo (VMC) 9−11 samples the real space electronic configurations of the considered system with a probability distribution given by the WF square, thus providing the efficient evaluation not only of the total energy but also of the expectation values of most commonly used many-body operators. Within VMC, it is possible to improve the description of the ground state (GS) WF by minimization of the total energy expectation value. The WF obtained can be used as it is or further improved by the diffusion Monte Carlo calculation (DMC) method. 10−13 This technique is a projection algorithm performed statistically using the information on the sign contained in the given WF, dubbed here as a guiding function. In this way, we can considerably improve the description of the GS, projecting on the lowest possible energy WF with the same signs of the guiding WF. In the ideal case of a guiding function that, for every configuration, has the same sign of the GS, the abovedescribed DMC algorithm provides the exact solution. 10,11 In the framework of QMC, different ansatzs are used to approximate the true GS WF, with the purpose to achieve an affordable compromise between the accuracy of the calculation and its computational cost. Though a good representation of the GS can sometimes be achieved with a simple and "cheap" WF, in most cases, the use of a very complicated and computationally demanding ansatz is necessary to get a correct answer.
Slater determinants (SDs) are the simplest fermionic WF used for QMC. They provide a single particle picture of the quantum many-body problem, preserving the Pauli principle, that is, the Fermi statistics for electrons. They can be obtained directly from mean-field calculations. Unfortunately in many situations of interest, it is not possible to give a good description of the system in terms of a single SD. 14 −17 In QMC, there are two possible strategies to overcome this problem: the use of a linear combination of different SDs 17−21 or different ansatzs with larger variational freedom. 22,23 The multideterminant WFs can be systematically improved and in principle can describe exactly every GS with a large enough number of SDs. Unfortunately the number of SDs that has to be taken into account scales exponentially with the number of electrons preventing the calculations on large systems. 24 The use of pairing function replaces the single particle description of the SD with a richer one in terms of electron pairs. The corresponding WF is a natural extension of the single SD ansatz and represents a direct and efficient implementation of the Anderson resonating valence bond (RVB) 25 theory of many-electron WFs. In particular, it provides a direct description of the singlet and triplet correlations that are absent in the SD. Depending on the definition of the pairing function, qualitatively different WFs can be obtained. They will be described in Section 2, where we will focus also on the technical details required for the calculation. We will introduce the symmetric antisymmetrized geminal power (AGP) 22,26,27 and the broken symmetry antisymmetrized geminal power (AGPu), 28 but we will mainly focus on the most general AGP. In the previous literature, 23 ,29 it has been indicated as Pfaffian WF, and people have been referring to the AGPs as AGP, but since the AGP (or Pfaffian WF) literally realizes the most general antisymmetrized geminal power, we dub this case with the shortest acronym, that is, AGP. It will be shown that this approach becomes very efficient in combination with an explicit correlation term, known as a Jastrow factor (JF), 14,30,31 that promotes or penalizes the bonds according to the electronic correlation. As it will be shown later, we have introduced a quite general JF, depending both on spin and electron charges. When it is applied to an AGP without definite spin, it allows its almost complete restoration, mimicking in this way a spin projection operation that, though approximate, is much cheaper than other approaches. 32,33 Even if the pairing functions cannot be improved systematically, these WFs have a much larger variational freedom than the SDs, nevertheless remaining with a similar computational cost.
If on one hand, for the multideterminant WF, the calculation can be computationally very expensive, on the other hand, for the pairing functions, the optimization of a large number of nonlinear variational parameters can become a serious limitation if not handled efficiently. Indeed, in order to exploit the full potential of these ansatzes, it has been fundamental to use the most recent techniques for the calculation of the derivatives and optimization strategies.
In a previous attempt, the AGP WF was used by exploiting only a very small fraction of the large variational freedom of the ansatz. 23,29 The results were not encouraging, and the energies obtained with this ansatz did not improve the ones of the AGPs that, in turn, has a lower computational cost. Despite the Pfaffian was no longer used in the electronic system to our knowledge, the experience with lattice models has shown that the AGP WF is able to improve considerably the description of magnetic and correlated systems. 34 Moreover, the introduction of a powerful JF and the recent results obtained in combination with the AGPs 14,35,36 encouraged us to look for the unexpressed potential of the full AGP WF.
In this paper, we will compare the results obtained with AGP WF with available state of the art VMC and DMC calculations. In particular, we benchmark our WF on the diatomic molecules with corresponding high spin atoms in the first row of the periodic table, i.e., carbon, nitrogen, and oxygen, and on the benzene. The first ones are systems that, despite their apparent simplicity, represent useful benchmarks for many highly correlated methods. 37−40 We will show that, with the use of our best WFs, even with a very compact basis set, we are able to achieve an accuracy comparable with the state of the art multideterminant WFs at a computational cost similar to the one of a single SD. Not only the total energies and the dissociation energies are extremely accurate but we also analyzed the magnetic proprieties of these molecules, unveiling the unexpected rich physics behind these systems. Finally we consider the benzene molecule, a system that represents the prototypical example of the RVB theory and thus a fundamental test case for our approach.

WFS AND PROCEDURES
For all calculations we present in this paper, we used the TurboRVB package for QMC calculations. 41,42 The WFs used for this work are factorized as the product of a fermionic mean field and an explicit bosonic correlation factor. Being Ψ(X) the WF of a given configuration X = (r 1 σ 1 , r 2 σ 2 , ..., r N σ N ) of N electrons of spins σ i and positions r i , we can write Ψ(X) as where Φ(X) takes into account the fermionic nature of the electrons, while J(X) is the JF: an exponential modulation of the WF that substantially improves the electronic correlation description for all types of WFs studied here. The fermionic term of the WF, dubbed as Φ(X) in eq 1, is the most important part, directly encoding the behavior of the electrons while imposing the antisymmetrization under particle exchange. In the following section, we will describe the basis set used, the definition of the AGPs, AGPu, and AGP after a brief introduction to the SD. Finally we will discuss the JF correlator.
2.1. Basis Set. We expand our ansatz in an atom-centered basis set of Gaussian orbitals for the calculation of the JF and a hybrid basis set for the fermionic part of the WF, as it will be discussed below. The Gaussian orbitals basis set is indicated as {ϕ I,ν (r)}, with each element being the νth orbital centered on the Ith atom at the position R I . The elements in the basis set have the form where Z ν is a numerical coefficient that describes how diffused the atomic orbital is around the atom, while Y l ν ,m ν is the spherical harmonic function with angular quantum numbers l ν and m ν , corresponding to the orbital type ν which is always assumed to be real. This basis set has been used without further contractions for the description of the JF. Instead, for the fermionic part of our WF, we have used hybrid atomic orbitals (HOs) 26,27 to expand them over a richer set of Gaussian orbitals and, by means of the contraction, considering only an affordable number of variational parameters. The HOs, indeed, are obtained as linear combinations of all the elements of the Gaussian basis set used for a given atom, labeled by I r r ( ) ( ) The above hybrid orbitals allow us to take into account the modification of the standard Slater orbitals corresponding to isolated atoms, to the case when they are instead placed in a complex environment. Therefore, we have chosen to use a number of hybrid orbitals equal to the single particle ones occupied in the absence of electron−electron interactions and including also all the ones corresponding to the same shell of degenerate one particle levels. The corresponding orbitals are the ones that should physically play a role in the considered electronic systems. Hence, in all of the first row molecules, we have considered the full hybridization of five atomic orbitals, coming from two s-wave and three p-wave ones, that can be corrected by several components with much higher angular momenta. This is because the full spherical symmetry is no longer satisfied even in a simple homonuclear molecule. For the sake of compactness, we indicate here all basis elements as {ϕ k (r)} combining the indices ω and I, and I and ν in a single index k for a lighter notation. Every time we refer to the AGPs, AGPu, and AGP, it is meant to be a basis of HOs.
The exponents Z ν have been chosen from the ccpVDZ or ccpVTZ basis set according to this criterium: the contraction are removed and all exponents with Z ν > 150 a.u −2 . are eliminated. This is possible because contracted orbitals containing very large exponents are necessary only with a pure Gaussian basis in order to satisfy the electron-ion cusp conditions. They are instead appropriately considered by the one-body term of our JF, as described in Section 2.4. The exponents chosen are then further optimized at molecular equilibrium distance and kept fixed in the corresponding atomic calculation (where the optimization of the exponents has an almost negligible effect) and the dispersion energy curves.
2.2. Slater Determinant. Here, we will provide a preparation description of the SD that is important both for the initialization of the pairing function and for comparing our results with the existing literature. 17 From theoretical and computational point of view, the simplest fermionic WF is the SD, called Jastrow SD (JSD) in the presence of a JF. The SD is built from the vacuum state by populating a number of orthogonal single particle molecular orbitals (MOs) equal to the number of electrons in the system. Henceforth, we omit the spin indices, by assuming that each spin component corresponds to a different SD. In our basis, the MO are in the form The MOs can be obtained directly from a density functional theory (DFT) or Hartree−Fock calculation, but they can also be further optimized with VMC. 43 It is well known that the antisymmetric product of these MO leads to the determinant of the matrix in which every MO is evaluated for each electron position i k j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z For weakly correlated systems, the JSD can often give reasonably good results with an affordable computational cost and a limited number of variational parameters. It is also a common choice to use a linear combination of SDs to improve the description of the WF, with ansatzs that take different names depending on the type and number of SDs considered. In this paper, we will compare directly the results of our WFs to the ones obtained with one of the most successful multideterminant WFs, the full valence complete active space (FVCAS) WF.
2.3. Pairing Function. The use of the pairing function in correlated WFs allows an electronic description that goes beyond the single particle picture of the SD. The building block of this WF has the following general form where all the elements of the matrix λ represent most of the WF variational parameters. They depend on the orbitals considered and on the spin σ 1 , σ 2 of the so called geminal function f. In principle, when we break the spin symmetry, the basis sets used for ↑ and ↓ electrons can be different, otherwise the chosen basis does not depend on the spin component. In order to set up a consistent many-body WF starting from the geminal, several choices are possible depending on the criteria adopted for the definition of the geminal. To highlight the different possibilities, we can recast eq 6 in a way in which the spin dependency is more explicit where f f f f f r r r r r r r r r r In order to satisfy the Pauli principle, we have f ± (r 1 ,r 2 ) = ±f ± (r 2 ,r 1 ) and f σ (r 1 ,r 2 ) = −f σ (r 2 ,r 1 ) for σ = ↑, ↓. Our WF is then obtained by antisymmetrizing the product over all electron pairs considered that, by definition, occupy the same pairing function. For simplicity, we will enumerate the spin up electrons from 1 to N ↑ and the spin down ones from N ↑ + 1 to N. As suggested by the name AGP, our goal is to define a WF that is literally the antisymmetrized product of the geminals and the unpaired orbitals (if present), namely where α is one of the possible way of distributing the N electrons between the p/2 pairs and the N − p unpaired orbitals Θ, and Sgn(α) is the sign of the corresponding permutation of the particles that is required to ensure the fermionic behavior. In particular, different choices of the pairing function, obtained by excluding one or more terms in the eq 8, lead to different ways to compute eq 9. These choices also impact quantitatively and qualitatively on the kind of physics that we can describe by means of this type of WF. Therefore, we will distinguish among three distinct cases: if we consider only the singlet term in eq 8 we obtain the AGPs, if we include the singlet and the S z = 0 triplet term we have the AGPu, while the most general case is just the definition adopted here for the AGP. 2.3.1. AGPs. Let us consider for the moment the unpolarized case N ↑ = N ↓ , the extension to the polarized cases will be straightforward and will be discussed later on. When no triplet correlations are allowed, we build our WFs using only singlet pairs, and the pairing function in eq 7 contains only the symmetric element f + f f r r r r ( , ) 1 2 ( ) (, ) In this case, we project a perfect singlet that we denote as AGPs. The λ matrix elements in eq 10 are nonzero only for σ 1 ≠ σ 2 , and they are symmetric for spin exchange. In order to calculate the AGPs, we can write all the possible combinations of pairs of opposite spin electrons in a matrix defined as i k j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z In this way, each row of the matrix corresponds to an electron of spin ↑ and each column an electron of spin ↓. The definition of the matrix F in this form is convenient because it allows the antisymmetrization requested by the eq 9 in a simple and efficient way. Indeed, it can be demonstrated 26 that the correct antisymmetrization of the pairs considered in this case is given by This is somehow intuitive because we want to sum all the possible products of N/2 matrix elements of F, where in all these factors, a column element or a row element is present only once, exhausting all the possible configurations of the system considered with an appropriate ± sign that, in this case, is just given by the one corresponding to the determinant of F. When the system is polarized and N ↑ ≠ N ↓ , we cannot build the solution using only the singlet terms because the matrix F written as in eq 11 is a rectangular matrix and its determinant cannot be computed. Supposing for simplicity that N ↑ > N ↓ , in this case, we have to add a number N ↑ − N ↓ of unpaired spinup MOs {Θ i (r)} not only for fulfilling the polarization required but also, most importantly, to turn the matrix F to a perfectly defined square matrix i k j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z Also in this case, a consistent antisymmetric WF can be again calculated as the determinant 26 of the matrix F exactly in the same way as the singlet pairing in eq 12.

AGPu.
For the AGPu, only the parallel spin term of the triplet component are omitted. This means that the spin symmetry is broken, and a magnetic order parameter can be directed along the z-quantization axis. This WF is called the broken symmetry AGP (AGPu), and the difference from the previous AGPs is the presence of the antisymmetric f − component in the definition of the pairing function in eq 7, that for this case is In order to define this pairing function, we break the spin symmetry in the opposite electron spin case with σ 1 ≠ σ 2 , by keeping equal to zero the σ 1 = σ 2 components of eq 6. With exactly the same procedure used in the case of the AGPs, depending on the polarization, we can build the same matrix F of eqs 11 or 13 that is no longer symmetric now. Even in this case, the correct antisymmetrized sum of these pairs is given by the determinant. 26 Thus, analogously to eq 12, we obtain that implements the simplest broken symmetry ansatz based on the pairing function. 2.3.3. AGP. The AGP (also known in the literature as Pfaffian WF 23 ) is in our opinion the most important pairing function, being the most general one and encoding new variational freedoms into the AGPs and the AGPu. We will show that it represents the most powerful description of the chemical bond within the paradigm developed in this work. This WF represents also the most general mean-field state, that is, the GS of a mean-field Hamiltonian containing BCS anomalous terms projected on a given number N of particles and total spin projection S z i N i tot 1 σ = ∑ = along the zquantization axis. In this case, the definition of the pairing function is exactly the one in eq 7, containing all terms including the parallel spin terms of the triplet. This means that now, when we build the AGP, we have to include also the parallel spin electron pairs in the WF. In this way, the AGP can also describe a magnetic order parameter in any direction of the space, and thus, it is also possible to rotate the spin component of the WF in any direction. This will allow us to break the symmetry along the spin-quantization axis and then rotate it. As we will explain later, this plays a crucial role when we use this WF in combination with our JF because it allows us to preserve the total S z of the molecules and include spin fluctuations.
Of course, we cannot create a WF using only pairs if the number of electrons in the system is odd, so, for the moment, let us assume N is even. The extension to the odd number of electrons is trivial and will be discussed immediately after. We will dub as W the matrix containing all possible pairs i k j j j j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z z z z where the matrix is antisymmetric for the fermionic commutation rules, and thus, the elements of the diagonal are set to zero. We can recast the W highlighting its different spin sectors as where W ↑↑ and W ↓↓ are respectively N ↑ × N ↑ and N ↓ × N ↓ antisymmetric matrices that take into account the parallel spin terms of the triplet, while W ↑↓ is a N ↑ × N ↓ matrix such that W ↑↓ = −W ↓↑ T , describing the remaining triplet and singlet contributions. In the case of AGPs and AGPu, we can also build a similar matrix where the matrices W ↑↑ and W ↓↓ are identically zero.
Analogously to the case of the AGPs and AGPu, we have to identify a way to calculate the antisymmetric product of all pairs considered. In this case, it is easy to identify the antisymmetrization procedure defined in eq 9 as the Pfaffian of the matrix W. After introducing this algebraic operation, the reason will be straightforward to the reader.
The Pfaffian is an algebraic operation acting on antisymmetric square matrices with an even number of rows and columns. Being N even, the matrix W satisfies these hypothesis. The usual definition of the Pfaffian requires the introduction of the concept of partition of the matrix W (18) where all i k and j k are different, i k < j k for each k and i 1 < i 2 < ... < i N . The sign(α) is given by the permutation that orders the vector of the indices {i 1 , j 1 , i 2 , j 2 , ..., i M , j M }. In this way, all indices are considered only once. The Pfaffian is then defined as where the sum over α is extended over all possible partitions. However, an alternative definition 44 of the Pfaffian can better clarify the correspondence to the eq 9. It can indeed be defined alternatively as where P now represents a generic permutation of the possible row and column indices of the matrix without any constraints, and the sign(P) is the parity of the permutation. In this definition, it is easy to recognize the antisymmetrized sum corresponding to the eq 9. Let us introduce now a further property of the Pfaffian that will be useful in the following section. We will indicate with 0 a m × m matrix containing only 0 and B a generic m × m matrix, we have For odd number of electrons, it is necessary to use a spindependent unpaired orbital Θ σ (r) so that we can accommodate the remaining electron that is not considered by the product of the pairs. The unpaired orbital introduces a supplementary row and column to the matrix W. Being Θ ↑ = (Θ ↑ (r 1 ), Θ ↑ (r 2 ), ..., Θ ↑ (r N↑ )) the vector containing the values of the unpaired orbital Θ ↑ at the ↑ electron positions and Θ ↓ = (Θ ↓ (r N ↑ +1 ), Θ ↓ (r N↑+2 ), ..., Θ ↓ (r N )) the one calculated for the ↓ electron ones, we modify the matrix in eq 17 as i k j j j j j j j j j j j j j j y Also in this case, the permutation sum implied by the Pfaffian leads to the correct antisymmetrization required from eq 9. The matrix W satisfies the hypothesis of the calculation having an even leading matrix dimension N̅ = N + 1. We can further notice that no assumption has been made on the polarization of the system, and so no unpaired orbital is required except for a single one in the case of odd N.
It is however possible in principle to introduce further pairs of unpaired orbitals, if, for example, we want to describe AGPs or AGPu with a full AGP WF. We define Θ iσ (r) as the set of the considered m unpaired orbitals and ) the vector containing the values of the unpaired orbital Θ i,↑ for the ↑ electron positions and Θ i ↓ = (Θ i,↓ (r N ↑ +1 ), Θ i,↓ (r N ↑ +2 ), ..., Θ i,↓ (r N )) the one calculated for the ↓ electron ones. We can modify the matrix in eq 17 as i k j j j j j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z z z z z that is a N̅ × N̅ matrix where N̅ = N + m. We can again antisymmetrize this product using the definition of the Pfaffian provided in eq 20. A careful reader could have noticed that, by applying the Pfaffian definition, we are antisymmetrizing not only over the electron indices but also over the orbital indices of the unpaired orbitals. This antisymmetrization, however, contains the one over the physical electrons, and leads therefore to a physically allowed electronic WF. Moreover, we can notice that, by using the previous definition, we can identify the AGPs and the AGPu as subcases of the general AGP. Indeed, by using the expressions of the pairing function and the unpaired orbitals of the AGPs and AGPu, we obtain W ↑↑ = 0, W ↓↓ = 0, Θ i ↓ = 0, and N̅ = 2N ↑ . By merging eqs 13 and 17, we can define 24) and this means that, by applying eq 21, we immediately obtain Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article W F Pf( ) det( ) = ± (25) where the sign only depends on the number of electrons and is constant, thus irrelevant. This shows in a straightforward way that the AGPs and AGPu defined in the previous subsection are nothing but particular cases of the most general AGP.
2.4. Jastrow Factor. Within QMC, it is easy to improve the quality of the WF by multiplying the WF with an exponential JF. This last one enriches the description of the GS by encoding explicitly the electronic correlation, while speeding up the convergence to the complete basis set limit. 7 Indeed, with an appropriate choice, the JF can satisfy exactly the electron−electron and electron-ion cusp conditions of the many-body WF, consequences of the Coulomb 1/r singularity at a short distance. In this paper, we introduce a new kind of JF that contains a richer dependence on the spin and that plays a fundamental role when used in combination with the AGP WF. The JF is defined as where U ei is a single body term that deals explicitly with the electron−ion interaction and U ee is a many-body term that properly accounts for the electronic correlation. The single body term is In eq 28, Z I is the atomic number of the atom I and b ei is a variational parameter, while g I (r i ) encodes the most general nonhomogeneous electron−ion one-body term, that is, depending explicitly on all nuclear and electron coordinates and not only on their relative distances, that is defined as where the summation is extended over all Gaussian orbitals in the JF basis set centered on the Ith atom. The electron− electron term instead is written as where the sum is extended over the pairs of different electrons and where with the 2 × 2 matrix b ee  In our expression, the first term in eq 31 deals explicitly with the electron−electron cusp conditions, the second term in eq 31 instead is a bosonic pairing function in the form (32) with the elements of the matrix ζ defining further variational parameters. Notice that both g I and g ee do not affect the cusp conditions because they are expanded over cuspless Gaussian orbitals. The g ee term has the same form of eq 6, but since the fermionic behavior is already encoded in the fermionic part of the WF, this term is symmetric under particle exchange. The use of a pairing function in the JF enriches the description of the charge and spin correlations of the system, noticeably improving the quality of the global WF. It is a common practice to adopt a simplified or even absent spin dependency in the function u of eq 31. This is often accurate for systems where the magnetic properties are not relevant. We will refer to it below with the prefix Js in the WF, in contrast with the prefix J used for the full spin-dependent JF. A perfect singlet remains as such after the multiplication of a spin-independent Jastrow, and so our spin-dependent JF is not appropriate if we do not want to break the spin symmetry. It is, instead, necessary if we want to recover, at least approximately, the singlet from a spin contaminated broken symmetry ansatz. A general spin-dependent u, as defined in eq 31, is therefore of fundamental importance for the AGPu or the AGP ansatzs.
Let us start with a simple example. We consider two atoms with opposite spins and break the spin symmetry by orienting the spins of the atoms along the z-quantization axis. In this case, the JF is not able to change the classical antiferromagnetic spin state because it acts as an irrelevant constant when applied to it. It is instead more physical to orient the spin moment of the atoms in a direction perpendicular to the quantization axis chosen for the JF. In this way, the JF can act on the electrons and the spins while the magnetic moment is free to fluctuate and recover its genuine quantum character. As previously mentioned, with the AGP it is possible to rotate the spin of the WF in every direction and orient the magnetic moment in any direction of the space. This works particularly well in combination with our JF that can suppress the unfavored triplet configurations with parallel spins generated by the rotation. This optimal spin orientation of the atoms, that is perpendicular to the JF one, is rigorously valid within the wellknown spin-wave theory of a quantum antiferromagnet. 34 In this case, the JF defined with a spin-quantization axis perpendicular to the magnetic moment of the atoms allows the description of the quantum fluctuations and the corresponding zero point energy (ZPE), even for a finite (as is our case) number of atoms. 34 2.5. Procedure. The first step to calculate and optimize our WF is to identify a reasonable starting point. We chose to start from a DFT calculation because of its flexibility. We have used LDA calculations for spin symmetric systems, while for the ones with opposite spin antiferromagnetic moments we have broken the symmetry with a LSDA calculation, by adopting an appropriate initialization of the WF. The SD obtained from a DFT calculation is mapped without loss of information into AGPs or AGPu and then in a second analogous step, we convert the AGPs and AGPu into a full AGP.
For the first conversion, let us consider eq 6. If we compute it in the basis set of the MOs obtained from the DFT, we have namely only the diagonal terms in the matrix λ̅ are present. Moreover, we can also remove the spin dependency if there is no symmetry breaking. For the polarized case, the unpaired orbitals are the last occupied MOs. By substituting the definition of the MOs with their expansion over a localized atomic basis set, as given in eq 4, we can recast the above equation exactly in the same form shown in eq 6 with a matrix λ When we convert AGPs or AGPu into an AGP WF, we already have an initialization for the sectors of the pairs with different spins W ↑↓ and W ↓↑ that can be obtained directly from the AGPs or AGPu pairing functions. The main challenge is to find a reasonable initialization for the two sectors W ↓↓ and W ↑↑ that are not described by the AGPs or AGPu.
There are two different procedures that we can follow: the first one is used for polarized systems, and the second one instead is preferred in the case of broken spin symmetry and in the presence of antiferromagnetism, namely, molecules well described by opposite atomic magnetic moments. If there is no antiferromagnetism and the polarization is such that |S tot z | < 1, the W ↓↓ and W ↑↑ are instead identically zero. This holds true not only for S tot z = 0 but also for S tot z = ±1/2, where the single unpaired MO used in eq 22 acquires also a spin dependency, not present in the AGPs and AGPu cases.
Obviously the atoms, and also the O 2 molecule, do not have antiferromagnetism, but, on the other side, they have a net polarization. We can build the W ↑↑ block of the matrix using the two unpaired orbitals Θ 1 and Θ 2 for the definition of the parallel spin matrices of the AGPs or AGPu in the following way where the presence of the minus sign guarantees the pairing function to be antisymmetric under particle exchange, while the λ̅ is an arbitrary scaling factor that has no influence on the final value of the WF. Once we map the unpaired orbitals in the desired basis set, we obtain the variational parameters of the matrix λ for the ↑↑ sector.
In the presence of opposite atomic magnetic moments, it is possible to rotate the spin component of the pairing function to initialize the W ↓↓ and W ↑↑ sectors. As we mentioned earlier, a further effect of this operation is to direct the atomic magnetic moments in a direction perpendicular to the spinquantization axis. It is worth mentioning that, within our method, the spin orientation with respect to the molecular axis is irrelevant because in a nonrelativistic Hamiltonian, the spin− orbit coupling is not present. In this case, we have chosen to work with the atomic magnetic moments perpendicular to the zaxis, hence we applied a rotation of π/2 around the yd irection. This operation maps ( ) a n d ( ) If we apply this transformation to the pairing function from eq 14, we obtain This transformation provides a meaningful initialization to our AGP WF that now has to be optimized to reach the best possible description of the GS. Indeed, within VMC, it is only thanks to the optimization that we can improve the description of the GS. So far, we have only converted the DFT WF from one ansatz to another, but the key for the success of this procedure is the optimization of all possible variational parameters. It is indeed crucial to optimize not only the ones corresponding to the matrix λ and the JF parameters but also the coefficients of the hybrid orbitals μ and the exponents of the Gaussian basis set Z ν . This is realized computationally in a very efficient way using a coding technique called adjoint algorithmic differentiation 45 that allows calculations of total energy derivatives with respect to all variational parameters involved in a given algorithm that computes only the energy. This is remarkably done by paying a very small slowing down of a factor ≈2−3 with respect to the latter algorithm. We have also used a state of the art optimization scheme 46,47 for a correct search of the energy minimum. Remarkably, even when there is some possible dependency among the many variational parameters considered in our ansatz, the stochastic reconfiguration technique remains stable and efficient, thanks to an appropriate regularization of the stochastic matrix S. 10 Once we calculate the variational minimum, the best description of the GS is then obtained with the DMC calculation.
Even considering that the number of variational parameters involved in the calculation may be quite large, the optimization has a very small impact over the total computational cost, that is indeed mostly given by the DMC for all cases reported in this work. In Table 1, we compare the computational cost of the DMC calculations for the different WFs considered. We notice that the JAGP is even less expensive than the JSD and JsAGPs WFs. The JAGP and JAGPu are so efficient because, in this case, the variance of the energy is considerably smaller, as we can see from Table 1. This implies that JAGP and JAGPu require a smaller number of DMC iterations to reach the desired accuracy because they have lower variance compared with the JsAGPs and the JSD, thanks to the spin-dependent JF.
For large number N of electrons, the DMC calculation should scale as N 4 for fixed total energy accuracy, and the main question, that we have not studied here, is whether the optimization remains computationally negligible because the number of variational parameter scales as N 2 . In this respect, we have experienced that an optimization technique performed with a slow but very stable method, that is, with a large number of "cheap" optimization steps, each one determined by a relatively small number of samples (even much smaller than the number of parameters) is very promising for future large N applications.
Finally we introduced a technique to deal with particularly unstable AGP WFs. Indeed, it is possible, after a very large number of optimization steps (>10,000), that some eigenvalues of the matrix λ become too small as compared with the largest reference eigenvalue. This creates some instabilities in the inversion of the matrix W required for the QMC fast updates. For this reason, by an appropriate use of the PFAPACK library, 49 we identified a procedure to map the diagonalization of a full skew-symmetric matrix λ to the one corresponding to a real tridiagonal symmetric matrix. After this mapping, we can use the most powerful and stable LAPACK routines for diagonalization. Indeed, most linear algebra packages cannot deal with antisymmetric matrices, and a general diagonalization tool was not available for this case. The introduction of this procedure, described in detail in the Appendix, allows us to describe the matrix λ in terms of eigenvalues and orthogonal orbitals playing the role of eigenvectors of an antisymmetric matrix. We will refer to them in the Appendix as MOs because it may be considered their formal definition, within the formulation introduced in this work. With this meaningful decomposition, we can finally regularize the matrix λ by replacing the too small eigenvalues with reasonable lower bounds and continue, if necessary, with the optimization of the variational parameters.
2.6. S 2 Operator. The basic concept of QMC relies on the real space configurations sampling of a general electronic system. All observables can be indeed calculated in the basis where the electron positions and their spins are defined. In particular, for the systems considered, it is interesting to estimate the spin observables in order to understand their magnetic properties and the quality of the corresponding WFs. If during the simulation the value of S z is fixed, when we break the symmetry the value of the S 2 is instead the result of the interplay between the JF and the AGP or the AGPu. The efficient computation of the expectation value of the S 2 operator has already been described in ref 50 for the JsAGPu and will be shown now for the JAGP.
Here, we show how to evaluate S 2 in a region of the space with a fast and computationally cheap approach based on the fast update algebra of the AGP and the spin-dependent JF. Let us consider the expectation value of the S 2 operator over a generic WF Ψ by direct application of its definition. Here, we use the completeness of the spatial configurations where the summation symbol implies a 3N-multidimensional integral over the electron coordinates. Assuming a fixed polarization, we can write the explicit expression of the total spin square as where the operators S i in the above equation act on the spin component corresponding to the electron position r i of the configuration X. We can notice that (40) and that, by using QMC sampling, we generate configurations according to the probability density p(X). Thus, we can evaluate the above multidimensional integral by directly sampling the estimator S 2 (X) that multiplies p(X) in eq 39, as follows The content of the former equation can be evaluated efficiently as we will explain in the following section. Indeed, the application of the operator S i + S j − to the configuration X generates only a configuration X ij = {(r 1 ↑), ..., (r i ↓), ..., (r j ↑), ..., (r N ↓)}. Considering X our sampled configuration and using the previously given definition of X ij , we can recast eq 41 as The only hard challenge of eq 42 is the calculation of the for i = 1, 2, ..., N ↑ and j = N ↑ + 1, N ↑ + 2, ..., N, that in our case read The configurations X and X ij differ for a spin flip of the electrons i and j, but we can also consider X ij as the configuration in which the electron i evolved to the position previously occupied by j and vice versa. We can then calculate the ratios in eq 44 using a fast algebra to update two positions for the AGP and for the JF with a direct evaluation based on Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article the Sherman−Morrison algebra and some simple manipulations, as discussed in detail later on. It is also possible to calculate the value S 2 (Λ) of the S 2 operator in a sub-region of the space Λ. For this quantity, we can obtain the similar expression to eq 42 (45) where N σ Λ (σ = ↑, ↓) is the number of σ-electrons in the region Λ, N Λ = N ↑ Λ + N ↓ Λ . The summation symbol over i ∈ {Λ, σ} indicates the sum for all σ-electron whose coordinate is in the region Λ. Therefore, we can use the same method as described below.
2.6.1. AGP Contribution. To calculate the AGP contribution to r ij , we were able to find a slim and fast algebra making an extensive use of the Pfaffian properties. 51 It was fundamental to find an efficient algebra to calculate the whole matrix of the ratios r with a computational cost that is O(N 3 ), by using mostly BLAS3 operations, thus avoiding that this computation could become the bottleneck of the whole procedure. In this way, we could ensure the evaluation cost of S 2 to be comparable with the one of a typical QMC cycle over all N electrons that is at most O(N 3 ). Before describing the fast updating rules for the position of two electrons with a single move, we need to introduce some quantities fundamental for the calculation.
Let us denote W −1 as the inverse of W. This inverse W −1 can be computed from scratch for each configuration used to sample the spin square. The electron coordinates r i are given for i = 1, ..., N, but since the corresponding spin can change with respect to the original choice (↑ for i ≤ N ↑ , and ↓ for i > N ↑ ) because of the spin flips mentioned in the previous subsection, we will consider explicitly the values of the spin here.
We then define the matrix θ as For the spin ↑ electrons, we can define the vectors i k j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z while for the spin ↓, we have instead i k j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z We can use these vectors to build the N × N matrix that allows us to define and finally Now, we have all ingredients that we need for our fast updating algebra, and upon application of Sherman−Morrison algebra, we arrive at the ratio We can notice that the preliminary calculation of the auxiliary matrices θ, V, U, and D, including the inversion of W, amounts to a total of O(N 3 ) operations, while the calculation of the ratios is O(N 2 ) once the matrices have been computed. 2.6.2. JF Contribution. In the JF that we introduced in the previous section, only the two-body term of eq 30 has a spin dependency, and thus, only this part contributes to the ratio. By simple substitution, it is easy to prove that where we have defined The whole operation has a O(N 2 ) computational cost and so does not limit the calculation in terms of performances.

RESULTS AND DISCUSSION
We apply this new approach for two types of systems: the first row high spin atoms (carbon, nitrogen, and oxygen) and their diatomic molecules and the benzene molecule. The first ones still represent useful benchmarks for the quantum chemistry approach, and a reasonable description of their properties and binding energies requires very expensive multireference methods. It is therefore very interesting to test our approach to find whether we are able to obtain a good description with a single Pfaffian ansatz. Benzene molecule on the other side is the most famous and important example of the RVB theory, so it represents a fundamental benchmark test for a method inspired by this theory. Here, we compare our results with exact available solutions, JSD WFs and with the JFVCAS multideterminant expansions for QMC. We will also show that our WF satisfies the size consistency both at VMC and DMC levels, a primary requirement if we want to use this approach for more challenging chemical studies. 3.1. Carbon. Carbon dimer is probably the most interesting example discussed in this study. A full understanding of the behavior of the carbon−carbon interaction is still missing, and the bond order of this molecule is still under debate. 52 The role of the spin fluctuations in this molecule has already been discussed, 53 but we believe that it is very instructive and represents the most important achievement of the JAGP WF. Indeed, it is only thanks to the spin fluctuations that we can have a correct description of its dimer bond.
The carbon atoms have spin triplet electronic configurations, and their mutual interaction leads to a singlet molecule. As we can see from Figure 1 and Table 3, the JAGP not only Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article improves the results of the JSD WF but remarkably also the description given by the JsAGPs and JAGPu. The huge difference between the multideterminant expansion JFVCAS and the JSD binding energies helps to quantify the effect of the multideterminantal nature of this molecule, and what makes this even more surprising is that the quality of the results obtained with a single JAGP WF, with a computational cost comparable to a SD, is already very close to the exact value. As already mentioned before, the explanation for the impressive improvement of the binding energy from JsAGPs and JAGPu to JAGP resides on the description of the strong spin fluctuations in this molecule. The JAGP gives a very accurate picture of its magnetic properties as we can see from Table 2, giving results very close to S 2 = 2 for the atom and S 2 = 0 for the molecule. Conversely, by using the JsAGP (the AGP without spin-dependent JF) and the JAGPu, we cannot recover the singlet from the broken symmetry initialization. Interestingly, as expected, the molecule does not have any magnetic moment on the z-direction because it is an almost perfect singlet. The atomic spins, localized around each atom, point in opposite directions in order to form the singlet molecular state. Because there is no magnetic moment along z, we can measure its magnetic moment only by separately evaluating the S 2 in the two semi-infinite regions, each one containing a single atom, separated by a plane perpendicular to the molecular axis and at the same distance from the two atoms. In Figure 2, we show that, even at bond distance, there is a very strong magnetic moment around the atoms and, in this way, we can explain the strong effect of the ZPE of the spin fluctuations described by the JAGP.
Moreover, Figure 2 shows that only with the JAGP WF, we have a size consistent solution with the molecule that recovers the energy of two independent atoms at large distance. This feature is fundamental if we want to use this WF to describe chemical reactions and perform large-scale simulations, with a size consistent behavior at large distances. The importance of the variational optimization of the WF is particularly evident in this small molecule. With the standard approach, by applying DMC to a SD taken by DFT (here obtained with Purdue and Zunger LDA 56 ), a level crossing in the occupation of the π MOs occurs at around 3 bohr distance, above which the π bonding orbitals are only partially occupied. This implies clear artifacts in the DMC energies. We have verified that this level crossing is reproduced with a standard DFT-LDA calculation by Gaussian 16 A.03 revision 57 and an almost converged basis set (the standard cc-pVQZ). The level crossing has also been observed in ref 58. In our variational optimization instead, we have verified that it is important to start at large distance with    DMC energy dispersion of the carbon dimer: only the JAGP allows the system to be size consistent at large distance, which means that it is able to recover the energy and the expectation value of the S 2 operator of two isolated atoms at bond distance; however, the carbon atoms maintain a large value of S 2 . The sharp change of the projected S 2 value at around 3 a.u. is probably due to an avoided crossing of two energy levels belonging to the same irreducible representation, in agreement with DMRG. 59 Within LSDA, this effect is reproduced by a discontinuous change in the occupation of the π orbitals in the corresponding SD. Lines are guides to the eye. the WF predicted by LSDA, otherwise a sizably higher energy is obtained. This effect is reflected also by the sharp change of the projected S 2 at around 3 bohr distance (see Figure 2), that could be compatible with an avoided crossing between two energy levels belonging to the same 1 Σ g + representation. 59

Nitrogen.
Nitrogen is in some sense similar to the carbon case: also its dimer is indeed a singlet formed by two large spin 3/2 atoms.
As we can notice from Figure 1 and Table 4, at the DMC level, the JsAGPu and JAGP are both exact within chemical accuracy. All our calculations compare with the exact result better than the JFVCAS solution. Surprisingly, at the VMC level, the binding energies calculated with JAGP, JsAGPs, and JAGPu are also very good.
We remark that a very powerful method, as the recently proposed Fermi net 5 (a neural network-based WF), cannot reach the same precision in the binding energy even if the total energies of the molecule and atom are the best available ones. This clearly shows that all our ansatzs allow a remarkable cancellation of errors, when computing the total energy differences between the molecule and the two independent atoms.
In this case, however, the difference between JAGP and JsAGPs/JAGPu is much smaller than in the previous case and should be related to a less important role of the spin fluctuations and also to a smaller magnetic moment of the atoms at equilibrium distance. By repeating the reasoning done for the carbon dimer, we can quantify the magnetic moment from the S 2 value in the semi-infinite region separated by a plane perpendicular to the axis of the molecule and equidistant from the atoms. As shown in Figure 3, at bond distance, the S 2 of the atom is much smaller than the one of an independent atom, and therefore, even if the nitrogen atom has a large spin, when it is forming a dimer it does not give rise to a strong antiferromagnetism.
Also in this case, it is important to notice that the JAGP solution is size consistent both in energy and spin. Despite the very good description at bond distance provided by the JsAGPs, we notice from Figure 3 that it is not perfectly size consistent. Within our approach, a fully consistent picture and a very accurate dispersion are possible only by means of the JAGP ansatz, that is able to work properly also in the strong correlation regime at large interatomic distance.
3.3. Oxygen. The oxygen is very different from the previous cases but nevertheless very interesting for different reasons. The oxygen dimer consists of two triplet atoms, but this time the molecule is a triplet. There are small atomic magnetic moments in the GS of the oxygen molecule, but the role of the magnetic interaction remains important, as shown by the application of the JAGP ansatz. In this case, it looks that the interaction of parallel spin electrons is particularly important, and this can be described by the JAGP ansatz more accurately than the corresponding JsAGPs and JAGPu ones, as discussed in the previous sections. Thus, we expect to recover with the JAGP some correlation that we miss when we simplify the ansatz by using the unpaired orbitals in the JsAGPs and in the JAGPu WFs.
As shown in Figure 1 and Table 5, at the DMC level, the energies obtained with the JAGP WF are extremely good even for the oxygen dimer. In this case, the correct description of     ansatz, appears to be fundamental. Indeed, the final result is so accurate that the binding energy is comparable to the one obtained with the multideterminant JFVCAS WF. It is even more surprising that the absolute energies of the atom and molecule are very close where not even better than the ones provided by the multideterminant expansion both at VMC and DMC levels. We have to point out, however, that within JFVCAS method, it is not possible to improve the JSD atom 17 and that the binding energy slightly better than the JAGP one derives from the poorer quality of the atom rather than a better description of the molecule. The problem of the size consistency for the oxygen dimer is absolutely nontrivial and even more complicated than the previous cases. Starting from bond distance, we have a molecule of spin one and fixed projection S z = 1, but we have to recover the behavior of two independent atoms. This means that, by keeping the projection S z = 1 constant, while separating the atoms far apart, we have to recover the correct atoms of spin one and thus we need to have one atom with the spin oriented in a direction perpendicular to the z-axis. This is impossible for the JsAGPs and the JAGPu but allowed by the JAGP, a remarkable and absolutely nontrivial feature of this WF. As we can see from Figure 4, at large distance, only with the JAGP the system recovers the energy and the spins of the independent atoms, showing that, by means of our advanced optimization tools, it is possible to dramatically change the WF up to the point of rotating completely the spin of an atom.
3.4. Benzene. The benzene molecule represents one of the most successful examples of the RVB theory with the carbon− carbon bonds resonating among several valence bond configurations, for example, Kekuléand Dewar. QMC methods are able to provide a very good description of this important molecule, 61,62 and thus, it is interesting to check whether, with our new approach, we can obtain a very accurate result. In particular, in Table 6, we compare the results obtained by JSD, JAGP, JAGPu, and JsAGPs WF, showing that all results obtained with a pairing function (from JsAGPs to JAGP) provide a very good estimate of the absolute energies, noticeably improving the results of the JSD. Moreover, the corresponding atomization energies are extremely accurate at the DMC level, whereas the JSD largely overestimates it. It is finally interesting to notice that, even if there is a sizeable gain in terms of absolute energy with our best ansatz, that is, the JAGP, it is not clear why this systematic improvement does not sizeably affect the atomization energy, likewise this could be almost converged to the exact value. This might be in principle explained because, at present, the accuracy of the state of the art "estimated exact" calculation is probably not enough to establish an energy difference ≪0.1 eV. For instance, the ZPE has been estimated by DFT 63 and some work is certainly necessary to clarify this issue, for example, by calculating the ZPE directly with QMC.
We remark here that the JsAGPs description of the benzene molecule is already very accurate and it is not improved by the JAGP. This is probably due to the lack of any sizeable spin moment around any atom composing this molecule. Indeed, the S 2 value calculated for the JAGP and JAGPu solutions are 0.032(1) and 0.0123 (7), respectively, proving that any local magnetic moment is almost completely melted during the optimization, despite its nonzero initialization. We conclude therefore that in the benzene molecule, the spin fluctuations are not relevant and the use of the Pfaffian leads only to a marginal improvement of the total energy while the molecule is correctly described by a perfect singlet RVB ansatz given by the JsAGPs, in agreement with the classical RVB picture by L. Pauling. 64

CONCLUSIONS
In this work we have proposed a new WF for QMC calculations given by the most general fermionic pairing function ansatzs in combination with a spin JF that provides a very rich description of the electronic correlation by means of a bosonic pairing function complementary to the fermionic one. With a computational cost comparable to a SD, we were able to improve not only the results achieved with a simple JSD but also with JsAGPs and JAGPu, reaching a level of accuracy comparable to the one obtained with the multideterminant JFVCAS WF. The powerful optimization techniques are probably the keys to explain the remarkable improvement we obtained with this WF, compared to previous attempts. 23,29 In particular, we have shown that the JAGP ansatz provides a very accurate description of high spin atoms and their dimers and that it is size consistent. This should increase the number of possible applications, providing a reasonably accurate and computationally feasible tool for studying chemical reactions. In the plot, also the expectation value of the projected S 2 operator on the atoms for the JAGP that recovers the value of two isolated atoms at large distance. Lines are guides to the eye. The triplet correlations have proven to be necessary to take into account correctly the ZPE of the spin fluctuations that we can now correctly describe thanks to a physical and accurate setup obtained by orienting the atomic magnetic moments of the AGP in the direction perpendicular to the spinquantization axis chosen for the JF. For this reason, we have obtained a very good description of the carbon and nitrogen dimers, remarkably even when the first molecule was found to be very poorly described by the JsAGPs and the JAGPu. Moreover, it is only thanks to the presence of the triplet correlations that we were able to improve the description of the oxygen dimer as a strongly correlated triplet molecule with a highly entangled spin interaction among the atoms. Comparison with other methods different from QMC is shown in App. A. Our QMC variational energy is much better than state of the art quantum chemistry methods that seem to be affected by strong basis set errors even when considering only energy differences. For instance, the total energy difference ΔE at R = 4.2 a.u. and R = 2.11 a.u. in Table 8 should be close to the estimated exact binding energy (i.e. ≃9.91 eV from ref 55) at most weakly corrected by the residual dispersive interaction. Both DMRG and MRCI clearly miss more than 1 eV with the DZ basis, that is, ΔE ≃ 8.49 eV. In order to show more clearly that the discrepancy between our DMC results and DMRG and MRCI is actually an artifact of the small basis, we have carried out UCCSD(T) calculation both for small (ΔE = 8.6 eV) and large (ΔE = 9.55 eV) basis set, and, as expected, our calculation (ΔE = 9.63) is much more in agreement with the most accurate large basis set calculation. In any event, our binding energy for N 2 (9.933 ± 0.006 eV) is surprisingly more accurate than the best state of the art calculation with CCSD(T) (9.73 eV from the Computational Chemistry Comparison and Benchmark Data-Base 65 ), implying that, most likely, our results should be considered to be the state art for the full dispersion curve of these small molecules.
Finally we demonstrated that for the benzene dimer, the JAGP is able to provide a very accurate atomization energy, though it is not clear in this case whether the triplet correlations are crucial for a highly accurate calculation. However, it is important to highlight that the accuracy in the binding energy is always much better than the accuracy in the total energy and that therefore there exists always a remarkable cancellation of errors in the total energy differences. This feature indeed is fundamental for a compact ansatz like the JAGP, and it challenges other very expensive highly correlated methods, even when these are able to achieve almost exact total energies as it was the case for the Fermi net approach to the nitrogen dimer.
The relatively low computational cost of QMC combined with powerful optimization techniques, allowing a reasonably large number of variational parameters, makes this approach ideal for studying systems even much larger than the ones considered in this work. Indeed, we believe that the paradigm presented in this paper could represent in the future a very powerful tool to investigate the electronic structure of interesting chemical compounds and physical systems where the spin interaction may play an important role, that in turn may be a number much larger than previously believed, as we have presented here the C 2 molecule as the very first and remarkable example of an antiferromagnetic chemical bond.

■ APPENDIX Energy Dispersion Comparison
Comparison between the energy dispersion calculated with DMC JAGP WF and unrestricted single reference coupled cluster (UCCSD-T) with ccpVDZ and ccpV5Z basis sets. We further compared the carbon energy dispersion with DMRG,  Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article heat-bath configuration interaction (HCI) and full configuration interaction (FCI) from the literature, 59,66,67 and the nitrogen dispersion with multireference coupled cluster (MRCC) and DMRG. 68 UCCSD(T) calculations were performed using Gaussian 16 A.03 revision with the counterpoise correction, with the frozen-core approximation and the full-core correlation. 57 Table 7 and Figure 5 show that there are significant discrepancies between different methods in the carbon dimer dispersion curve at large distances. However, one has to consider that even in a quadruple zeta basis, the binding energy D e = 6.22 eV 67 is about 6 mH lower than the estimated exact one and therefore if we reference all curves at the bond length minimum energy, as reported in the mentioned figure, a method that is supposed to be weakly dependent on the basis, as our DMC, should be slightly higher in energy at large distance, provided it remains close to the exact dispersion energy curve. Moreover, there may be sizeable corrections due to the frozen-core approximation employed by DMRG, HCI, and FCI. We have indeed verified that they are non-negligible in the UCCSD-T calculation, implying that core−valence interaction can lead to a further nonparallelity error of about 3 mH (see Figure 5). Core−valence interaction is considered in our DMC calculations simply because, within this technique, it is not possible to employ the frozen-core approximation. Nevertheless, it is clear that our results may have some error, but it is remarkable that if we use the corresponding energy values for computing the ZPE of the dimer, we find excellent agreement with the experimental value, given by 0.1146 eV. 69 Indeed, the ZPE calculated values, using a standard fit with a quartic polynomial close to the equilibrium distance, are 0.1153(6), 0.108, 0.106, 0.112, 0.114, and 1133(3) eV for our DMC, UCCSD-T full-core, UCCSD-T frozen-core, DMRG, HCI, and FCI, respectively. In summary, by taking into account all possible sources of error, we believe that our results are in reasonable agreement with the expected "exact result" converged in the complete basis set limit and with full core− valence interaction taken into account. Indeed, we believe that only a more direct comparison with experiments or a full-core FCI/DMRG or HCI extrapolated to the complete basis set limit can further improve the accuracy of the dispersion curve.

Diagonalization of a Skew Symmetric Generally Complex Matrix λ
Here, we discuss a general procedure to transform a generic complex antisymmetric matrix into a canonical Youla form that represents the equivalent of the standard diagonalization of Hermitian matrices. This is obtained by means of an appropriate unitary matrix U defined by an orthonormal set of states that we will address in the following MOs. Given a N̅ × N̅ antisymmetric matrix λ, our goal is to identify a set of p paired states {(ϕ j 1 , ϕ j 2 )} of orthonormal MOs, such that p ≤ N̅ is even and 56) where the LHS of the above equations indicate standard matrix vector products, with shorthand notations adopted also in the remaining part of this Appendix. In this basis, we can write any skew-symmetric matrix λ in the canonical Youla form i k j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z using only p/2 strictly positive parameters a j . These ones play the same role of the eigenvalues for an ordinary Hermitian matrix and henceforth we will use this name for them, even if the matrix λ MO is not diagonal but represents the simplest nonvanishing skew-symmetric matrix. The transformation of the original matrix λ to the corresponding canonical Youla form by means of an appropriate unitary transformation λ = U*λ MO U † provides us also a very simple way to regularize the matrix λ, as discussed in the main text. In the case of odd N̅ , it will be shown later that there exists always an eigenvector of λ with a vanishing eigenvalue, but the decomposition remains possible as λ MO will contain at least one vanishing row and corresponding column. Here, we define that an eigenvector is singular if it corresponds to a vanishing eigenvalue as in the odd N̅ case.
It would be ideal for this calculation to use a very robust and stable diagonalization routine to maintain machine accuracy for the MOs. Unfortunately these routines are not commonly available for antisymmetric matrices and thus several mathematical transformations are necessary to map our task to a sequence of more commonly used or at least easily available algorithms.
A generic N̅ × N̅ antisymmetric matrix is written in the following way i k j j j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z z z The first step is to transform λ in a tridiagonal antisymmetric real matrix. This operation is implemented in the subroutine zsktrd (dsktrd) contained in the PFAPACK library. 49 The use of the Householder algorithm allows us to decompose the generic matrix λ as U U 1 Tr 1 λ λ = * † (59) where U 1 is the unitary transformation matrix output of the algorithm, while λ Tr is a tridiagonal real antisymmetric matrix written in the standard tridiagonal form i k j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z Thus, we can multiply the matrix λ Tr for the imaginary unit i, yielding a more conventional tridiagonal Hermitian matrix λ iH , defined by purely imaginary matrix elements.
We highlight that it is possible to map the matrix λ iH into a real Hermitian matrix via a unitary transformation and use the appropriate LAPACK routine for its fast diagonalization. This procedure is well known and will be discussed later.
At this point, we can use the spectral theorem for Hermitian matrices to decompose the matrix λ iH = ψλ diag ψ † , where λ diag is a diagonal matrix containing in its diagonal part the real eigenvalues a i of λ iH , and ψ is the unitary matrix, where each column is given by the eigenvector in the principle complex, corresponding to each eigenvalue, in the chosen order. This decomposition implies: However, because the matrix ψ is generally complex and ψ † ≠ ψ T , some manipulation is necessary if we want to satisfy the skew-symmetry property of λ in an easy and transparent way.
If we consider one eigenvector ψ ̅ j associated to an eigenvalue a j > 0, we have where we have used that both λ Tr and the eigenvalues a j are real. This means that if ψ ̅ j is an eigenvector of λ iH relative to the eigenvalue a j , then ψ ̅ j * is an eigenvector corresponding to the eigenvalue −a j and thus orthogonal to ψ ̅ j because of the orthogonality between eigenvectors of an Hermitian matrix corresponding to different eigenvalues ± a j . We can thus easily verify, by using the relations given in eqs 62 and 63, the following simple equations Once we have identified all pairs corresponding to all positive eigenvalues a j > 0, we can write the unitary matrix ψ ̅ , that is now real, by adding the remaining eigenvectors (that can be also chosen real, as shown in subsection Singular Eigenvectors) with vanishing eigenvalues in the remaining rightmost columns. In this way, we can finally define a unitary real matrix ψ ̅ yielding λ Tr = −iλ iH = ψ ̅ λ MO ψ ̅ T , where λ MO is defined in eq 57 and therefore by using eq 59 which represents the desired decomposition because the product of two unitary matrices U* = U 1 *ψ ̅ remains a unitary matrix and its transponce U † coincides with ψ ̅ T U 1 † , yielding λ = U*λ MO U † , that is the purpose of this Appendix.

Triangular Hermitian Matrices: A Mapping from Imaginary to Real
In order to use the LAPACK routines for diagonalization, we have to map the tridiagonal fully imaginary Hermitian matrix λ iH , defined only (the diagonal elements are zero to fulfill Hermitianity) by its upper diagonal elements ib j with b j real for j = 1, 2, ..., N̅ − 1, into a tridiagonal real symmetric matrix λ R . We can implement this mapping by applying a unitary transformation to the matrix λ iH . For this purpose, we introduce the following transformation described by the matrix The matrix U 2 is a complex diagonal matrix defined as i k j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z The explicit calculation of the right-hand side of the eq 69 gives i k j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z that therefore implies that λ R , with the above definition, is a real symmetric matrix. At this point, we can diagonalize the matrix λ R by means of a real unitary matrix U R , that is the output of a standard LAPACK diagonalization routine of tridiagonal real matrices (e.g., dstevx for double precision arithmetic). In this way, λ iH can be diagonalized as λ iH = , where λ diag is a diagonal matrix containing the corresponding eigenvalues of the LAPACK diagonalization.
Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article

Singular Eigenvectors
Within this formulation, it is also particularly easy to compute all the real singular eigenvectors of λ iH corresponding to the possible vanishing eigenvalues. They were used in this Appendix to complete the columns of the unitary real matrix ψ ̅ . From the outcome of the previous subsection, any eigenvector ϕ k j of λ iH can be obtained by applying the diagonal matrix U 2 to a real eigenvector ϕ̅ k j of λ R , namely, ϕ k j = ϕ̅ k j exp(−iπ/2(k − 1)), implying that even k-components are purely imaginary and odd k-components are purely real. Then, it is obvious to realize that ϕ̅ k j corresponds to a singular eigenvector of λ R , and also ( ) ] correspond to singular eigenvectors or at most null vectors (not both) of λ iH because this matrix is purely imaginary, and the complex conjugation of a singular eigenvector is again a singular eigenvector by eqs 62 and 63 with a j = 0.
Then, it follows that all orthogonal eigenvectors ϕ̅ j (output of dstevx) corresponding to the zero eigenvalues of the matrix λ R can be used to define the real singular eigenvectors corresponding to the matrix λ iH , that is that are explicitly real and orthogonal to each other because ∑ k ϕk j ϕk l = ∑ k ϕ̅ k j ϕ̅ k l = δ l,j . They are also orthogonal to all the other pairs of non-singular eigenvectors because of the orthogonality property of eigenvectors of an Hermitian matrix λ iH that we have already used in the previous section.