A cautionary tale: Problems in the valence-CASSCF description of the ground state (X¹Σ⁺) of BF

The full optimized reaction space^1–3 (FORS) and valence complete active space self-consistent field^4–6 (vCAS) methods are powerful and popular means of describing the electronic structure of molecules that require more than a single configuration wave function. In both the FORS and vCAS methods, an active space is defined consisting of the union of the valence orbitals of the atoms, all possible configurations using the active space orbitals are constructed, and the orbitals and configuration coefficients are optimized self-consistently. FORS and vCAS wave functions account for all of the non-dynamical correlation in the atoms and molecules, but they also include an unspecified amount of dynamical correlation.⁷ An ensuing FORS+1+2 or vCAS+1+2 calculation is then used to account for the remaining dynamical correlation. Since the orbitals in the FORS+1+2 or vCAS+1+2 wave functions are taken directly from the FORS or vCAS wave function, it is critical that these wave functions and orbitals be appropriate for the FORS+1+2 and vCAS+1+2 calculations. However, this is not guaranteed.

Once the active orbitals for the FORS/vCAS calculations have been selected, the calculations can be run as “black boxes” and are often performed in this way. However, optimization of the orbitals and configuration coefficients may not lead to a consistent description of the molecule as the geometry is varied. We have, in fact, encountered such inconsistencies in a number of our studies with full vCAS and vCAS+1+2 calculations on the ground state of BF being the “poster child” for the various types of problems that can arise. Although, the closely related BH molecule is well described by the FORS/vCAS wave function and presents no problems, substitution of a fluorine atom for the hydrogen atom in BH dramatically complicates the FORS/vCAS and FORS/vCAS+1+2 calculations. In this article, we will illustrate how the FORS/vCAS descriptions of BH and BF differ and how these differences impact the subsequent FORS+1+2/vCAS+1+2 calculations. We shall also show that a spin-coupled generalized valence bond^8–10 (SCGVB) inspired wave function eliminates these inconsistencies and provides a theoretically consistent description of the electronic structure of both BH and BF.

In Sec. II, we briefly describe the SCGVB and vCAS wave functions for BH and BF. In Sec. III, we present and discuss the results of the SCGVB, vCAS, and vCAS+1+2 calculations on these two molecules. In Sec. IV, we present the conclusions drawn from the studies described in Sec. III. In the remainder of this article, we will use the notation vCAS(n_e, n_a) to represent a complete active space self-consistent field wave function for n_e electrons in n_a active orbitals. The notation vCAS(n_e, n_a)+1+2 represents a vCAS(n_e, n_a) wave function plus all single and double excitations from the active orbitals plus any doubly occupied valence orbitals.

The traditional SCGVB wave functions for the ground (X¹Σ⁺) states of BH and BF are

In Eq. (1), {ϕ_c1} refers to the B1s core orbital and {φ_a1, φ_a2, φ_a3, φ_a4} refers to the active valence orbitals—one orbital for each electron, i.e., a SCGVB(4, 4) wave function. In the separated atom limit, these orbitals become the B2s₋, B2s₊, B2p_z, and H1s atomic orbitals (see below for the definition of the B2s₋ and B2s₊ lobe orbitals). In Eq. (2), {ϕ_c1, ϕ_c2} refers to the F1s and B1s core orbitals, {φ_a1, φ_a2, φ_a3, φ_a4} refers to the active valence orbitals that become the B2s₋, B2s₊, B2p_z, and F2p_z atomic orbitals in the separated atom limit, and {ϕ_v1, ϕ_v2, ϕ_v3} refers to the doubly occupied valence orbitals that become the F2s, F2p_x, and F2p_y orbitals in the separated atom limit.

In both equations (1) and (2), the B2s lobe orbitals, which would be described by a doubly occupied orbital in the Hartree–Fock (HF) wave function, are described by a pair of non-orthogonal orbitals with the spins of the electrons in the two orbitals singlet coupled as they are in the HF wave function. These two orbitals are given by

