Extraction of a One-Particle Reduced Density Matrix from a Quantum Monte Carlo Electronic Density: A New Tool for Studying Nondynamic Correlation

: In this work, we present a method to build a ﬁrst order reduced density matrix (1-RDM) of a molecule from variational Quantum Monte Carlo (VMC) computations by means of a given correlated mapping wave function. Such a wave function is modeled on a Generalized Valence Bond plus Complete Active Space Conﬁguration Interaction form and ﬁts at best the density resulting from the Slater-Jastrow wave function of VMC. The accuracy of the method proposed has been proved by comparing the resulting kinetic energy with the corresponding VMC value. This 1-RDM is used to analyze the amount of correlation eventually captured in Kohn-Sham calculations performed in an unrestricted approach (UKS-DFT) and with different energy functionals. We performed test calculations on a selected set of molecules that show a signiﬁcant multireference character. In this analysis, we compared both local and global indicators of nondynamic and dynamic correlation. Moreover, following the natural orbital decomposition of the 1-RDM, we also compared the effective temperatures of the corresponding Fermi-like distributions. Although there is a general agreement between UKS-DFT and VMC, we found the best match with the functional LC-BLYP.


Introduction
The first measure of electron correlation energy was introduced by Löwdin as the energy difference between the exact solution of the non-relativistic Schrödinger equation and the Hartree-Fock (HF) approximation [1]. Although this is not a rigorous measure, considering that some correlation, namely Fermi correlation, is included in HF theory, it accounts for almost all the so-called Coulomb correlation originated by the two-body electron-electron interaction in the non-relativistic Hamiltonian. The correlation energy can further be decomposed in two terms denoted dynamic and nondynamic (or static) correlation. As the HF description based on a single Slater determinant, the dynamic correlation enters when the ground state of the system is sufficiently well described by a single determinant. On the other hand, nondynamic correlation arises when, in order to reproduce accurately the electronic structure of the system, two or more near degenerate determinants are necessary. Usually, the so-called post Hartree Fock methods (such as perturbation expansions or coupled clusters methods) are able to recover most of the dynamic correlation of the system; however, the recover of nondynamic correlation results more tricky and requires employing multiconfigurational methods. Valence Bond (VB) theory, as well as Multiconfiguration Self Consistent Field (MCSCF) theory, are able to capture the resonance between different configurations and then they can be used to estimate nondynamic correlation. Worth mentioning, in this regard, are the seminal works on the dissociation of the hydrogen molecule of Coulson and Fischer [2] and of Heitler and London [3]. Among the various multiconfiguration techniques, the Complete Active Space SCF (CASSCF) is the most appropriate to fix a reference for the estimate of the nondynamical correlation energy. In such a case, the full valence active space should be taken into account in order to consider all possible arrangements of electrons on the different molecular orbitals. Unfortunately, such large active spaces require computation resources rapidly increasing with the number of active electrons and, consequently, this kind of CASSCF cannot be performed on medium-large molecular systems. Moreover, despite the CASSCF wave function offers a proper treatment of nondynamic correlation, it is usually not able to reproduce dynamic correlation effects unless a large amount of orbitals are included in the calculations. From a computational point of view, density functional theory (DFT) could be an efficient alternative to multiconfigurational methods because it assures a good balance between accuracy and computational cost. However, although DFT usually is considered able to recover most of dynamic correlation, common density functional approximations (DFA) are generally unable to take into account nondynamic correlation effects or they can reproduce them in an unspecified manner [4]. This limitation is related to the intrinsic single determinant shape of the Kohn-Sham (KS) electronic density. Indeed, in order to preserve the spin symmetry of the electronic Hamiltonian, the wave function employed in HF or DFT computations has to be constructed in a restricted framework. However, this formulation results inadequate for the accurate reproduction of ground state energy of systems with a multiconfigurational character. One way to overcome this problem is to relax such restriction by allowing electrons to occupy different orbitals for different spins. Such unrestricted approach allows a lowering of the final energy mimicking in this case the effects of nondynamic correlation. The resulting determinant will not be anymore an eigenfunction of the spin operatorŜ 2 . However, it is possible to evaluate the spin contamination between different spin states through the evaluation of the expectation values <Ŝ 2 >. It was shown that this approach allows including correlation for some particular systems, but its results are inadequate for others (see, for example, [5][6][7]). In the case of spin contamination we speak of symmetry breaking and in KS theory we can interpret the difference between restricted and unrestricted energy as an estimate of nondynamic correlation energy, namely To evaluate the extent of correlation that is included in a computational method, several tools based on an analysis of the wave function have been developed. Worth mentioning in this context are the T1 [8][9][10], D1 [11] and D2 [12] diagnostics which were specifically designed to identify nondynamic correlation signatures in Coupled Cluster wave functions. Other indicators of correlation are based on the analysis of the first order reduced density matrix (1-RDM), such for example the ones based on the quantum entanglement and its connection to entropy [13][14][15]. Particularly relevant indicators are the so-called I D and I ND which are obtained through a decomposition of the Coulomb hole into dynamic and nondynamic contributions [16]. These indicators admit a global [17] and a local [18] expression and can be evaluated by computing the natural orbitals in the framework of the method employed. Therefore these indicators can be easily computed for each method which has a non-diagonal 1-RDM. Further indicators based on natural orbitals have been developed by Tischenko et al. [19] and by Grimme and Hansen [20] in the context of Finite Temperature DFT (FTDFT) [20]. Such a thermal DFT approach treats statistical ensemble of quantum states describing the thermodynamic equilibrium of an electronic system [21]. It has been shown that this approach leads to fractional occupation numbers that can be used to define a sound measure of static electronic correlation. In principle, none of these indicators could be employed in DFT because the corresponding 1-RDM comes from a single determinant wave function and therefore it is diagonal on the basis of the molecular orbitals. On the other hand, if in the unrestricted framework a breaking in the spin symmetry occurs, through the Pulay localization [22] it is possible to recover a non-diagonal 1-RDM for a limited number of orbitals and studying the correlation effects.
However, in order to study these effects it is necessary to have a proper reference which appropriately takes into account the correlation. Among all the correlated methods in computational chemistry, the Quantum Monte Carlo (QMC) and in particular its variational variant (VMC) is one of the most efficient for the inclusion of both dynamic and nondynamic correlation effects. VMC is based on a stochastic integration over the expectation value of a many-body wave function, which contains several variational parameters optimized through stochastic techniques. The advantage of this method with respect to standard highly accurate wave function methods, such as Full Configuration Interaction (FCI), is that it guarantees a more favorable computational scaling with respect to the size of the molecule. Indeed it has been largely used in the last 30 years for the evaluation of energies and molecular properties (see, for example, [23]). The correlation effects of a VMC calculation method could be roughly estimated by comparing the energy of the system with the HF and CASSCF results; however, for a deeper analysis, it is necessary to compute the 1-RDM. The direct computation of the 1-RDM needs the performance of N particle integrals (where N is the number of electrons) of the wave function involving non local operators. The peculiar shape of the VMC wave function makes the computations of the 1-RDM by stocastic integration very challenging. In the past, Kent et al. [24] recovered the 1-RDM of silicon by resorting to a QMC wave function in which the determinantal part was a single determinant of local density approximation orbitals.
In this work, we present a route to estimate a correlated 1-RDM from a VMC electronic density through a stochastic fit procedure. The electronic density is one of the most relevant properties which could be reproduced through the VMC method [25], indeed it has been widely shown that accurate electronic densities could be computed by an efficient sampling of the positions of the electrons in the space employing different kinds of density estimators [26,27]. The trial 1-RDM employed for the fit needs to have a proper set of guess orbitals {φ i } and occupation numbers {n i } which are extracted from a test wave function able to preserve the N-representability conditions [28]. Indeed, it has been shown that is fundamental to satisfy this constraint in order to define a functional of 1-RDM [29][30][31]. The reference wave function has the skeleton's core orbitals coming from a single determinant calculation (DFT/B3LYP), while the highest occupied and the lowest virtual orbitals will be parameterized following different multiconfigurational strategies (GVB and CASSCF like wave functions). This wave function will be proved to be flexible enough to recover most of the correlation (dynamic and nondynamic) effects of the VMC electronic density. Further information on the method proposed will be reported in Section 2.1. The resulting 1-RDM will be used to analyze the correlation effects of some selected systems for which nondynamic correlation is predominant. Once the reliability of this method has been established, the results obtained will be employed as a reference for studying nondynamic correlation of unrestricted DFT computations.
The outline of the paper is as follows. In Section 2, we illustrate the methodology followed in this work based on the construction of the 1-RDM from VMC and from unrestricted Kohn-Sham computations and the subsequent comparison and analysis including a study on nondynamic correlation effects made by finite temperatures extracted from a Fermi-like distribution. In Section 3, we show the details of the calculations performed in this study and, in Section 4, we present the results obtained on a set of illustrative examples with the related discussion. Finally, in Section 5, we discuss the conclusions and possible developments for future work.

