Front. Chem.Frontiers in ChemistryFront. Chem.2296-2646Frontiers Media S.A.10.3389/fchem.2019.00225ChemistryOriginal Researchλ-Density Functional Valence Bond: A Valence Bond-Based Multiconfigurational Density Functional Theory With a Single Variable Hybrid ParameterYingFuming123ZhouChen123ZhengPeikun123LuanJiamin123SuPeifeng123*WuWei1231Fujian Provincial Key Laboratory of Theoretical and Computational Chemistry, Xiamen University, Xiamen, China2The State Key Laboratory of Physical Chemistry of Solid Surfaces, Xiamen University, Xiamen, China3College of Chemistry and Chemical Engineering, Xiamen University, Xiamen, China
Edited by: Stephane Humbel, Aix-Marseille Université, France
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
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.
valence bond (VB) methodmulti-configurationdensity functional theorymulti-reference characterstrong correlation215031722157317621733008National Natural Science Foundation of China10.13039/501100001809Introduction
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),
Ψ=∑KCKΦK,
where ΦK and CK 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
ΦK=ÂΩ0ΘK,
where  is an antisymmetrizer for electron indices, Ω0 is an orbital product,
Ω0=ϕ1(1)ϕ2(2)⋯ϕN(N),
and ΘK 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 kp is for unpaired electrons.
Coefficients {CK} in Equation 1 can be obtained by solving the secular equation:
HC=EMC,
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 {CK} and VB orbitals {φi} are optimized simultaneously to minimize the total energy E. VB orbitals are usually expanded as linear combinations of basis functions.
ϕi=∑μTμiχμ.
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,
WK=∑LCKMKLCL.
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
Edc-DFVB[ρ]=minΨ{〈Ψ|T+Vext+Wee|Ψ〉+EC[ρ]},
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:
EMC1H=minΨ{〈Ψ|T+Vext+λWee|Ψ〉+EHXCλ[ρ]},
where T, Vext, and Wee are the kinetic energy, external potential, and electron–electron interaction operators, respectively; EHXCλ[ρ]is the complement λ-dependent Hartree–exchange–correlation density functional for electronic density ρ, and λ is a coupling parameter. EHXCλ[ρ]is defined as (Sharkas et al., 2012):
EHXCλ[ρ]=(1-λ)(EH[ρ]+EX[ρ])+(1-λ2)EC[ρ],
where EX[ρ] and EC[ρ] are the exchange and correlation functionals, respectively. EH is the Hartree energy, given as
EH=12∬ρ(r)ρ(r′)|r-r′|drdr′.
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-C02 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, FA, is defined as
FA=VA-∑B,B≠AOAB,
where VA and OAB are the total valence of atom A and the bond order between atoms A and B, respectively, defined as (Mayer, 2003)
VA=∑μ∈A2(DS)μμ-∑μ,ν∈A(DS)μν(DS)νμ,
and
OAB=∑μ∈A∑ν∈B[(DS)μν(DS)νμ+(PsS)μν(PsS)νμ],
In Equations 13 and 14, S is the overlap matrix in terms of basis functions, D = Pα + Pβ is the total density matrix, and Ps = 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
K=∑AFA∑AVA.
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 OAB = 0 and FA = VA.
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
λ=K1/4.
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
When a molecule composed of two fragments/atoms A and B is fully dissociated, λ = 1.0. Then, the total energy can be expressed as
EABλ-DFVB=EABVBSCF+EC[ρLD].
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 EHXCλ-DFVB[ρ] by Equation 18.
Compute the operator vHXCλ-DFVB[ρ], which is defined as
vHXCλ-DFVB[ρ]=δEHXCλ-DFVB[ρ]δρ.
Optimize the λ-DFVB wave function with the following equation:
ελ-DFVB=〈Ψ|T+Vext+λWee+∫drvHXCλ-DFVB[ρ]n(r)|Ψ〉.
where n(r) is the density operator, ρ(r) = 〈Ψ| n(r) |Ψ〉. The contribution of vHXCλ-DFVB[ρ] 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:
Eλ-DFVB=ελ-DFVB+EHXCλ-DFVB[ρ]-〈Ψ|∫drvHXCλ-DFVB[ρ]n(r)|Ψ〉.
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 DiscussionsThe λ 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.
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.
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
Ecorr = Eλ−DFVB-EVBSCF.
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.
The curves of λ as functions of bond distances for diatomic molecules.
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).
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.
The barriers of the D-A and Menshutkin chemical reactions (in kcal/mol).
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.
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).
The Fe(II)–porphyrin complex.
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.
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.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2019.00225/full#supplementary-material
ReferencesAnderssonK.MalmqvistP. Å.RoosB. O. (1992). Second-order perturbation theory with a complete active space self-consistent field reference function. 96, 1218–1226. 10.1063/1.462209AndraeD.HäußermannU.DolgM.StollH.PreußH. (1990). Energy-adjusted ab initio pseudopotentials for the second and third row transition elements. 77, 123–141. 10.1007/BF01114537AquilanteF.AutschbachJ.CarlsonR. K.ChibotaruL. F.DelceyM. G.De VicoL.. (2016). Molcas 8: New capabilities for multiconfigurational quantum chemical calculations across the periodic table. 37, 506–541. 10.1002/jcc.2422126561362Benavides-RiverosC. L.LathiotakisN. N.MarquesM. A. L. (2017). Towards a formal definition of static and dynamic electronic correlations. 19, 12655–12664. 10.1039/C7CP01137G28474027BobrowiczF. W.GoddardW. A.III. (1977). “The self-consistent field equations for generalized valence bond and open-shell Hartree-Fock wave functions,” in , ed H. F. Schaefer III (New York, NY: Plenum), 79–127.BryndaM.GagliardiL.RoosB. O. (2009). Analysing the chromium–chromium multiple bonds using multiconfigurational quantum chemistry. 471, 1–10. 10.1016/j.cplett.2009.02.006CarlsonR.TruhlarD.GagliardiL. (2015). Multiconfiguration pair-density functional theory: a fully translated gradient approximation and its performance for transition metal dimers and the spectroscopy of Re2Cl82-. 11, 4077–4085. 10.1021/acs.jctc.5b00609CaseyS. M.LeopoldD. G. (1993). Negative ion photoelectron spectroscopy of chromium dimer. 97, 816–830. 10.1021/j100106a005CembranA.SongL.MoY.GaoJ. (2009). Block-localized density functional theory (BLDFT), diabatic coupling, and their use in valence bond theory for representing reactive potential energy surfaces. 5, 2702–2716. 10.1021/ct900289820228960ChenZ.SongJ.ShaikS.HibertyP. C.WuW. (2009). Valence bond perturbation theory. 113, 11560–11569. 10.1021/jp903011j19569658ChenZ.YingF.ChenX.SongJ.SuP.SongL.. (2015). XMVB 2.0: A new version of Xiamen valence bond program. 115, 731–737. 10.1002/qua.24855ChirgwinB. H.CoulsonC. A. (1950). The electronic structure of conjugated systems. 201, 196–209. 10.1098/rspa.1950.0053CooperD. L.GerrattJ.RaimondiM. (1991). Applications of spin-coupled valence bond theory. 91, 929–964. 10.1021/cr00005a014EssD. H.JohnsonE. R.HuX.YangW. (2011). Singlet–triplet energy gaps for diradicals from fractional-spin density-functional theory. 115, 76–83. 10.1021/jp109280y21141988FilatovM.ShaikS. (1999). Application of spin-restricted open-shell Kohn-Sham method to atomic and molecular multiplet states. 110, 116–125. 10.1063/1.477941FromagerE.RéalF.WåhlinP.WahlgrenU.JensenH. J. A. (2009). On the universality of the long-/short-range separation in multiconfigurational density-functional theory. II. Investigating f0 actinide species. 131:054107. 10.1063/1.318703219673551FromagerE.ToulouseJ.JensenH. J. A. (2007). On the universality of the long-/short-range separation in multiconfigurational density-functional theory. 126:074111. 10.1063/1.256645917328597GordonM. S.SchmidtM. W. (2005). “Advances in electronic structure theory: GAMESS a decade later,” in , eds. C. E. Dykstra, G. Frenking, K. S. Kim, and G. E. Scuseria (Amsterdam: Elsevier), 1167–1189. 10.1016/B978-044451719-7/50084-6GräfensteinJ.CremerD. (2000). Can density functional theory describe multi-reference systems? . 2, 2091–2103. 10.1039/a909905kGrafensteinJ.CremerD. (2005). Development of a CAS-DFT method covering non-dynamical and dynamical electron correlation in a balanced way. 103, 279–308. 10.1080/00268970512331318858GunerV.KhuongK. S.LeachA. G.LeeP. S.BartbergerM. D.HoukK. (2003). A standard set of pericyclic reactions of hydrocarbons for the benchmarking of computational methods: the performance of ab initio, density functional, CASSCF, CASPT2, and CBS-QB3 methods for the prediction of activation barriers, reaction energetics, and transition state geometries. 107, 11445–11459. 10.1021/jp035501wGusarovS.MalmqvistP.-Å.LindhR.RoosB. O. (2004). Correlation potentials for a multiconfigurational-based density functional theory with exact exchange. 112, 84–94. 10.1007/s00214-004-0568-1Head-GordonM. (2003). Characterizing unpaired electrons from the one-particle density matrix. 372, 508–511. 10.1016/S0009-2614(03)00422-6HibertyP. C.FlamentJ. P.NoizetE. (1992). Compact and accurate valence bond functions with different orbitals for different configurations: application to the two-configuration description of F2. 189, 259–265. 10.1016/0009-2614(92)85136-XHibertyP. C.HumbelS.ByrmanC. P.Van LentheJ. H. (1994). Compact valence bond functions with breathing orbitals: application to the bond dissociation energies of F2 and FH. 101, 5969–5976. 10.1063/1.468459HibertyP. C.ShaikS. (2002). Breathing-orbital valence bond method - a modern valence bond method that includes dynamic correlation. 108, 255–272. 10.1007/s00214-002-0364-8HibertyP. C.ShaikS. (2008). . Hoboken, NJ: John Wiley.HohenbergP.KohnW. (1964). Inhomogeneous electron gas. 136, B864–B871. 10.1103/PhysRev.136.B864HuangJ.YingF.SuP.WuW. (2014). VBEFP/PCM: a QM/MM/PCM approach for valence-bond method and its application for the vertical excitations of formaldehyde and acetone in aqueous solution. 57, 1409–1417. 10.1007/s11426-014-5192-xHuangZ.WangH.KaisS. (2006). Entanglement and electron correlation in quantum chemistry calculations. 53, 2543–2558. 10.1080/09500340600955674JanssenC. L.NielsenI. M. B. (1998). New diagnostics for coupled-cluster and Møller–Plesset perturbation theory. 290, 423–430. 10.1016/S0009-2614(98)00504-1Johnson IIIR. D. (2018). . Available online at: http://cccbdb.nist.gov/. (accessed April 4, 2019).KeppK. P. (2017). Trends in strong chemical bonding in C2, CN, CN-, CO, N2, NO, NO+, and O2. 121, 9092–9098. 10.1021/acs.jpca.7b0820129112409KochW.HolthausenM. C. (2001). . Weinheim: Wiley-VCH.KohnW.ShamL. J. (1965). Self-consistent equations including exchange and correlation effects. 140, A1133–A1138. 10.1103/PhysRev.140.A1133KozlowskiP. M.SpiroT. G.BércesA.ZgierskiM. Z. (1998). Low-lying spin states of iron(II) porphine. 102, 2603–2608. 10.1021/jp973346dKurzweilY.LawlerK. V.Head-GordonM. (2009). Analysis of multi-configuration density functional theory methods: theory and model application to bond-breaking. 107, 2103–2110. 10.1080/00268970903160597LeeT. J.TaylorP. R. (1989). A diagnostic for determining the quality of single-reference electron correlation methods. 36, 199–207. 10.1002/qua.560360824LeiningerM. L.NielsenI. M. B.CrawfordT. D.JanssenC. L. (2000). A new diagnostic for open-shell coupled-cluster theory. 328, 431–436. 10.1016/S0009-2614(00)00966-0LeiningerM. L.SherrillC. D.AllenW. D.IiiH. F. S. (1998). Benchmark configuration interaction spectroscopic constants for X1Σg+ C2 and X 1 Σ + CN+. 108, 6717–6721. 10.1063/1.476087LiX.PaldusJ. (2008). Electronic structure of organic diradicals: Evaluation of the performance of coupled-cluster methods. 129, 174101. 10.1063/1.299956019045327LieG. C.ClementiE. (1974). Study of the electronic structure of molecules. . 60, 1288–1296. 10.1063/1.1681193LinstromP. J.MallardW. G. (2011). . Gaithersburg, MD: National Institute of Standards and Technology.ManniG.CarlsonR.LuoS.MaD.OlsenJ.TruhlarD.. (2014). Multiconfiguration pair-density functional theory. 10, 3669–3680. 10.1021/ct500483tManniG. L.AlaviA. (2018). Understanding the mechanism stabilizing intermediate spin states in Fe(II)-porphyrin. 122, 4935–4947. 10.1021/acs.jpca.7b12710MayerI. (2003). . Springer: New York.MiehlichB.StollH. S. (1997). A correlation-energy density functional for multideterminantal wavefunctions. 91, 527–536. 10.1080/002689797171418MoY.LinZ.WuW.ZhangQ. (1996). Bond-distorted orbitals and effects of hybridization and resonance on C-C bond lengths. 100, 11569–11572. 10.1021/jp953433aMoY.WuW.ZhangQ. (1994). The treatment of big basis sets in valence bond method. 15, 899–902.NakanoH. (1993). Quasidegenerate perturbation theory with multiconfigurational self-consistent-field reference functions. 99, 7983–7992. 10.1063/1.465674Olivares-AmayaR.HuW.NakataniN.SharmaS.YangJ.ChanG. K.-L. (2015). The ab-initio density matrix renormalization group in practice. 142, 034102–034102. 10.1063/1.490532925612684ParrR. G.YangW. (1989). . New York, NY: Oxford University Press. 16592663PaunczR. (1979). . New York, NY: Plenum Press.Pérez-JiménezÁ. J.Pérez-JordáJ. M.IllasF. (2004). Density functional theory with alternative spin densities: Application to magnetic systems with localized spins. 120, 18–25. 10.1063/1.163002115267256Pou-AmérigoR.MerchánM.Nebot-GilI.WidmarkP.-O.RoosB. O. (1995). Density matrix averaged atomic natural orbital (ANO) basis sets for correlated molecular wave functions. 92, 149–181. 10.1007/BF01114922PradhanA. D.PatridgeH.BauschlierC. (1994). The dissociation energy of CN and C2. 101, 3857–3861. 10.1063/1.467503Ramos-CordobaE.MatitoE. (2017). Local descriptors of dynamic and nondynamic correlation. 13, 2705–2711. 10.1021/acs.jctc.7b0029328520420Ramos-CordobaE.SalvadorP.MatitoE. (2016). Separation of dynamic and nondynamic correlation. 18, 24015–24023. 10.1039/C6CP03072F27523386RapacioliM.SpiegelmanF.ScemamaA.MirtschinkA. (2010). Modeling charge resonance in cationic molecular clusters: combining DFT-tight binding with configuration interaction. 7, 44–55. 10.1021/ct100412f26606217Rodriguez-MayorgaM.Ramos-CordobaE.Via-NadalM.PirisM.MatitoE. (2017). Comprehensive benchmarking of density matrix functional approximations. 19, 24029–24041. 10.1039/C7CP03349D28832052RoosB. O.LindhR.MalmqvistP. A.VeryazovV.WidmarkP. O. (2005). New relativistic ANO basis sets for transition metal atoms. 109, 6575–6579. 10.1021/jp058112616834004RoosB. O.TaylorP. R.SiegbahnP. E. M. (1980). A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach. 48, 157–173. 10.1016/0301-0104(80)80045-0RuipierezF.AquilanteF.UgaldeJ.InfanteI. (2011). Complete vs restricted active space perturbation theory calculation of the Cr2 potential energy surface. 7, 1640–1646. 10.1021/ct200048zSchmidtM. W.BaldridgeK. K.BoatzJ. A.ElbertS. T.GordonM. S.JensenJ. H.. (1993). General atomic and molecular electronic structure system. 14, 1347–1363. 10.1002/jcc.540141112SearsJ. S.SherrillC. D. (2008a). Assessing the performance of density functional theory for the electronic structure of metal–salens: The d2-metals. 112, 6741–6752. 10.1021/jp802249n18593130SearsJ. S.SherrillC. D. (2008b). Assessing the performance of density functional theory for the electronic structure of metal–salens: The 3d0-metals. 112, 3466–3477. 10.1021/jp711595wSharkasK.SavinA.JensenH. J.ToulouseJ. (2012). A multiconfigurational hybrid density-functional theory. 137:044104. 10.1063/1.473367222852594SiegbahnP. E. M.AlmlöfJ.HeibergA.RoosB. O. (1981). The complete active space SCF (CASSCF) method in a Newton–Raphson formulation with application to the HNO molecule. 74, 2384–2396. 10.1063/1.441359SiegbahnP. E. M.HeibergA.RoosB. O.LevyB. (1980). A comparison of the super-CI and the Newton-Raphson scheme in the complete active space SCF method. 21:323. 10.1088/0031-8949/21/3-4/014SmithJ. E. T.MussardB.HolmesA. A.SharmaS. (2017). Cheap and near exact CASSCF with large active spaces. 13, 5468–5478. 10.1021/acs.jctc.7b0090028968097SongL.WuW.ZhangQ.ShaikS. (2004). A practical valence bond method: a configuration interaction method approach with perturbation theoretic facility. 25, 472–478. 10.1002/jcc.1038214735567SuP.WuW. (2013). Ab initio nonorthogonal valence bond methods. 3, 56–68. 10.1002/wcms.1105TishchenkoO.ZhengJ.TruhlarD. G. (2008). Multireference model chemistries for thermochemical kinetics. 4, 1208–1219. 10.1021/ct800077r26631697van LentheJ. H.Balint-KurtiG. G. (1980). The valence-bond scf (VB SCF) method: synopsis of theory and test calculation of oh potential energy curve. 76, 138–142. 10.1016/0009-2614(80)80623-3van LentheJ. H.Balint-KurtiG. G. (1983). The valence-bond self-consistent field method (VB-SCF): theory and test calculations. 78, 5699–5713. 10.1063/1.445451WebbS. P.GordonM. S. (1999). Solvation of the Menshutkin reaction: a rigorous test of the effective fragment method. 103, 1265–1273. 10.1021/jp983781nWuQ.ChengC. L.Van VoorhisT. (2007). Configuration interaction based on constrained density functional theory: a multireference method. 127:164119. 10.1063/1.280002217979331WuW.SongL.CaoZ.ZhangQ.ShaikS. (2002). Valence bond configuration interaction: a practical ab initio valence bond method that incorporates dynamic correlation. 106, 2721–2726. 10.1021/jp0141272WuW.SuP.ShaikS.HibertyP. C. (2011). Classical valence bond approach by modern methods. 111, 7557–7593. 10.1021/cr100228r21848344WuW.WuA.MoY.LinM.ZhangQ. (1998). Efficient algorithm for the spin-free valence bond theory. I. New strategy and primary expressions. . 67, 287–297. 10.1002/(SICI)1097-461X(1998)67:5<287::AID-QUA2>3.0.CO;2-RYamanakaS.NakataK.TakadaT.KusakabeK.UgaldeJ.YamaguchiK. (2006). Recent development of multireference density functional theory. 35, 242–247. 10.1246/cl.2006.242YingF.SuP.ChenZ.ShaikS.WuW. (2012). DFVB: A density-functional-based valence bond method. 8, 1608–1615. 10.1021/ct200803h26593654ZanardiP. (2002). Quantum entanglement in fermionic lattices. 65:042101. 10.1103/PhysRevA.65.042101ZeischeP.GunnarsonO.JohnW.BeckH. (1997). Two-site Hubbard model, the Bardeen-Cooper-Schrieffer model, and the concept of correlation entropy. 55, 10270–10277. 10.1103/PhysRevB.55.10270ZhouC.ZhangY.GongX.YingF.SuP.WuW. (2017). Hamiltonian matrix correction based density functional valence bond method. 13, 627–634. 10.1021/acs.jctc.6b0114427992721