with λ = 0.427 (the bond axis is the z-axis) and the orbitals (φ_B2s, $φB2px$⁠) are the atomic orbitals of the boron atom. The wave function for the (B2s₋ and B2s₊) pairs in Eqs. (1) and (2) can be alternately written as a linear combination of 2s² and 2p_x² configurations. Thus, this pair of orbitals accounts for the 2s–2p near-degeneracy effect in the boron atom. However, Eqs. (1) and (2) include only one of the two 2p² configurations describing the 2s–2p near-degeneracy effect—it is missing the 2p_y² configuration. Furthermore, the wave functions in Eqs. (1) and (2) do not have pure ¹Σ⁺ symmetry. Both of these problems can be addressed by using a ¹Σ⁺-projected SCGVB wave function, e.g., for the BH molecule at R = ∞, this wave function is

The SCGVB wave function in Eq. (4) describes four electrons in five orbitals—the φ_B2s, $φB2px$⁠, $φB2py$⁠, $φB2pz$⁠, and φ_H1s orbitals (note that the $φB2px$ and $φB2py$ orbitals need not be equivalent to the $φB2pz$ orbital). We refer to this wave function as a SCGVB(4,5) wave function. This wave function can be used for all R values by optimizing the orbitals and spin coupling coefficients in $Θ0,04$ in the first or second line of Eq. (4) or (c₁, c₂) in the bottom line of Eq. (4). The SCGVB(4,5) wave function for BF has the same general form as the SCGVB(4,5) wave function for BH, and the CASVB module^11–14 in Molpro^15,16 has been used to compute the SCGVB(4,5) wave functions for both molecules. In this paper, we will only report the results for the SCGVB(4,5) wave functions for the BH and BF molecules.

The active space for full vCAS calculations on the ground state (X¹Σ⁺) of BH consists of the (B2s, B2p) orbitals and the H1s orbital. In C_2v symmetry, this means that there are three orbitals of σ(a₁) symmetry and one orbital each of π_x(b₁) and π_y(b₂) symmetry in the active space, i.e., a vCAS(4,5) wave function with four electrons in five orbitals. The configurations in the vCAS(4,5) wave function include the configurations required to describe the dissociation of BH into atoms as well as all of the configurations needed to describe the 2s–2p near-degeneracy effect in the boron atom.¹⁷ There are additional configurations in the vCAS(4,5) wave function whose coefficients approach zero as R → ∞ but are otherwise non-zero. These configurations describe a small subset of the dynamical correlation effects in BH.

For the ground state (X¹Σ⁺) of BF, the full valence active space includes the (B2s, B2p) and (F2s, F2p) orbitals, resulting in an active space consisting of four orbitals of σ(a₁) symmetry and two each of π_x(b₁) and π_y(b₂) symmetry, i.e., a vCAS(10,8) wave function. This valence active space can also describe the dissociation of the BF molecule into atoms as well as the full boron atom 2s–2p near-degeneracy effect. The use of “can” in the last sentence was purposeful—the vCAS wave function will seek to describe those correlation effects that minimize the energy; for further details, see the discussion in Sec. III. Just as in BH, some of the configurations included in the vCAS(10,8) wave function can describe dynamical correlation effects.

An alternate set of active orbitals can be constructed for the BF molecule inspired by the SCGVB(4,5) wave function. In the SCGVB(4,5) wave function, the F2s, F2p_x, and F2p_y orbitals are doubly occupied and, therefore, are excluded from the active space. This set of active orbitals describes the dissociation of the BF molecule into a HF description of the fluorine atom plus the full 2s–2p near-degeneracy effect in the boron atom. This wave function is a vCAS(4,5) wave function,¹⁸ which has the same basic structure as the vCAS(4,5) wave function for BH. This wave function would be expected to primarily describe non-dynamical correlation effects, accounting for less dynamical correlation energy than is possible in the vCAS(10,8) wave function.

The vCAS+1+2 calculations included all single and double excitations from the vCAS wave function and account for dynamical correlation effects. The results of coupled cluster calculations including singles and doubles¹⁹ plus perturbative triples²⁰ [CCSD(T)] were included in the study for comparison.