Methodology
Through the analysis of the one-particle reduced density matrix (1-RDM), it is possible to identify the amount of nondynamic correlation included in the calculation. In this work, we follow two alternative ways to reconstruct a correlated 1-RDM based on different methods which are both able to account for a certain amount of nondynamic correlation effects. In the first of the two routes, we have generated numerically a correlated electronic density by means of a large set of electronic configurations obtained from a VMC sampling. This density has been transcribed in terms of an N-representable 1-RDM that employs both occupied and virtuals KS orbitals through a multi determinantal wave function form.
Thanks to the quality of the VMC method, we expect that the resulting 1-RDM includes a considerable amount of both dynamic and nondynamic correlation contributions. The second route explored in this study, instead, involved unrestricted DFT (UDFT) computations. In this case we assumed that it is possible to include nondynamic correlation effects by breaking the spin symmetry of the system. Thus, the projection of the resulting electron density matrix on a natural orbital basis set will allow the study of the effects of correlation. Several density functional approximations (DFAs) have been considered to analyze how nondynamic correlation is treated by unrestricted UDFT in the different rungs of the Jacob's ladder [32]. The correlation effects have been studied by computing the global indicators of correlation developed by Ramos-Cordoba and coworkers [17], namely (3) where n σ i is the occupation number of the spin-natural orbital i with spin σ. Finally, the calculated natural orbital occupancies will be related to canonical KS single particle energies through a Fermi-like distribution.

