λ-Density Functional Valence Bond: A Valence Bond-Based Multiconfigurational Density Functional Theory With a Single Variable Hybrid Parameter

A new valence bond (VB)-based multireference density functional theory (MRDFT) method, named λ-DFVB, is presented in this paper. The method follows the idea of the hybrid multireference density functional method theory proposed by Sharkas et al. (2012). λ-DFVB combines the valence bond self-consistent field (VBSCF) method with Kohn–Sham density functional theory (KS-DFT) by decomposing the electron–electron interactions with a hybrid parameter λ. Different from the Toulouse's scheme, the hybrid parameter λ in λ-DFVB is variable, defined as a function of a multireference character of a molecular system. Furthermore, the EC correlation energy of a leading determinant is introduced to ensure size consistency at the dissociation limit. Satisfactory results of test calculations, including potential energy surfaces, bond dissociation energies, reaction barriers, and singlet–triplet energy gaps, show the potential capability of λ-DFVB for molecular systems with strong correlation.


INTRODUCTION
One of the major interests in quantum chemistry is the methodology development for electronic correlation energy calculation with an affordable computational cost. The basic multiconfigurational wave function methods, the multiconfigurational self-consistent field (MCSCF) method Siegbahn et al., 1980), and the valence bond analog, valence bond self-consistent field (VBSCF) method Balint-Kurti, 1980, 1983), mainly consider static electron correlation, which is not covered in a single configuration-based wave function. Based on multiconfigurational wave function, the perturbative theory (PT), coupled-cluster (CC), or configuration interaction (CI) can be employed to cover dynamic correlation. With these post-self-consistent field (SCF) techniques, many high-level multiconfigurational wave function methods are proposed, including the multireference perturbation theory complete active space second-order perturbation theory (CASPT2) (Andersson et al., 1992), multireference second-order Møller-Plesset perturbation theory (MRMP2) (Nakano, 1993), valence bond second-order perturbation theory (VBPT2) (Wu et al., 1998;Chen et al., 2009), multireference configuration interaction (MRCI) (Siegbahn et al., 1981), valence bond configuration interaction (VBCI) (Wu et al., 2002;Song et al., 2004), breathing-orbital valence bond (BOVB) (Hiberty et al., 1992(Hiberty et al., , 1994, and so on. The computational costs of these post-SCF multiconfigurational wave function methods increase rapidly with the size of active space. Meanwhile, owing to its high efficiency for dynamic correlation calculation, the Kohn-Sham density functional theory (KS-DFT) is the most widely used electronic structure method (Hohenberg and Kohn, 1964;Kohn and Sham, 1965;Parr and Yang, 1989;Koch and Holthausen, 2001). Thus, the development of the multireference wave function-based DFT (MRDFT) method, in which the dynamic correlation is considered by DFT functionals and the static correlation is covered by a multiconfigurational wave function method, is promising because of its economical computational cost. Since the 1990s, many MRDFT schemes have been proposed (Lie and Clementi, 1974;Miehlich and Stoll, 1997;Filatov and Shaik, 1999;Gräfenstein and Cremer, 2000;Grafenstein and Cremer, 2005;Head-Gordon, 2003;Gusarov et al., 2004;Pérez-Jiménez et al., 2004;Yamanaka et al., 2006;Wu et al., 2007;Cembran et al., 2009;Kurzweil et al., 2009;Rapacioli et al., 2010;Sharkas et al., 2012;Ying et al., 2012;Manni et al., 2014;Zhou et al., 2017), most of which use the MCSCF/ complete active space self-consistent field (CASSCF) as the multireference wave function. Recently, two valence bond wave function-based MRDFT (DFVB) methods are presented: the first one is the dynamic correlation-corrected density functional valence bond (dc-DFVB) method (Ying et al., 2012), and the second one is the Hamiltonian matrix correction-based density functional valence bond (hc-DFVB) method (Zhou et al., 2017). These two methods are capable of providing satisfactory accuracy with relatively cheaper computational costs, compared to the currently existing post-VBSCF methods. However, these two methods still suffer from double counting error (DCE).
DCE is one of the key issues for MRDFT, because it is impossible to separate the static and dynamic correlations exactly. The range-separated scheme is considered to be helpful for MRDFT to avoid double counting error (Fromager et al., 2007(Fromager et al., , 2009). In the range-separated MRDFT, electron-electron interaction operator is decomposed into two components: the long-range term considered by the wave function method and the short-range term described by a density functional approximation. Recently, a multiconfigurational hybrid density functional theory, which is called multiconfigurational oneparameter hybrid (MC1H) approximation, is proposed by Sharkas et al. (2012). MC1H is based on a linear decomposition of electron-electron interactions with a hybrid parameter λ. It is shown that the accuracy of this multiconfigurational hybrid scheme matches that of the range-separated multiconfigurational hybrid method (Sharkas et al., 2012).
To remove the DCE from the DFVB methods, this paper presents a new VB-based MRDFT method, named λ-DFVB, by utilizing the MC1H scheme. Different from the MC1H scheme, which suggests setting λ as 0.25, the value of λ is variable in λ-DFVB, defined as a function of multireference character. Furthermore, the energy expression is modified to consider the dynamic correlation energies of dissociated fragments/atoms, ensuring the size consistency at the dissociation limit.

METHODOLOGY
In VB theory, the many-electron wave function Ψ is expressed as a linear combination of VB functions (Hiberty and Shaik, 2008;Wu et al., 2011;Su and Wu, 2013), where Φ K and C K are VB functions corresponding to a specific structure and its coefficient, respectively. In spin-free quantum chemistry, VB function, which is an eigenfunction of spin and antisymmetric with respect to permutations of electron indices, is of the form where Â is an antisymmetrizer for electron indices, Ω 0 is an orbital product, and K is a spin eigenfunction (Pauncz, 1979), defined as In Equation 4, spin pairs (k 1 , k 2 ), (k −3 , k 4 ), etc., correspond to covalent bonds in structure K, and k p is for unpaired electrons. Coefficients {C K } in Equation 1 can be obtained by solving the secular equation: where H, M, and C are Hamiltonian, overlap, and coefficient matrices, respectively. In a similar fashion to molecular orbital methods, there are various ab initio classical VB methods Balint-Kurti, 1980, 1983;Hiberty et al., 1992Hiberty et al., , 1994Hiberty and Shaik, 2002;Wu et al., 2002;Song et al., 2004;Chen et al., 2009). Among them, the valence bond self-consistent field (VBSCF) method is the basic one. In VBSCF, both VB structure coefficients {C K } and VB orbitals {ϕ i } are optimized simultaneously to minimize the total energy E. VB orbitals are usually expanded as linear combinations of basis functions.
VB orbitals may be taken as strictly localized hybrid atomic orbitals (HAOs), semilocalized bond-distorted orbitals (BDOs) (Mo et al., 1994(Mo et al., , 1996, or delocalized overlap-enhanced orbitals (OEOs) (Bobrowicz and Goddard, 1977;Cooper et al., 1991), according to the specific purpose of a study. Analogous to the MCSCF/CASSCF method, VBSCF mainly considers static correlation. VB structural weights can be evaluated by the Coulson-Chirgwin formula (Chirgwin and Coulson, 1950), which is an equivalence of the Mulliken population analysis, VB structural weights are typically used to compare the relative importance of individual VB structures and can be helpful in the understanding of the correlation between molecular structure and reactivity.
In dynamic correlation-corrected density functional valence bond (dc-DFVB), the total energy can be expressed as (8) where E C [ρ] is obtained from a pure correlation functional with electronic density computed from the VBSCF wave function. The dc-DFVB method improves the VBSCF results, but it suffers from the double counting error due to the fact that it simply includes the total E C energy in the total energy.
In the MC1H approximation by Sharkas et al. (2012), the energy of the MRDFT can be determined by minimizing the following expression: where T, V ext , and W ee are the kinetic energy, external potential, and electron-electron interaction operators, respectively; E λ HXC [ρ]is the complement λ-dependent Hartree-exchangecorrelation density functional for electronic density ρ, and λ is a coupling parameter. E λ HXC [ρ]is defined as (Sharkas et al., 2012): where E X [ρ] and E C [ρ] are the exchange and correlation functionals, respectively. E H is the Hartree energy, given as Equation 10 turns to the KS-DFT formula when λ = 0.0, while it becomes wave function theory (WFT) if λ = 1.0. As suggested by Toulouse, the value of λ is approximately taken as 0.25. It is clear that the parameter λ indicates the hybrid extent of the WFT and the KS-DFT. Based on the fact that multiconfiguration-based WFT is suitable for molecules with multireference character, while KS-DFT is a good tool for molecules with single-reference character, it is more reasonable to allow the value of λ to be different with different molecules. In this paper, we use a variable parameter for λ, which represents the multireference character of a molecule instead of a fixed value. There have been various indices for estimating multireference character or static correlation character (Zeische et al., 1997;Janssen and Nielsen, 1998;Leininger et al., 2000;Zanardi, 2002;Huang et al., 2006;Sears and Sherrill, 2008a,b;Tishchenko et al., 2008;Ramos-Cordoba et al., 2016;Benavides-Riveros et al., 2017;Ramos-Cordoba and Matito, 2017;Rodriguez-Mayorga et al., 2017), A large diagnostic value indicates a strong multireference character. These diagnostics include the T1 and D1 diagnostics in coupled-cluster wave functions (Lee and Taylor, 1989;Janssen and Nielsen, 1998;Leininger et al., 2000), the M diagnostic (Tishchenko et al., 2008), the 1-C 2 0 diagnostic in CASSCF wave function (Sears and Sherrill, 2008a,b), the S 2 diagnostic (Zeische et al., 1997;Zanardi, 2002;Huang et al., 2006); and the I ND diagnostic (Ramos-Cordoba et al., 2016;Ramos-Cordoba and Matito, 2017), etc. In this paper, the concept of free valence is utilized to diagnose the multireference character of a molecule and is further used to determine the value of λ.
The free valence of an atom A, F A , is defined as where V A and O AB are the total valence of atom A and the bond order between atoms A and B, respectively, defined as (Mayer, 2003) and In Equations 13 and 14, S is the overlap matrix in terms of basis functions, D = P α + P β is the total density matrix, and P s = P α − P β is the spin polarization density matrix for basis functions, where P α and P β are α and β density matrices, respectively. The molecular free valence index K is defined as It is clear that K ranges from 0 to 1. The value of K is small at the equilibrium geometry because all the atoms are bonded. For example, at the equilibrium distance, the Mayer bond order of H 2 is 0.953 at VBSCF/cc-pVTZ. Then, F H1 = V H1 -0.953 = 0.047 (subscript H1 denotes the first hydrogen atom); K = (0.047 + 0.047)/(1.0 + 1.0) = 0.047. At the dissociation limit, K = 1.0, as O AB = 0 and F A = V A . Table S1 shows the comparisons of various diagnostics for diatomic molecules H 2 , HF, F 2 , N 2 , C 2 , and Cr 2 in their equilibrium geometries, respectively. A large diagnostic value indicates a strong multireference character. Although the various diagnostic values are much different, their trends are in good agreement, showing the validation of K.
In this paper, the hybrid parameter λ is expressed as a function of the free valence index K. Based on some numerical investigations, shown in the Supporting Information, the function is defined as Clearly, the λ also ranges from 0.0 to 1.0, which satisfies the requirement of the hybrid parameter. In Equation 10, the factor (1-λ 2 ) of E C [ρ] arises from the fact that WFT covers the λ fraction of correlation and thus should be deducted from the functional E C [ρ]. In λ-DFVB, VBSCF wave function is used for the WFT part, and thus, only the static correlation of valence electrons, which results from the use of multideterminants, is covered by WFT. Therefore, the corresponding E C functional should be approximately expressed as is the E C correlation energy determined by the electronic density of the leading determinant (LD), which is a single determinant with the largest coefficient in the VBSCF wave function. At a short distance, if dynamic correlation energy is dominating, E C [ρ]-E C [ρ LD ] is small. At a long distance, the difference can be large because static correlation is largest at the bond dissociation limit. At last, the λ-DFVB energy is expressed as where When a molecule composed of two fragments/atoms A and B is fully dissociated, λ = 1.0. Then, the total energy can be expressed as At the infinite distance, there is no overlap between fragments/atoms A and B; thus, the density of the leading determinant can be expressed as the sum of the densities (ρ A + ρ B ) of two fragments/atoms. As such, at the dissociation limit, E C [ρ LD ] takes the dynamic correlations of dissociated fragments/atoms into account. A λ-DFVB computation can be performed with the following steps: (1) Compute a VBSCF calculation to obtain the λ value and the VBSCF density ρ.

HXC
[ρ] is set into the VB Hamiltonian matrix. The computation is consistently iterative until convergence is achieved. (5) Obtain the λ-DFVB energy based on the wave function optimized in step 4:

COMPUTATIONAL DETAILS
The λ-DFVB method has been implemented in the Xiamen Valence Bond (XMVB) package (Su and Wu, 2013;Chen et al., 2015). All the VB calculations are performed by XMVB, while all KS-DFT calculations are carried out by General Atomic and Molecular Electron Structure System (GAMESS) (Schmidt et al., 1993;Gordon and Schmidt, 2005 Test calculations involve the potential energy surfaces of H 2 , HF, F 2 , N 2 , C 2 , and Cr 2 , the reaction barriers of the Diels-Alder (D-A) and Menshutkin reactions, and the λ-DFVB energy gaps of carbon atom, oxygen atom, carbene (CH 2 ), trimethylenemethane (TMM), and organometallics Fe(II)porphyrin. The geometry of Fe(II)-porphyrin is optimized at the B3LYP level. The geometries of TMM, carbene, and the reactants and the transition states for the D-A reaction and the Menshutkin reaction are taken from previous papers (Ying et al., 2012;Huang et al., 2014;Zhou et al., 2017).

RESULTS AND DISCUSSIONS
The λ Values With Various Truncation Levels of VB Wave Function First, the parameter λ is examined with various truncation levels of VBSCF wave function, including covalent structures only (denoted as COV, which is the 0th truncation level shown in Figure 1), covalent structures plus the 1st ∼ nth order ionic structures, . . . , and all structures (denoted as CAS). Clearly, the CAS levels for N 2 and C 2 are the 3rd and 4th truncation levels, respectively. The numbers of VB structures in various truncation levels of N 2 and C 2 are listed in Table S2. Figure 1 displays the λ values of diatomic molecules N 2 and C 2 using OEOs, which are fully delocalized over the whole molecules, with the various truncation levels of VB wave function. As can be seen, the curves are almost flat, showing that the value of λ is not sensitive to the truncations of wave function. Meanwhile, including the ionic structures tends to reduce the λ value. In general, C 2 has the larger λ value compared to N 2 because C 2 has the stronger multireference.   Table 1 displays the VBSCF and λ-DFVB energies of H 2 , F 2 , HF, N 2 , C 2 , and Cr 2 with the CAS and COV levels of the VB wave function. The number of active electrons and active orbitals (n e , n o ) are listed in the second column. The dynamic correlation correction of λ-DFVB, which is defined as the energy difference between VBSCF and λ-DFVB, is listed in the last column. The corresponding data for polyatomic molecules CH 3 -CH 3 , C 2 H 4 , and C 2 H 2 are shown in Table S3. In general, the λ values of CAS are smaller than those of COV. The molecules with strong correlation tend to have the large λ values. Among them, H 2 has the smallest dynamic correlation correction while Cr 2 has the largest one.
The λ-DFVB calculations in the next are carried out at the CAS level of the VB wave function. Potential Energy Surfaces of H 2 , HF, F 2 , N 2 , C 2 , and Cr 2 The calculation of potential energy surface for bond breaking is one of the most rigorous tests for electronic structure methods. Figure 2 shows the curves of the λ values of diatomic molecules along the potential energy surfaces. As can be seen, the λ value goes up with the increase in bonding distance and reaches to 1.0 at the dissociation limit. For example, the λ value of H 2 is 0.465 at the equilibrium bond distance (0.74 Å), while it is 1.0 at the dissociation limit. The other curves share similar behavior. It can be found that the molecules with strong correlation have large λ values in the short distances, showing the large portions of electron-electron interaction energy computed by the VB wave function methods. Based on the λ values shown in Figure 3, the bond dissociation curves of H 2 , HF, F 2 , N 2 , C 2 , and Cr 2 by λ-DFVB are plotted in Figure 3. For comparison, the PES curves with various WFT and KS-DFT methods are also shown.
It can be seen from Figure 3A that the restricted B3LYP and BLYP calculations for the H 2 molecule go to the wrong dissociation limits, as expected. Both VBSCF and dc-DFVB predict the dissociation correctly. However, VBSCF gives a higher energy at the equilibrium geometry due to the lack of dynamic correlation, while the dc-DFVB curve is underneath the full CI one because of the double counting error. Encouragingly, the performance of λ-DFVB is excellent, virtually overlapped with the full CI and CASPT2 ones along the whole curve. If one zooms  (Casey and Leopold, 1993) a Values in parentheses are the IPEA shifts used in CASPT2 calculations. Expt 33.0 (Webb and Gordon, 1999) 23.3 ± 2 (Webb and Gordon, 1999;Guner et al., 2003) in the figure, it can be found that the λ-DFVB curve is the closest to that of the full CI near the equilibrium geometry. For the HF molecule, shown in Figure 3B, the λ-DFVB curve is quite close to those by dc-DFVB, MRCI, and CASPT2. Interestingly, although the dc-DFVB curve virtually overlaps with the others around the equilibrium geometry, it deviates from the others at about 1.7 Å and converges again at the dissociation limit. For the F 2 results in Figure 3C, the VBSCF curve is the highest, while that of the dc-DFVB is the lowest after 1.60 Å. The λ-DFVB curve almost coincides with those of the MRCI and CASPT2 after 1.60 Å. However, the λ-DFVB curve is lowest around the equilibrium distance. For the N 2 curves, shown in Figure 3D, the performance of λ-DFVB is very close to MRCI and CASPT2, better than those of VBSCF and dc-DFVB. For the C 2 curves in Figure 3E, the three VB methods (VBSCF, dc-DFVB, and λ-DFVB) and CASPT2 predict the bonding dissociation of C 2 in a similar fashion. It is shown that the λ-DFVB curve is slightly higher than that of the VBSCF around the equilibrium bond length, indicating that the dynamic correlation does not play a key role in the relative energy surface.
It is well-known that Cr 2 , which has sextuple bond with a small bonding energy, is a challenging molecule to quantum chemical methods due to its strong multireference character. Traditional KS-DFT calculations are unable to provide the potential energy surface properly (Brynda et al., 2009). For CASPT2, different values of the ionization potential and electron affinity (IPEA) shift lead to different results (Ruipierez et al., 2011). Figure S1 displays the CASSCF curves with the Stuttgart RSC 1997 ECP basis set (abbreviated as ECP) and the ANO-RCC-VTZP basis set (abbreviated as ANO), and the VBSCF curve with the ECP. It is found that the three curves are quite similar, and all of them predict a minimum around 3.0 Å, far away from the experimental equilibrium distance. In Figure 3F, ECP is used for VBSCF and λ-DFVB, and ANO is used for CASPT2. As expected, CASPT2 is sensitive to the value of the IPEA shift (Ruipierez et al., 2011;Manni et al., 2014). The CASPT2 curve with a value of 0.45 a.u. is close to that of the experimental. It is interesting that λ-DFVB successfully predicts the global minimum around 1.60 Å. Meanwhile, the barrier around 2.8 Å is in agreement with the curves computed by the modified generalized valence bond (MGVB) method and the curve by the MC-PDFT method (Manni et al., 2014).

The Bond Dissociation Energies of Diatomic Molecules
The computed bond dissociation energies (BDEs; D e ) for the six diatomic molecules at their optimized geometries are shown in Table 2. The BDEs of BLYP and B3LYP are computed as the energy difference between the equilibrium bond distances and the sum of atomic energies. Generally speaking, CASPT2 and MRCI perform well, showing the small  (Li and Paldus, 2008) deviations from the experimental data. With a large IPEA shift of 0.45 a.u., CASPT2 provides a satisfactory result for Cr 2 . Both BLYP and B3LYP are unable to provide the proper descriptions for C 2 and Cr 2 molecules, as mentioned in literature (Carlson et al., 2015;Kepp, 2017).
Because of the lack of dynamic correlation, VBSCF is unable to provide satisfactory BDEs for H 2 , HF, F 2 , and N 2 . Analogous to CASSCF, VBSCF is unable to describe the Cr-Cr bonding properly, but it predicts the bonding of C 2 quite well. BOVB uses different orbitals for different VB structures to consider the dynamic electron correlation. It cannot be employed in the molecules with large active spaces (for example, C 2 and Cr 2 in this work). In the BOVB calculations, the active orbitals are HAOs, while the remaining orbitals are OEOs. It is shown that the BOVB results are better than those of the VBSCF. By incorporating E C functional into VBSCF, the BDEs of dc-DFVB are larger than their corresponding VBSCF values and are mostly overestimated for N 2 and C 2 compared to the experimental data, because of the DCE. Similar to VBSCF, dc-DFVB is also incapable of predicting the stable Cr-Cr bonding.
For λ-DFVB, in general, its computational results are excellent, close to the MRCI and CASPT2 values and superior to the BOVB and dc-DFVB ones. For N 2 , the λ-DFVB value is 224.3 kcal/mol, which is close to the experimental data of 228.5 kcal/mol and even better than the CASPT2 one. The λ-DFVB value of C 2 , 137.4 kcal/mol, is close to the MRCI result of 137.8 kcal/mol and that of ICMRCI+Q/cc-pVTZ, 138.8 kcal/mol (Pradhan et al., 1994). For Cr 2 , the BDE value of λ-DFVB is 38.7 kcal/mol, close to the experimental data, 33.9 kcal/mol, and the CASPT2 result of 32.8 kcal/mol with the IPEA shift of 0.45 a.u. and cc-pVTZ-DK basis set by Truhlar (Carlson et al., 2015), better than the MC-PDFT result, 13.8 kcal/mol, at the tPBE/cc-pVTZ-DK level (Manni et al., 2014). Table 3 lists the reaction barriers for the Diels-Alder reaction and the Menshutkin reaction. The λ values for the reactants and transition states of the two reactions are shown in Table S4. It can be seen that VBSCF overestimates the values of barriers due to the lack of dynamic correlation, while KS-DFT underestimates the reaction barriers, as expected. The dc-DFVB results are much improved from those of the VBSCF. The performance of λ-DFVB is very good. For the Diels-Alder reaction, the barrier of 24.6 kcal/mol is very close to those of the CASPT2 (23.5 kcal/mol) and the experimental (23.3 kcal/mol). For the Menshutkin reaction, the λ-DFVB value of 32.2 kcal/mol is even better than that of CASPT2, 40.5 kcal/mol, compared to the experimental value of 33.0 kcal/mol.

The Excitation Energy Gaps
The singlet-triplet energy gaps of carbon atom, oxygen atom, carbene (CH 2 ), and trimethylenemethane (TMM) by various methods are shown in Table 4, the corresponding λ values for the λ-DFVB calculations are shown in Table S4. As expected, the KS-DFT results show very large deviations from the experimental data and the CASPT2 results for all atoms and molecules except porphyrin. Meanwhile, VBSCF predicts the excitation and transition energies quite well. Compared to VBSCF and dc-DFVB, the performance of λ-DFVB is much improved, particularly for molecules. As can be seen, the deviation values from the experimental data are 1.6 and 0.3 kcal/mol for CH 2 and TMM, respectively, even smaller than those of CASPT2.
Fe(II)-porphyrin, shown in Figure 4, is an important organometallic compound comprising the active center of several important biological proteins. Experimental studies show that the triplet state is lower than the quintet state. However, it is a challenge to describe the relative stability of the quintet and triplet states correctly with multireference wave function methods. Manni and Alavi (2018) and Smith et al. (2017) found that only the CASSCF calculations with very large active spaces are able to predict the triplet-quintet gaps properly. For the VB and CASSCF calculations, the active space includes the six valence electrons of the metal center and its five 3d orbitals. Our computed triplet-quintet gaps for Fe(II)-porphyrin with various methods are displayed in Table 5. It is found that the VBSCF, dc-DFVB, CASSCF, and CASPT2 results predict that the quintet state is more stable than the triplet state. Only The λ-DFVB with the PW91 functional is also performed for BDEs, chemical barriers, and singlet-triplet gaps. In general, the results are close to those of λ-DFVB with BLYP. Details are shown in Tables S5-S7. The structure weights and orbitals of N 2 are displayed in Figures S2, S3 showing that the wave function of λ-DFVB is similar to that of the VBSCF.

CONCLUSION
A new hybrid multireference density functional theory method based on the VB theory, named λ-DFVB, is presented in this paper. Based on the MC1H approximation presented by Sharkas et al. (2012), λ-DFVB combines VBSCF and KS-DFT with a linear decomposition for electron-electron interactions. In λ-DFVB, the hybrid parameter λ is variable, ranging from 0.0 to 1.0, and defined as a function of the free valence index K, which diagnoses the multireference character for a given system. Furthermore, an additional correlation term, E C (ρ LD ), is introduced to consider the correlation energies of fragments/atoms in the dissociation limit, which ensures that the λ-DFVB method is size consistent.
The λ-DFVB method was carefully examined by performing test calculations for various chemical properties, including potential energy surfaces, bond dissociation energies, chemical reaction barriers, and singlet-triplet energy gaps. The performance of λ-DFVB is promising, close to those of CASPT2 and MRCI. Especially, the proper descriptions of the Cr 2 bonding and the triplet-quintet gap of the model molecule Fe(II)-porphyrin in their equilibrium geometries, which are challenging both to the WFT and DFT methods, show the capability of λ-DFVB for strong correlation systems.
The λ-DFVB method shares its dual advantage. On the one hand, λ-DFVB improves the accuracy of the VBSCF method by incorporating dynamic correlation; on the other hand, it overcomes some problems with KS-DFT that result from the use of a single determinant. Though the current strategy of the λ-DFVB is applied to the GGA functionals in this paper, it will be extended to more general functionals, such as hybrid functionals and meta-GGA functionals, which will be discussed in the near future.

AUTHOR CONTRIBUTIONS
FY: the implementation of λ-DFVB and numerical tests; CZ: the implementation of λ-DFV and numerical tests; PZ and JL: numerical tests; PS and WW: methodology and formula derivations.

ACKNOWLEDGMENTS
This project is supported by the Natural Science Foundation of China (No. 21503172, 21573176, and 21733008) and the New Century Excellent Talents in Fujian Province University. PS thanks Prof. Emmanuel Fromager of Strasbourg University and Prof. Andreas Savin of Sorbonne University for the helpful discussions.