All of the SCGVB(4,5), vCAS, and CCSD(T) calculations were performed with the Molpro suite of computational chemistry programs^15,16 (version 2010.1). The SCGVB(4,5) calculations used the CASVB module^11–14 in Molpro. Augmented correlation consistent basis sets of quadruple zeta quality (aVQZ) were used for the hydrogen, boron, and fluorine atoms.^21,22

The total energies (E_e), equilibrium bond lengths (R_e), dissociation energies (D_e), and fundamental frequencies (ω_e) for the SCGVB, vCAS, vCAS+1+2, and CCSD(T) calculations on the ground state (X¹Σ⁺) of BH are listed in Table I. The results of the corresponding calculations on the ground state (X¹Σ⁺) of BF are listed in Table II. Also included in the tables are the results of recent computational studies of BH^23,24 and BF^25,26 along with the available experimental data.^27–29 

The configurations in the SCGVB(4,5) wave function of BH are a significant subset of the configurations in the vCAS(4,5) wave function. Thus, it is not surprising that the SCGVB(4,5) calculations predict a bond length (R_e) that is just 0.0005 Å shorter, a dissociation energy (D_e) that is just 0.26 kcal/mol less, and a fundamental frequency just 7 cm⁻¹ more than the spectroscopic constants predicted by the vCAS(4,5) calculations; see Table I. In fact, the energies of the SCGVB(4,5) and vCAS(4,5) wave functions at their respective R_e’s differ by only −0.42 millihartrees (mE_h).

The vCAS(4,5) wave function predicts a R_e that is 0.0158 Å longer, a D_e that is 6.54 kcal/mol smaller, and an ω_e that is 75 cm⁻¹ smaller than the spectroscopic constants predicted by the vCAS(4,5)+1+2 calculations. These differences are significant and are a result of the dynamical correlation included in the vCAS(4,5)+1+2 wave function. Furthermore, the vCAS(4,5)+1+2 calculations predict spectroscopic constants that are in good agreement with the present CCSD(T) calculations—ΔR_e = 0.0003 Å, ΔD_e = −0.20 kcal/mol, and Δω_e = −2 cm⁻¹—as well as with prior calculations and experimental measurements.

The potential energy curves for BH from the SCGVB(4,5), vCAS(4,5), vCAS(4,5)+1+2, and CCSD(T) calculations are plotted in Fig. 1. The curves from the SCGVB(4,5) and vCAS(4,5) calculations are almost indistinguishable in the plot. However, there is some variation in the curves: at R = 0.7337 Å, the difference is −0.56 mE_h with the difference reaching a maximum of −0.94 mE_h at R = 2.1337 Å before decreasing to zero at R = ∞. The potential energy curves for the vCAS(4,5)+1+2 and CCSD(T) calculations are also nearly indistinguishable in Fig. 1, deviating only slightly for R > 2.7 Å where the CCSDS(T) calculations begin to break down because they are based on the RHF wave function. In summary, the SCGVB(4,5) and vCAS(4,5) provide excellent zero-order descriptions of the BH molecule, yielding semi-quantitative predictions of the spectroscopic constants, while the vCAS(4,5)+1+2 wave function provides nearly quantitative descriptions of the BH molecule for all R.

A detailed examination of the orbitals and configuration coefficients in the vCAS wave function shows that the configurations corresponding to the 2s–2p near-degeneracy effect in the boron atom are important at all R; see Fig. 2(a). If the configurations in the vCAS(4,5) wave function are classified as (n_2σn_3σn_4σ n_1πxn_1πy), where n_i is the occupation of the active valence orbital i, the behaviors of the coefficients of the configurations with occupancies of (200 2 0), (020 2 0), and (002 2 0) (and the three corresponding 1π_y configurations where n_1πy = 2) in the region R = 1.8 Å–2.4 Å reflect the fact that the B2s-like σ orbital and the BH σ bond orbital smoothly interchange in this region. This leads to an “avoided crossing” in the plot of the occupation numbers of the 2σ and 3σ orbitals; see Fig. 2(b). At large R values, 2σ is the B2s-like orbital and 3σ is the BH bonding orbital, while the opposite is true for smaller separations including R_e. For intermediate R, as the two orbitals switch character, both orbitals become delocalized over the atoms. The orbitals then re-localize after the interchange. The 4σ orbital, on the other hand, is always the antibonding orbital. Despite this orbital interchange, the SCGVB(4,5) and vCAS(4,5) wave functions, both of which include the 2σ and 3σ orbitals in a symmetrical fashion, fully account for the boron 2s–2p near-degeneracy effect at all R values. The vCAS(4,5) orbitals for BH at R_e, at an intermediate R value where the orbitals are undergoing the interchange, and at a large R value are plotted in Fig. 3.