1-RDM from a Correlated Electron Density
A correlated electron density can be generated from an accurate wave function obtained through a QMC calculation. In the present study, such a wave function has a multiconfiguration Slater-Jastrow form [33] and has been optimized at VMC level. More informations on the shape of the wave function are reported in Section 1 of the Supporting Information (SI). The optimized wave function Ψ QMC ( R), in which R is a 3N-dimensional vector which includes the coordinates of the N electrons, is used to efficiently sample the configuration space through the probability density Through a Metropolis importance sampling [34] it is possible therefore to generate an arbitrary number of electronic configurations. The electron density has been reconstructed on a grid employing an appropriate Gaussian smearing technique applied to a sufficiently large number of such configurations, namely where M is the number of configurations and ρ k is the value of the electron density for the k-th configuration. The expression of ρ k has been chosen as a combination of peaked Gaussian functions centered in the position of each electron and simulating as many delta functions, namely We expect that, with a sufficiently large number of configurations, the electronic density is well reproduced by Equation (6) and that the resulting density brings the information on the correlation effects of the VMC computation.

The Fit Procedure
The information on correlation could be recovered from an appropriate 1-RDM which is able to reproduce the correlated density of Equation (6). Such a 1-RDM has to be flexible enough to capture the basic features of the VMC electron density. We decided to construct a trial N-representable 1-RDM from a multiconfigurational reference wave function. In our construction we used a well given set of orbitals and we varied the coefficients of each determinant of such a function in order to fit the VMC density on diagonal. Therefore, the expression of the trial 1-RDM would be then where the matrix elements ρ ij depend on the coefficients of the determinants C k , namely the parameters of this fit. As many orbitals we include in Equation (8) as more accurate will be the fit itself. In this work, we used N el orbitals, where N el is the number of electrons. Finally, the 1-RDM has been parameterized combining a Complete Active Space Configuration Interaction (CASCI) approach for the π-orbitals included in the sample and a Generalized Valence Bond (GVB) one for the rest. In principle, several trial wave function could be employed for the parameterization of the 1-RDM; however, we believe that this kind of function is flexible enough to capture the main correlation contributions for the systems considered in this work. The single particle π-orbitals were taken from a CASSCF computation, while the others have been selected from a set of localized orbitals (bonding and antibonding type) obtained from a KS-DFT computation. After the fit of Equation (6) through Equation (8), we diagonalized ρ ij to obtain a set of natural orbitals {φ NO i } and the corresponding occupancies {n i } which will be used to analyze nondynamic correlation effects. The details of the procedure are illustrated in Section 3.

