λ-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 E C 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 (Roos et al., 1980; Siegbahn et al., 1980), and the valence bond analog, valence bond self-consistent field (VBSCF) method (van Lenthe and 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, 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, 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 one-parameter 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 Φ and C 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 Θ is a spin eigenfunction (Pauncz, 1979), defined as

In Equation 4, spin pairs (k1, k2), (k−3, k4), etc., correspond to covalent bonds in structure K, and k is for unpaired electrons.

Coefficients {C} 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 (van Lenthe and Balint-Kurti, 1980, 1983; Hiberty et al., 1992, 1994; Hiberty 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} and VB orbitals {φ} 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, 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

where EC[ρ] 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 EC 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, Vext, and Wee are the kinetic energy, external potential, and electron–electron interaction operators, respectively; is the complement λ-dependent Hartree–exchange–correlation density functional for electronic density ρ, and λ is a coupling parameter. is defined as (Sharkas et al., 2012):

where EX[ρ] and EC[ρ] are the exchange and correlation functionals, respectively. EH 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- diagnostic in CASSCF wave function (Sears and Sherrill, 2008a,b), the S2 diagnostic (Zeische et al., 1997; Zanardi, 2002; Huang et al., 2006); and the IND 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, is defined as

where V and O 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 = 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 H2 is 0.953 at VBSCF/cc-pVTZ. Then, FH1 = VH1-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 = 0 and F = V.

Table S1 shows the comparisons of various diagnostics for diatomic molecules H2, HF, F2, N2, C2, and Cr2 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 EC[ρ] arises from the fact that WFT covers the λ fraction of correlation and thus should be deducted from the functional EC[ρ]. 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 EC functional should be approximately expressed as EC[ρ]–EC[ρLD], where EC[ρLD] is the EC 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, EC[ρ]–EC[ρ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, EC[ρLD] takes the dynamic correlations of dissociated fragments/atoms into account.

A λ-DFVB computation can be performed with the following steps:

Compute a VBSCF calculation to obtain the λ value and the VBSCF density ρ.

Compute by Equation 18.

Compute the operator , which is defined as

Optimize the λ-DFVB wave function with the following equation:

where n(r) is the density operator, ρ(r) = 〈Ψ| n(r) |Ψ〉. The contribution of is set into the VB Hamiltonian matrix. The computation is consistently iterative until convergence is achieved.

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). Free valence and Mayer's bond order are computed with the VBSCF wave function. The MOLCAS 8.0 program (Aquilante et al., 2016) was used for CASSCF, MRCI, and CASPT2 calculations. The Davidson correction is considered for MRCI calculations. Two Generalized Gradient Approximation (GGA) functionals, Becke88 and Lee-Yang-Parr (BLYP) and Perdew-Wang 91 (PW91), are employed for the λ-DFVB calculations. The results of λ-DFVB with BLYP are shown in the main text, while those with PW91 are shown in the Supporting Information. In the dc-DFVB calculations, Lee-Yang-Parr (LYP) functional is used. For comparison, the corresponding results of B3LYP, BLYP, CASSCF, CASPT2, BOVB, and dc-DFVB are also provided.

Test calculations involve the potential energy surfaces of H2, HF, F2, N2, C2, and Cr2, the reaction barriers of the Diels–Alder (D-A) and Menshutkin reactions, and the energy gaps of carbon atom, oxygen atom, carbene (CH2), 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).

The cc-pVTZ (CCT) basis set was used for the potential energy surfaces of H2, HF, F2, N2, C2, and the energy gaps (except Fe(II)–porphyrin). 6-31G* was used for the two chemical reactions. For Fe(II)–porphyrin, Lanl2DZ is for the iron atom and 6-31G* is for the C, H, and N atoms. For Cr2, two basis sets, Stuttgart Royal Society of Chemistry (RSC) 1997 ECP (Andrae et al., 1990) and ANO-RCC-valence triple-zeta with polarization (VTZP) (Pou-Amérigo et al., 1995; Roos et al., 2005), were used.

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 N2 and C2 are the 3rd and 4th truncation levels, respectively. The numbers of VB structures in various truncation levels of N2 and C2 are listed in Table S2. Figure 1 displays the λ values of diatomic molecules N2 and C2 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, C2 has the larger λ value compared to N2 because C2 has the stronger multireference.

Figure 1

The λ values for N2 and C2 with various truncation levels of VB wave function.

Table 1 displays the VBSCF and λ-DFVB energies of H2, F2, HF, N2, C2, and Cr2 with the CAS and COV levels of the VB wave function. The number of active electrons and active orbitals (ne, no) 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 CH3-CH3, C2H4, and C2H2 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, H2 has the smallest dynamic correlation correction while Cr2 has the largest one.

Table 1

The λ-DFVB energies of H2, F2, HF, Cr2, N2, and C2 energies with variable λ values at their equilibrium geometries (a.u.).

	Active space		EVBSCF	Eλ−DFVB	λ	Ecorra
N2	(6,6)	COV	−75.589398	−75.924588	0.764	−0.335190
		CAS	−75.637301	−75.952588	0.728	−0.315287
C2	(8,8)	COV	−109.065922	−109.524697	0.560	−0.458775
		CAS	−109.120064	−109.540330	0.546	−0.420266
H2	(2,2)	COV	−1.151417	−1.173454	0.465	−0.022037
		CAS	−1.151419	−1.172552	0.465	−0.021133
F2	(2,2)	COV	−198.828556	−199.505642	0.736	−0.677086
		CAS	−198.828556	−199.506120	0.727	−0.677564
HF	(2,2)	COV	−100.081614	−100.450733	0.437	−0.369119
		CAS	−100.081618	−100.450755	0.436	−0.369137