For the BF molecule, the situation is quite different; see Table II. However, again, the spectroscopic constants from the SCGVB(4,5) and vCAS(4,5) calculations are quite similar, ΔR_e = −0.0002 Å, ΔD_e = −0.05 kcal/mol, and Δω_e = +2 cm⁻¹. In addition, ΔE_e = −0.09 mE_h at their respective R_e’s. The potential energy curves from the SCGVB(4,5) and vCAS(4,5) curves are plotted in Fig. 4. The two curves are almost indistinguishable, with the difference in the curves varying from 0.06 mE_h at R = 0.8705 Å and reaching a maximum of 1.30 mE_h at R = 1.8705 Å before decreasing to zero at R = ∞. In summary, the SCGVB(4,5) wave function clearly captures the essential features of the vCAS(4,5) wave function for BF.

Comparing the predictions from the vCAS(4,5) wave function and the vCAS(10,8) wave function, we see that the predicted R_e and ω_e are very similar, ΔR_e = +0.0016 Å and Δω_e = −21 cm⁻¹, but D_e, which is calculated as the difference between the energies at R_e and ΔR = 10 Å (ΔR = R − R_e and R_e = 1.2705 Å), is much larger for the vCAS(10,8) wave function than that for the vCAS(4,5) wave function: ΔD_e = −27.78 kcal/mol. In fact, the dissociation energy calculated with the vCAS(10,8) wave function is slightly larger than that from the present CCSD(T) calculations, a surprising result. For the vCAS+1+2 calculations, we also see reasonable agreement for (R_e, ω_e) between the vCAS(4,5)+1+2 and vCAS(10,8)+1+2 results: ΔR_e = −0.0051 Å and Δω_e = +20 cm⁻¹. However, the calculated vCAS(10,8)+1+2 D_e, 194.36 kcal/mol, is more than 10 kcal/mol larger than the experimental value.

The unusual results from the vCAS(10,8) and vCAS(10,8)+1+2 calculations on BF led us to carefully examine the details of the calculations. We found a number of problems with the vCAS(10,8) calculations, problems that also impacted the vCAS(10,8)+1+2 calculations. The first set of problems involve mixing or swapping of the core and active orbitals, while the second problem involves a change in the nature of the 2π orbitals, the orbitals that nominally arise from the B2p_x,y atomic orbitals in the vCAS(10,8) wave function.

1. At short R values, the core F1s orbital mixes with the active F2s-like orbital in the vCAS(10,8) wave function with the amount of mixing increasing as R decreases. Although this mixing lowers the energy of the vCAS(10,8) wave function unexpectedly, it increases the energy of the vCAS(10,8)+1+2 wave function, i.e., the resulting active orbital is not optimal for the vCAS(10,8)+1+2 calculations. This is the reason that the energy of the vCAS(10,8)+1+2 wave function is 3.6 mE_h above the energy of the vCAS(4,5)+1+2 wave function at their respective R_e’s; see Table II.