1-RDM from Unrestricted KS Computations
An unrestricted DFT calculation consists of solving the KS equations using a different set of orbitals for electrons with different spins. Thus, the KS determinant of the model system will not be anymore an eigenfunction of the spin operatorŜ 2 bringing, in principle, spin contamination in the calculation. Nevertheless, the resulting electron density should be again an acceptable N-representable density for the real interacting system. Therefore, the unrestricted KS determinant will have enough flexibility to include nondynamic correlation effects which are, otherwise impossible to be treated by restricted DFT. We can analyze the importance of these contributions recovering the 1-RDM of this computation appropriately computed following the Pulay localization scheme. The derivation of this 1-RDM is shown below. The electronic density of the system can be split into the sum of the spin up (α) and spin down (β) contributions where φ α i ( r) are the orbitals assigned to spin up electrons, while φ β i ( r) are the corresponding ones with spin down component.
The φ σ i ( r), as usual, are obtained from the SCF procedure in a basis of atomic orbitals: To compute the natural orbitals we project this density into the atomic orbital basis orthogonalized by the Löwdin procedure [35] denoted as {χ k }. The corresponding matrix expression of the 1-RDM would bē Using Equation (10), Equation (11) becomes where S ij are the elements of the overlap matrix constructed from the initial atomic orbitals, namely S ij = χ i |χ j . The diagonalization of Equation (12) leads to natural orbitals φ NO i ( r), and fractional occupation numbers n i . Thus, the electronic density can be rewritten as in which the occupation numbers n i lie between 0 and 2. Because these occupation numbers depend on the amount of spin contamination, we expect to observe different results for different DFAs and we will compare the occupancies with the one obtained through the QMC computations. Concerning the selection of the DFA, it is widely known that hybrid functionals with higher percentage of explicit HF exchange show bigger spin contamination. Indeed GGA functionals and hybrids functionals with a low percentage of explicit exchange are expected to bring less spin contamination with respect to range separated (RS) functionals. On the other hand, UHF calculations tend to largely overestimate symmetry breaking effects. In this work, we tested several density functional approximations focusing in particular on the different percentage of explicit HF exchange.

Fit of the Fermi Distribution
The occupation numbers previously introduced can be related to the corresponding canonical single particle energies through a Fermi-like distribution. The parameterization of such Fermi distributions should be made, in principle, by resorting to an effective temperature T e f f , by taking as chemical potential (µ e f f ) the average between the HOMO and LUMO energies, namely µ e f f = H + L 2 . In this work, we decided to include the chemical potential as a further parameter in order to increase the accuracy of the fit. Therefore, the function used for the fit is and the fit has been carried out on the occupation numbers calculated in both of the two approaches described in the previous sections. The effective temperature resulting from the fit is commonly considered as a further measure of nondynamic correlation. The single particle energies (˜ i ) have been evaluated as the expectation values of a generalized Fock operator over the set of natural orbitals so far calculated, namelỹ where the operatorF isF HF =ĥ +Ĵ −K (16) in whichĥ is the monoelectronic operator, containing the kinetic energy and the interaction with the external potential, andĴ andK are, respectively, the Coulomb and the exchange operators defined in terms of the 1-RDM. Equation (15) takes the form in which the integrals are evaluated in the natural orbital basis set and their expressions are:

Construction of the Electronic Density
To fit numerically the electron density, we interpolated the expression of ρ( r) on a regular cubic grid. The value of the density at each point of the grid has been taken as in which the exponent w of each Gaussian function has been considered as depending on the cell of definition. More precisely, w is proportional to the number of sampled electron positions in the given cell as follows where N c ( r α ) is the number of sampled configurations, and d is the side of the cubic cell. To make the electron density profile more homogeneous, the interpolation has been performed a few times up to convergence. At each iteration, the value of the density has been updated through a weighted average. Therefore, the electronic density at the iteration k (ρ k ) will have the following expressionρ for a density function ρ( r) extracted from a QMC run performed within the iteration itself.