Cr2	(12,12)	COV	−172.514493	−173.471875	0.911	−0.957382
		CAS	−172.598336	−173.588538	0.809	−0.990202

E.

The λ-DFVB calculations in the next are carried out at the CAS level of the VB wave function.

Potential Energy Surfaces of H2, HF, F2, N2, C2, and Cr2

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 H2 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 H2, HF, F2, N2, C2, and Cr2 by λ-DFVB are plotted in Figure 3. For comparison, the PES curves with various WFT and KS-DFT methods are also shown.

Figure 2

The curves of λ as functions of bond distances for diatomic molecules.

Figure 3

The PES curves of diatomic molecules with various methods: (A) H2; (B) HF; (C) F2; (D) N2; (E) C2; and (F) Cr2. The numbers in the brackets after “CASPT2” denote the IPEA shift.

It can be seen from Figure 3A that the restricted B3LYP and BLYP calculations for the H2 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 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 F2 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 N2 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 C2 curves in Figure 3E, the three VB methods (VBSCF, dc-DFVB, and λ-DFVB) and CASPT2 predict the bonding dissociation of C2 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 Cr2, 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; De) 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 deviations from the experimental data. With a large IPEA shift of 0.45 a.u., CASPT2 provides a satisfactory result for Cr2. Both BLYP and B3LYP are unable to provide the proper descriptions for C2 and Cr2 molecules, as mentioned in literature (Carlson et al., 2015; Kepp, 2017).

Table 2

The computed De for diatomic molecules (in kcal/mol).

	CASPT2	MRCI	BLYP	B3LYP	VBSCF	BOVB	dc-DFVB	λ-DFVB	Expt
H2	106.1	109.4	109.5	110.3	95.3	95.3	118.9	109.1	109.5 (Linstrom and Mallard, 2011; Johnson, 2018)
HF	133.8	135.8	139.4	138.0	113.4	124.4	139.0	142.9	141.3 (Linstrom and Mallard, 2011; Johnson, 2018)
F2	34.0	34.9	52.6	40.1	16.8	33.9	35.7	38.9	38.2 (Linstrom and Mallard, 2011; Johnson, 2018)
N2	215.6	219.9	242.3	230.2	204.1	238.6	263.2	224.3	228.5 (Linstrom and Mallard, 2011; Johnson, 2018)
C2	149.5	137.8	137.5	121.3	137.3	–	176.2	137.4	148.0 (Leininger et al., 1998)
Cr2	28.4(0.25)a 37.8(0.45)a	–	–	–	–	–	–	38.7	33.9 (Casey and Leopold, 1993)

Values in parentheses are the IPEA shifts used in CASPT2 calculations.

Because of the lack of dynamic correlation, VBSCF is unable to provide satisfactory BDEs for H2, HF, F2, and N2. Analogous to CASSCF, VBSCF is unable to describe the Cr–Cr bonding properly, but it predicts the bonding of C2 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, C2 and Cr2 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 EC functional into VBSCF, the BDEs of dc-DFVB are larger than their corresponding VBSCF values and are mostly overestimated for N2 and C2 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 N2, 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 C2, 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 Cr2, 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).

Chemical Reaction Barriers

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.

Table 3

The barriers of the D-A and Menshutkin chemical reactions (in kcal/mol).

	Menshutkin reaction	Diels–Alder reaction
	H3N + CH3Cl → H3N…CH3+…Cl− → H3NCH3+ + Cl− Reactant TS Product
CASPT2	40.5	23.5
BLYP	27.0	17.7
B3LYP	29.9	21.3 (Zhou et al., 2017)
VBSCF	41.5	41.7
dc-DFVB	38.5	34.5
λ-DFVB	32.4	24.6
Expt	33.0 (Webb and Gordon, 1999)	23.3 ± 2 (Webb and Gordon, 1999; Guner et al., 2003)

The Excitation Energy Gaps

The singlet–triplet energy gaps of carbon atom, oxygen atom, carbene (CH2), 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 CH2 and TMM, respectively, even smaller than those of CASPT2.

Table 4

The singlet–triplet energy gaps of C, O, carbene (CH2), and trimethylenemethane (TMM) (in kcal/mol).

		VBSCF	CASPT2	BLYP	B3LYP	dc-DFVB	λ-DFVB	Expt
C	3P → 1D	34.5	30.0	39.2	40.3	24.6	25.3	29.1 (Ess et al., 2011)
O	3P → 1D	50.3	46.3	60.9	62.4	40.3	41.1	45.4 (Ess et al., 2011)
CH2	3B1 → 1B1	39.2	26.9	7.3	29.8	25.9	31.3	32.9 (Ess et al., 2011)
TMM	3A2 ^′→1A1	21.5	20.1	34.7	43.9	14.7	18.4	18.1 (Li and Paldus, 2008)

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 λ-DFVB and B3LYP describe the gap correctly. The λ-DFVB result of 2.40 kcal/mol is close to the reference data of ca 3.0 kcal/mol by stochastic-CASSCF (32,34) (Manni and Alavi, 2018) and 2.0 kcal/mol by heat-bath configuration interaction with semistochastic perturbation theory at the active space (44,44) (Olivares-Amaya et al., 2015).

Figure 4

The Fe(II)–porphyrin complex.

Table 5

The triplet–quintuplet energy gap of Fe(II)–porphyrin by various methods (in kcal/mol).

Method	Energy gap
VBSCF	−26.2
dc-DFVB	−17.0
λ-DFVB	2.4
MCSCF (6,5)	−26.0
CASPT2 (6,5)	−6.7
B3LYP (Kozlowski et al., 1998)	6.2
SHCI (44,44) (Smith et al., 2017)	1.9
Stoch-CAS (32,34) (Manni and Alavi, 2018)	3.1

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 N2 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, EC(ρ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 Cr2 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.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.