2. The nature of the orbitals in the vCAS(10,8) wave function depend on the direction of the scan: from R_e to ΔR = 10 Å or from ΔR = 10 Å to R_e (ΔR = R − R_e, R_e = 1.2705 Å).

   • Starting at R_e and increasing R, the 2π orbitals in the vCAS(10,8) wave function describe radial correlation of the two F2p_x,y-like pairs with the 2π orbitals resembling F3p_x,y-like orbitals. In addition, the core F1s-like orbital and active B2s-like orbital swap places at ΔR ≈ 3.9 Å. These changes again result in lower vCAS(10,8) energies unexpectedly, but at ΔR = 10 Å, the vCAS(10,8)+1+2 energy is 30.5 mE_h higher than the vCAS(4,5)+1+2 energy because of the change in the active orbitals (B2s replaced by F1s); see Table III.

   • A vCAS(10,8) calculation at ΔR = 10 Å results in 2π orbitals that describe correlation of the B2s-like pair (2s–2p near-degeneracy effect). This remains the case as R is decreased until ΔR ≈ 1.65 Å. At this point, the 2π orbitals switch to being F3p_x,y-like orbitals that describe radial correlation of the F2p_x,y-like pairs. In addition to this change, the core B1s and active F2s orbitals swap places at ΔR = 10 Å and remain this way until ΔR ≈ 2.0 Å. At ΔR = 10 Å, the energy of this vCAS(10,8)+1+2 wave function is 40.4 mE_h higher than the vCAS(10,8)+1+2 wave function that uses the 2π orbitals to correlate the F2p_x,y pairs; see Table III.

As can be seen, the problems in the vCAS(10,8) calculations on BF are wide ranging, involving both core-active orbital mixing and/or swapping. This is further complicated by changes in the basic nature of the 2π orbitals as a function of R. These problems lead to much higher vCAS(10,8)+1+2 asymptotic energies, discontinuities, and other irregularities in the vCAS(10,8) and vCAS(10,8)+1+2 potential energy curves for BF.

We explored a number of options for addressing the above problems in the vCAS(10,8) calculations. Our initial focus was on restricting the valence active space. In this regard, we sought vCAS(n_e, n_a) wave functions with n_e ≤ 10 and n_a ≤ 8 that satisfied the following two criteria:

1. The vCAS(n_e, n_a) wave function provides a theoretically consistent description of the electronic structure of BF for all R, without any mixing or swapping of the core and active orbitals.

2. The vCAS(n_e, n_a) wave function uses the 2π orbitals to correlate the B2s pair in the separated atom limit since these orbitals give the lowest vCAS+1+2 energy for this limit.

We found only one vCAS(n_e, n_a) wave function that satisfied both of these criteria—the SCGVB(4,5)-inspired vCAS(4,5) wave function. Some of the other vCAS wave functions used the 2π orbitals to correlate the F2p-like pairs at all R values. Other vCAS wave functions used the 2π orbitals to correlate the F2p-like pairs at short R values and then discontinuously switched to using a different set of 2π orbitals to correlate the B2s-like pair at large R values. Many of these wave functions also exhibited core-active orbital mixing and/or orbital swapping. In summary, only the vCAS(4,5) and vCAS(4,5)+1+2 wave functions were found to provide a theoretically consistent methodology for describing the electronic structure of both the BH and BF molecules, exhibiting none of problems associated with the vCAS(10,8) and vCAS(10,8)+1+2 wave functions for BF. The potential energy curves from the vCAS(4,5) and vCAS(4,5)+1+2 calculations for BF are shown in Fig. 4, along with the potential energy curve from the present CCSD(T) calculations.

The last option that we explored was to combine the vCAS(4,5) and vCAS(10,8) approaches to determine if the “correct” behavior of the vCAS(4,5) wave function would address the problems in the vCAS(10,8) calculations. At each value of R, we computed the vCAS(4,5) wave function and then used this wave function as the initial guess for the vCAS(10,8) calculations, freezing the B1s and F1s orbitals to be those from the vCAS(4,5) wave function. We refer to the wave functions from these calculations as the vCAS(10,8)^* and vCAS(10,8)^*+1+2 wave functions. This approach did, indeed, solve many of the above problems, giving an energy that was 4.3 mE_h below that of the vCAS(4,5)+1+2 wave function at R_e and leading to values for the spectroscopic constants for BF—R_e = 1.2654 Å, D_e = 180.16 kcal/mol, and ω_e = 1405 cm⁻¹—that are in excellent agreement with the present CCSD(T) calculations; see Table II. However, this wave function used the 2π orbitals to describe the correlation of the two F2p-like pairs at short R values, and then, in the region ΔR ≈ 1.4 Å, the 2π orbitals switched to correlating the B2s-like pair. This inconsistent description of BF as a function of R resulted in a discontinuity in the vCAS(10,8)^* potential energy curve; see Fig. 4. Fortunately, a discontinuity is not present in the curve from the vCAS(10,8)^*+1+2 calculations, although the behavior of the curve at large R values is unusual.