Description of the Fitting Procedure
The fit has been performed with the least squares method, minimizing the quantity ∆ 2 , which is the square of the difference between the VMC electron density ρ QMC interpolated by Equation (19) and the diagonal of the trial 1-RDM (ρ T ) integrated over the whole space, namely The correct parameterization of the 1-RDM is of paramount importance for an accurate success of this fit. To reproduce appropriate values for the occupation numbers and kinetic energy it is necessary to impose precise N-representability constraints on this 1-RDM. For this reason, we decided to generate the 1-RDM by starting from a GVB wave function and, in a second stage, by modifying the parameterization of a reduced set of orbitals to include some larger active space taken from a CASSCF wave function. The starting 1-RDM has been created using bonding and antibonding orbitals obtained by the localization of B3LYP Kohn-Sham orbitals. We choose B3LYP orbitals for two main reasons, (i) they result from a computation in which electronic correlation is subsumed in the KS potential and (ii) B3LYP is often used to provide orbitals for standard QMC calculations. However, we believe that there is some flexibility in choosing the energy density functional for this purpose. This picture is appropriate for two-electron two-centre bonds while for the delocalized π electrons we improved the description by resorting to a complete active space properly defined. We chose a 1-RDM from a GVB wave function because it has a multiconfigurational character and the density matrix is diagonal for the representation in terms of bonding and anti-bonding orbitals. Moreover, the GVB wave function uses N el active orbitals, enough to capture the effects of nondynamic correlation. In the case of a perfect pairing, the form chosen in this study, the number of these configuration functions would be N * = 2 n p where n p is the number of pairs [36]. When the GVB geminals are written in terms of bonding and antibonding orbitals (see, for example, [37]), the GVB 1-RDM could be expressed in a diagonal form with the occupation numbers which depends strictly on 2 n p parameters, namely where the configuration coefficients C k depend on the order of excitation and the φ i are the bonding and antibonding orbitals. We apply such GVB scheme to the skeleton of the molecules considered in this work while we use a different active space for the remaining electrons. The two relevant subspaces are strong-orthogonal and the final 1-RDM can be written as a sum of two contributions, namely where ρ N a CAS is the 1-RDM parameterized as in a CASCI wave function with N a orbitals, while ρ N el −Na GVB is the remaining part of the 1-RDM. The expression of ρ Na CAS is indeed The 1-RDM ρ Na a ij has to be diagonalized in order to obtain the set of natural orbitals and occupancies. The minimization of ∆ 2 has been achieved by exploring the full space of parameters with a Metropolis sampling in order to find the deepest minimum and finally by a steepest descent to find the minimum.

Methods and Softwares
In this work, within the QMC method, we used VMC to optimize the wave function taken in a Slater-Jastrow form for all the systems. In all cases we resorted to the Burkatzki, Filippi and Dolg (BFD) pseudopotentials [38] with their associated valence triple-zeta (VTZ) Gaussian basis set. The Jastrow factor included always e − n and e − e two-body terms and e − e − n three-body terms [33]. The determinantal part of the wave function had a basic CASSCF structure. We used the suite of programs CHAMP [39] for all QMC calculations. The program GAMESS-US [40] is also used for the setup of the QMC wave function. All DFT and UDFT single point calculations have been performed using Gaussian 16 [41] employing the same pseudopotentials and basis sets used for the QMC calculations. Several functionals from the Jacob's ladder of DFT [32] have been tested, these functionals are: BLYP [42], B3LYP [43], BH&HLYP [44], M06-2X [45], CAM-B3LYP [46] and LC-BLYP [47]. The interpolation of the densities on the grids, the construction of the multiconfigurational parameterization, the numerical fitting of the density as well as the extraction of single particle canonical energies from UKS orbitals have been performed through homemade computer programs.

Fit of the Electronic Density
To test our approach, we considered a set of eight molecules which present a multiconfigurational character. This set includes two triatomic system: the ozone molecule and the Be + H 2 complex in its transition state structure. Moreover, we included five cyclic compounds: the cyclobutadiene molecule, studied in both D2h and D4h geometries, the benzyne ortho-, meta-and para-isomers. Finally, we included in this set the retinal protonated Schiff base chromophore with four conjugated π bonds (PSB4). A scheme with all the molecules employed in this study has been collected in Figure 1. These molecules display a different multiconfigurational character and therefore a different amount of nondynamic correlation. The geometries for the first seven molecules have been taken from the work of Grimme and Hansen [20] while, for the PSB4, we have taken the geometry from the work of Zulfikri et al. [48]. We report the results of the energies computed with three different wave function methods: HF, CASSCF, in which the active space has been chosen in order to include all the orbitals and electrons of the valence shell, and VMC. These total energies are collected on Table 1. From this table, we notice that the CASSCF energy is always between the HF and VMC values recovering from 7 to 18 percent of correlation in all cases except for the transition state of BeH 2 for which, instead, CASSCF and VMC are very close (98 percent). In these cases, the multireference character of the system is strongly related to the fraction of VMC correlation energy recovered by CASSCF. In Figure S1 of the SI we reported the electronic densities resulted from the QMC sampling. In Table 2, we report the results of the density fit in order to assess the reliability of the method presented in this work. In this Table, we compare the kinetic energy derived from the occupation numbers of the fitted one-particle reduced density matrix with the single particle KS kinetic energy and the correlated Pandharipande-Bethe kinetic energy derived from QMC calculation. The expression of the kinetic energy functional is known exactly and it is explicitly written in terms of the 1-RDM. As expected, the K f it value is in between KS and QMC recovering from about 20 to about 90 of the difference. The fit is not able to recover exactly the QMC value but goes beyond the single-particle kinetic energy of KS. This result is not trivial and we believe that it supports the use of the fitted one particle density matrix in the present study. Obviously there is room for improvement especially by switching from GVB to a more flexible functional form for the mapping of the QMC wave function. Finally, in the last column of Table 2 we report the final values of ∆ 2 which show that the accuracy is somewhat dependent on the size of the system. Table 2. Report on the accuracy of the fit and comparison between the kinetic energies. The number of digits for the kinetic energy depends on the accuracy of the method of calculation and reflects the statistical error of VMC and the standard error of the fit.

Molecule
Kinetic At this point, our fitted 1-RDM can be used to compute global and local indicators of correlation introduced by Ramos-Cordoba and Matito [17,18]. The values of global indicators for CASSCF and the fitted occupancies are reported in Table S1. For further information to the reader, we have also calculated the local indicators of correlation, these are displayed in Figures S2 and S3 of the SI. In Figure 2, we show a comparison between the results obtained from the QMC fit and the CASSCF calculations. More precisely, we compare the two fractions I ND /I TOT and I D /I TOT for the two sets of calculations. From this Figure, we can see that the eight systems cover a wide range of the intervals [0:1] of these fractions. Moreover, the QMC fit recovers essentially the nondynamical CASSCF correlation and adds some more dynamical correlation.

Unrestricted Computations
The breaking of the spin symmetry is strongly system and functional dependent. For the systems in which nondynamic correlation is predominant (Be + H 2 , cyclobutadiene (D4h), m-benzyne and p-benzyne) all the functionals considered in this study were able to break the spin symmetry leading to a lower unrestricted KS energy. However for systems in which nondynamic correlation effects are moderate (O 3 , cyclobutadiene (D2h), o-benzyne and PSB4) the results differ from case to case. In particular, for O 3 and PSB4 none of the density functionals tested was capable of breaking the spin symmetry. The same happens for the molecules cyclobutadiene (D2h) and o-benzyne employing UBLYP and UM06-2X. Instead, the GGA and RS functionals which have a BLYP core approximation broke the symmetry for all the other six molecular systems. UHF calculations introduced the largest spin contamination for all the systems except for O 3 molecule. Finally, we wanted to compare the correlation indicators of these unrestricted calculations with the corresponding values coming from the QMC fit. Pulay localization can be used to construct the unrestricted natural orbitals (UNO) which can be used to compute I ND . All global indicators relevant to unrestricted computations are reported in Table S3 of the SI. In Table 3 Table 3, (b) we confirm what already observed, indeed, also for relative values UB3LYP performs poorly while the other functionals give optimal results. Interestingly, we can observe that for UHF the relative indicators I ND /I TOT excellently agree with the reference results even if the global values of I ND are very overestimated. Overall, ULC-BLYP gives the best performance among the cases treated in this work.