The valence orbitals from the vCAS(4,5) and vCAS(10,8)^* wave functions for BF at R_e and R_e + 4.0 Å (essentially the separated atom limit) are plotted in Fig. 5. Of particular interest are the 1π and 2π orbitals. For the vCAS(4,5) wave function, these orbitals correspond to the F2p_x,y-like (1π) and B2p_x,y-like (2π) orbitals at all R values. However, for the vCAS(10,8)^* wave function, these orbitals correspond to the F2p_x,y (1π) and B2p_x,y (2π) orbitals in the separated atom limit and F2p_x,y-like (1π) and F3p_x,y-like (2π) orbitals at R_e.

As a check on the performance of the vCAS and vCAS+1+2 wave functions at large R values, the last row in Table III lists the sum of the vCAS(3,4) and HF energies (−24.563 454 h, −99.409 209 h) as well as the vCAS(3,4)+1+2 and HF+1+2 energies (−24.601 200 h, −99.642 289 h) of the boron and fluorine atoms. The vCAS(4,5) and vCAS(10,8) energies for BF at ΔR = 10 Å are slightly below the sum of the atomic energies because the vCAS(3,4) calculations on the boron atom yield orbitals that are pure symmetry functions (s, p), whereas the vCAS(4,5) and vCAS(10,8) calculations yield orbitals that have only C_2v symmetry, e.g., the 2p_x,y orbitals are slightly different than the 2p_z orbitals. This additional flexibility in the vCAS BF calculations at large R values yields a slightly lower energy (−0.2 mE_h). The differences in the vCAS+1+2 calculations at ΔR = 10 Å are a result of multiple factors. The difference in the vCAS(4,5)+1+2 energy and the corresponding sum of the atomic energies is a result of the fact that the vCAS+1+2 energy is not size consistent. As a result, the sum of the atomic energies differs by −6.4 mE_h from that of the vCAS(4,5)+1+2 wave function at ΔR = 10 Å. The errors in the vCAS(10,8)+1+2 energies are much larger, −36.9 mE_h (F2p) and −77.3 mE_h (B2s), because of the interchange of core and active orbitals as well as size inconsistency in energies from these wave functions.

The problems in the full vCAS(10,8) and vCAS(10,8)+1+2 descriptions for BF are a direct result of the extraordinary variational freedom implicit in the vCAS wave function. Optimization of both the orbitals and the configuration coefficients in the full configuration space associated with the active orbitals enables the vCAS wave function to describe non-dynamical and dynamical correlation effects in a potentially unpredictable way. Although the full vCAS wave functions for both BH and BF are able to describe the 2s–2p near-degeneracy effect in the boron atom, only for BH is this found to be the case at all internuclear separations. In BF, radial correlation of the electrons in the doubly occupied F2p-like orbitals, a dynamical correlation effect, is energetically more favorable at intermediate to short internuclear distances, and the 2π orbitals are used to describe this effect, not the B2s–2p near-degeneracy effect. In addition, because of orbital mixing and core-active orbital swapping, the active orbitals from the full vCAS calculations are often not optimal for the vCAS+1+2 calculations, leading to an unreasonably large D_e’s.

The scope of the vCAS and vCAS+1+2 calculations could, of course, be increased by including additional orbitals in the active space. Although this extension is beyond the scope of the current work, the above results suggest that it will still be critical to monitor the nature of the orbitals as a function of R to ensure that a consistent description of the BF molecule is obtained for all values of R.