Fit of the Fermi Distribution
We performed the fit of the Fermi distribution of the canonical single particle energies for the QMC results and for the different functionals following Equation (14). The results highlight that there is a good correlation between the energies and the occupation numbers for both the QMC and UDFT calculations. The fitted distributions obtained for the LC-BLYP functional for the molecules o-, m-and p-benzyne have been reported in Figure 3, while the other plots have been reported in in Figures S4-S10 of the SI. In the figure n i are the natural occupancies of the UNOs, whileε i are the canonicalized single particle orbitals. In black is reported the Fermi distibution calculated using Equation (14) and the parameters collected in Table S4 of the SI.
The parameters of the fit and the standard deviations have been reported in Table S4 of the SI. In Table 4, we compare the resulting effective temperatures T e f f coming from each functional with respect to QMC-fit reference. All the tested functionals give similar results for the systems for which nondynamic correlation is predominant, with a discrepancy in the range between 0 and 15% with respect to the reference T QMC . However these methods differ and give worse results for the systems in which nondynamic correlation is low. In particular, for o-benzyne all the resulting effective temperatures are significantly lower than the reference. UHF generally overestimate T e f f for all the systems tested. Furthermore, in this case, the best results are obtained with the ULC-BLYP approach (MAPE 11.1 and MAE 2.9 × 10 3 ).

Conclusions
The objective of this study has been to present a method able to recover a 1-RDM derived from an accurate wave function and then studying the nondynamic correlation effects included by unrestricted KS-DFT calculations. As accurate wave function, we mean a function of at least multireference configuration interaction quality. For this purpose, we have chosen a VMC wave function which shows a better scaling with the size of the system under study.
We propose a novel method to reconstruct the 1-RDM of a molecule at QMC level avoiding computing the challenging N particle integrals of the QMC wave function. This method consists of fitting the QMC density, properly defined through a large enough number of electronic configurations, using the density produced by a correlated wave function of which the 1-RDM is easily obtainable. In this way, in some sense, we move from the ideal non-interacting KS system toward the physical one along a path in which the density is the same but the reference wave function is modified to include, although not exhaustively, the e-e interaction. As mapping wave function, we chose a function that has a GVB structure for the skeleton of the molecule and a CASCI structure for the π electrons. The parameters of this fit are the coefficients of the Slater determinants of the mapping function. The subsequent analysis of the 1-RDM obtained in this way showed that it is possible to capture the main contributions related to both dynamic and nondynamic electronic correlation of the QMC wave function. To assess the reliability of the this procedure, we compared the resulting kinetic energy with the corresponding VMC value.
The 1-RDM extracted from the QMC approach above has been considered as a reference for comparison with the 1-RDM derived from unrestricted KS-DFT calculations. The spin broken symmetry that occurs in unrestricted calculations of an independent particle model system is due in part to static correlation. The associated 1-RDM can thus be obtained from a natural bond analysis such as the one we describe in Section 2.2. We applied this study on a selected set of molecules which have an electronic ground state with a non negligible multireference character. We noticed that broken symmetry KS computations, performed with density functionals with an high percentage of explicit exchange, are able to reproduce the global nondynamic correlation index I ND [17] computed through the QMC fitted 1-RDM. In particular, the best agreement occurs with the functional ULC-BLYP.
As a further test, we reconstructed a Fermi like distribution for each system and approach based on occupation numbers and canonicalized orbital energies. In each one of these distributions we have calibrated the effective temperature and chemical potential needed for the best fit. Although it is necessary a deeper investigation to understand a real meaning of such parameters, we notice that there is a general agreement between the effective temperature resulting from UDFT calculation with the QMC reference for the systems in which nondynamic correlation is important. For the other systems, the discrepancies are more evident. Again, the best functional among the ones included in the set is ULC-BLYP.
As for future work, we intend to improve the form of the mapping function to derive the 1-RDM from the QMC calculation. The need for a large active space, essentially a full valence active space, leads to a method that does not scale exponentially in the computation time with the size of the system under study. In this regard, the Density Matrix Renormalization Group (DMRG) method (see, for example, Ref. [49]) is a valid alternative to our GVB plus CASCI approach. DMRG method, in fact, is able to achieve a polynomial scaling under well defined conditions of calculation. As a final remark, we think that the present suggestion, namely the comparison of QMC and KSDFT through modeled 1-RDM, could be useful in the design of new functionals for DFT and FTDFT calculations.