Although it is tempting to use the complete active space self-consistent field (CAS) method as a “black box”—select a set of active orbitals and run the calculation—the full valence-CAS (vCAS) calculations on BF presented here serve as a cautionary tale. Selecting the active orbitals as the union of a set of atomic orbitals does not imply that these orbitals maintain their basic identities for the calculations. In BH, the 1π orbitals kept their basic identities as boron 2p-like orbitals in the full vCAS calculations, describing the 2s–2p near-degeneracy effect in the boron atom. In the full vCAS wave function of BF, on the other hand, the corresponding 2π orbitals are F3p-like orbitals because the correlation of the two fluorine 2pπ pairs is energetically more important than the correlation of the B2s pair (2s–2p near-degeneracy effect) at R_e. The full vCAS calculations are further complicated by the fact that there is mixing between the core F1s and active F2s orbitals at short R, swapping of the core B1s and active F2s orbitals and between the core F1s and active B2s orbitals at large R, with changes in the nature of the orbitals depending on how R is scanned (from R_e → ∞ or ∞ → R_e). Although these changes lower the full vCAS energy, they do not necessarily lead to the lowest vCAS+1+2 energy.

Although a larger active space will always provide a better vCAS energy than a smaller active space of which it is a superset, it does not follow that the vCAS+1+2 energies will do the same. This is the case for BF where the vCAS(4,5) wave function gave a better vCAS(4,5)+1+2 energy than the vCAS(10,8)+1+2 calculation. This rather unsettling finding is a consequence of the unusual behavior of the vCAS(10,8) wave function where core-active orbital mixing and swapping cause the active orbitals to be suboptimal for the vCAS(10,8)+1+2 wave function. We found only one computationally sound means of addressing the problems in the vCAS(10,8) and vCAS(10,8)+1+2 calculations on BF—the use of the SCGVB(4,5)-inspired vCAS(4,5) wave function for the vCAS+1+2 calculations. This approach provides very good predictions of (R_e, D_e, and ω_e) in comparison to the results from the present CCSD(T) calculations: ΔR_e = −0.0034 Å, ΔD_e = −3.32 kcal/mol, and Δω_e = +15 cm⁻¹. Furthermore, the potential energy curve exhibits no discontinuities or other troublesome features. In comparison, the differences in these same quantities for the BH molecule are ΔR_e = 0.0003 Å, ΔD_e = −0.20 kcal/mol, and Δω_e = −2 cm⁻¹.

We found that the use of the F1s and B1s orbitals from the vCAS(4,5) wave function in the vCAS(10,8) calculations, yielding the CAS(10,8)^* wave function, eliminates the problems with the mixing and swapping of the core and active orbitals. Furthermore, the spectroscopic constants predicted by vCAS(10,8)^*+1+2 calculations are in excellent agreement with those from the present CCSD(T) calculations, ΔR_e = −0.0034 Å, ΔD_e = −0.61 kcal/mol, and Δω_e = 14 cm⁻¹. Unfortunately, at large R, there was a discontinuous change in the nature of the 2π orbitals in the vCAS(10,8)^* calculations, from F3p_x,y-like orbitals at R ≲ 2.7 Å to B2p_x,y-like orbitals at R ≳ 2.7 Å, resulting in a discontinuity on the vCAS(10,8)^* potential energy curve. Although the curve from the corresponding vCAS(10.8)^*+1+2 calculations did not exhibit such a discontinuity, it still exhibited strange behavior as it approached the asymptote at large R.

We have found similar problems with full vCAS calculations on other molecules, including molecules that contain other atoms in the latter part of the first row. Although the multiconfiguration vCAS method is a powerful means of describing the electronic structure of molecules that are not well described by a single configuration method, the tremendous flexibility of the vCAS method is a mixed blessing. As shown here, it is critical to pay close attention to both the orbitals and configuration coefficients in the vCAS calculations to ensure that the calculations provide a theoretically consistent description of the electronic structure of molecules for all geometries of interest. Although this check was relatively easy to perform for BH and BF, it will be far more challenging, yet undoubtedly just as important, for vCAS calculations on larger molecules. The problems with the full vCAS wave function of BF can be avoided by using the SCGVB(4,5)-inspired vCAS(4,5) and vCAS(4,5)+1+2 wave functions, and a similar approach could be used for other molecules.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This research was supported, in part, by the Center for Scalable Predictive Methods for Excitations and Correlated phenomena (SPEC), which is funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, as part of the Computational Chemical Sciences Program at Pacific Northwest National Laboratory. We wish to thank Professor Kirk A. Peterson for his reading of and comments on this manuscript and Professor David L. Cooper for his advice on the use of the CASVB module in Molpro for the SCGVB(4,5) calculations.