Binding and excitations in Si_xH_y molecular systems using quantum Monte Carlo

Quantum Monte Carlo (QMC) methods belong to a rapidly developing family of many-body approaches that address the central challenge of accurate electronic structure calculations posed by electron–electron correlations. Several flavors of QMC with different sampling strategies are being developed to achieve a better accuracy and to be used in a wider variety of applications, including model Hamiltonians as well as real molecules and solids. Beyond perhaps the most established QMC in real space of particle coordinates,^1–3 there are several successful alternative approaches (for example, auxiliary-field QMC, full configuration interaction QMC, just to name a few^4,5). In this work, we use the diffusion Monte Carlo (DMC) algorithm that relies on the fixed-node (FN) approximation to deal with the fermion sign inefficiencies.⁶ We simplify the problem by using a valence-only Hamiltonian based on recently derived correlation consistent effective core potentials (ccECPs). Our ccECPs provide a new level of accuracy and fidelity to the original all-electron Hamiltonian as demonstrated in both the atomic and molecular systems.^7–11 

One of our goals is to test and benchmark these recently established ccECPs^8–11 in realistic systems beyond just the training problems. This part of work also involves cross-checks with existing all-electron calculations that harness state-of-the-art alternative many-body wave function methods. The second focal point is to better understand the fixed-node errors in simple molecular systems where high accuracy can be verified by alternative approaches and combined into reference quality data. The third goal is to understand the impact of fixed-node errors on excited states, which is an area where much less is known overall. Indeed, more data on the applicability of QMC methods to excited states of various types (different spin channels, adiabatic vs vertical, Rydberg-like, etc) would be fruitful so as to expand upon a few previous studies.^12,13

Considering these aims, we have opted for a few molecular Si_xH_y systems where extensive basis sets and alternative many-body wave function expansion methods can be used to the full extent. The reason for choosing the particular element of Si was motivated by the fact that its valence electronic structure exhibits some of the smallest fixed-node errors in QMC calculations of the atom, molecule, and also solid systems.¹⁴ Surprisingly, our previous results suggested that this is true even though only single-reference trial wave functions have been employed.¹⁴ Therefore, another interesting issue we wanted to elucidate was the degree of accuracy that can be achieved for excitations using the same single-reference trial function setting. We have included several types of excitations such as singlets vs triplets, as well as low-lying states in the range of an eV, and also very high excitations (≈10 eV) with Rydberg-like character. The idea behind studying excitations was motivated by an effort to build a dataset that can be used for assessment of fixed-node errors in larger systems where expansions in excited states within either Configuration Interaction (CI) or Coupled Cluster (CC) methods are very limited or not feasible.

We employ our recently generated ccECPs for Si⁹ and H.¹⁰ Previously, extensive ccECP transferability tests on atomic spectra and on SiO, Si₂, H₂ molecular binding energies, equilibrium geometries, and vibrational frequencies have been verified to provide very high accuracy, in some cases significantly beyond widely used ECPs of previous generations.¹⁵ Our ccECPs and the corresponding basis sets are available in Ref. 7.

Several quantum chemical methods such as coupled cluster with double excitations and perturbative triples [CCSD(T)], as well as triples with perturbative quadruples [CCSDT(Q)], configuration interaction with double excitations (CISD), and also triple excitations (CISDT), and finally, configuration interaction using a perturbative selection made iteratively (CIPSI) have been used. We use aug-cc-pVnZ basis sets throughout this work since we observed that the augmentations have been important in the description of molecular bonds and some of the excitations. The basis sets are TZ − 6Z quality for CCSD(T) and DZ − 5Z quality for CI/CIPSI calculations. Using these sets enabled us to recover the total energies and differences within about half of the chemical accuracy threshold (1 kcal/mol ≈ 0.0434 eV). Whenever possible, Hartree–Fock (HF) and correlation energies are extrapolated to the complete basis set (CBS) limit as¹⁶ 

where n is the basis set cardinal number. CI/CIPSI energies with multiple roots were also extrapolated with Eq. (1). In some cases, CIPSI energies were extrapolated also in regard to PT2 corrections (using the last two data points), shown as exFCI. CCSDT(Q) calculations were only feasible for small basis sets such as DZ or TZ; therefore, we estimated the CCSDT(Q) CBS limit value by evaluating the energy decrease from CCSD(T) to CCSDT(Q) at a specific basis set and adding that to the CCSD(T) CBS value,

This estimation was observed to be reasonable on the Si pseudo-atom where CCSDT(Q) was carried out for DZ-6Z basis sets.¹⁷ Additionally, Ref. 17 demonstrates that CCSDT(Q) energies are very close to FCI values (CBS values are the same within the fitting errors).

On the QMC side, we have carried out FN-DMC calculations for each system using single-reference and in some cases CI wave functions with an increasing number of determinants (N_det). Unless otherwise specified, we used orbitals from density functional theory (DFT) that are obtained at an augmented QZ basis set level. We have tested how the ground and excited states depend on the basis set for DZ-6Z levels, and we have found that both states were essentially converged at the QZ basis set level. We have also tried a doubly augmented basis set for SiH₄, and we did not observe noticeable DMC energy gain for the ground state and for the Rydberg-like singlet excitation. In all DMC calculations, we use time step τ = 0.001 hartree⁻¹ with a T-move algorithm,^18,19 which makes the DMC energies rigorously variational.

For each system, we evaluated the ground state and the first singlet and triplet excited state energies. Ground and triplet states can be expressed with a single determinant wavefunction making it suitable for calculation in CC, CI, and QMC methods. However, for singlet excited states, the correct space-spin symmetry typically leads to a multi-determinant description that is generally challenging for CC methods available in quantum chemistry packages. Therefore, for these cases, we have not attempted to calculate the single-reference CC data. In QMC, it is straightforward to build such wave functions, in the simplest cases as a linear combination of two spin-up × spin-down determinantal products with a single promoted orbital in one or the other spin channel. Clearly, it is also straightforward to employ CIS, CISD, or higher-order excitation wave functions in QMC.

In all cases, we have probed for relaxation effects in excited state orbitals when constructing the mentioned two-determinant wave functions. To do so, we employed PBE/PBE0 functionals in GS SCF calculation or even used orbitals from different spin symmetry. The most optimal orbitals obtained are indicated where necessary, and we also include corresponding HF reference DMC values for comparison. For the fixed-node DMC and other many-body calculations that use a given set of self-consistent orbitals, we employ notation “Method/Orbitals” throughout this paper.

Note that the DMC method can access the excited states of the given Hamiltonian because the constructed nodal surface constrains it to converge to the ground state of that symmetry. This setting corresponds to all the excitations considered here so that the DMC method provides variational results¹.

All molecular geometries in this work are obtained from experimental values.^20–22 Table I provides these geometries. Software packages MOLPRO²³ and MRCC²⁴ were used for CC calculations. PySCF,²⁵ GAMESS,²⁶ and QUANTUM packages²⁷ were employed to generate DFT orbitals and CI expansions. Subsequent DMC calculations were performed using QWalk²⁸ and QMCPACK^29,30 packages.

In this subsection, we provide the ground state (GS) and excited state (EX) total energies and excitation gaps using various methods as specified above. The results for the Si atom are taken from Ref. 17. Table II shows the total energies as well as the singlet and triplet gaps of the SiH₂ molecule where we opted to calculate both the vertical and adiabatic excitations since there is a significant geometry relaxation in the adiabatic case according to the experiments.²⁰ The total energies and vertical gaps for SiH₄ and Si₂H₆ are listed in Tables III and IV, respectively.

Our calculations show solid consistency throughout various methods. For each of the molecules, single-reference DMC calculations provide remarkably good accuracy for both the total energies and excitation gaps of singlet and triplet states. In SiH₂ and Si₂H₆, single-reference DMC excitation gaps also agree very well with experimental values that are known reliably. We list CCSD(T) and CCSDT(Q) calculations to illustrate the exact or nearly exact energies obtained by the CBS extrapolations. General agreement between the values from DMC/DFT and CC extrapolations is in the range of or less than ≈0.1 eV. For the singlet excitations, DMC/DFT gaps are compared with the CI method with the caveat of less accurate total energies but reasonably good estimations of the gaps. In SiH₂, we also provide CIPSI+PT2/NatOrb CBS energies with a large number of determinants showing that single-reference DMC gaps are accurate to almost the chemical accuracy level. In SiH₄, we observe higher DMC/PBE excitation gap discrepancies with some experiments for the vertical states. Inspection of the experiments³⁷ shows that the singlet peak is very broad, and considering the high excitation energies in optical absorption, we conjecture that the measured spectrum might have involved contributions from the adiabatic singlet (which is roughly 0.5 eV lower). Here, we believe that the genuine vertical excitation is close to 9.4 eV.

Similar considerations apply to the disilane triplet excitation where we conjecture that the measured value probably involves some degree of adiabacity. Hence, we believe that our values for vertical triplet clustered around 6.7 eV represent a valid prediction for this state. This requires a further study including cleaner and more specific experiments. Considering the overall accuracy of our calculations and consistency of diversified methods, we assume that our result is closer to the true vertical excitation energy.

Table V presents the Si_xH_y computational and experimental atomization energies. Included are values from extensive CCSD(T) calculations with large basis sets up to 6Z and limited CCSDT(Q) calculations to probe for convergence due to higher-order excitations. This enabled us to check the extrapolation results more carefully and provided what we believe are very accurate, nearly exact estimations of the corresponding energy differences.

Included also are results from published high-accuracy studies^42–44 that we use for comparison and also as a cross-check with all-electron results. We find very close agreements with the results of these independent sets of calculations. Indeed, the differences are at the level of about ≈0.07 eV or smaller, thus showing a very clear consistency between different methods. In fact, the differences between these calculations appear to be overall smaller than are the differences with experiments; see Table V. Our high-accuracy calculations together with results quoted from references, therefore, suggest that experiments mildly underestimate the atomization energies. This appears to be true also after all the corrections for core-valence, zero-point motion, and relativity are taken into account as outlined in the references.

We believe that our results provide a strong argument in favor of the accuracy of the employed ccECPs. Below, we further analyze the excitations; however, we expect them to be less sensitive to ccECP’s quality due to their more delocalized nature and decreased electron density in the core regions. Therefore, the main source of errors in excitations is expected to be due to the fixed-node bias.

In order to simplify the analysis of nodal biases, we probed for the impact of basis set sizes using silane as an example (Table VI). We found essentially monotonous improvements in total energy with basis set size, although beyond the quadruple zeta basis the improvements have become marginal at a 0.1 mHa level. Quite surprisingly, a better basis does not automatically imply better nodes (see counter-examples of Bressanini et al.⁴⁷). Note that the kinetic energies also follow a similar pattern with monotonous improvements but with slower convergence. In addition, total energy variances decrease with an increase in basis, signaling more accurate trial wave functions.

Let us now analyze the obtained fixed-node errors in more detail (Table VII). In particular, we have identified the following:

• As suggested in our previous study, the lowest percentage of fixed-node errors is observed for closed-shell states in systems with a tetrahedral arrangement of single bonds. This also holds in extended systems since almost the same error is observed in the Si crystal (to be published elsewhere).

• Somewhat larger errors are observed in open-shell systems and in excitations (Table VII) where we see, for example, an increase in percentage by a factor of 3 between the silane ground and its lowest triplet excited states.

• The largest percentage of errors does vary across the systems. For example, the two largest errors are found in the vertical excited singlet state of SiH₂ and the vertical excited triplet state of silane.

Overall, however, it is quite remarkable that in all these systems and states, the fixed-node errors are between 1% and 3.6%, especially when we consider that single-reference wave functions have been used to approximate the nodal hypersurfaces.

It has been conjectured and can be proved in a few cases that the correct number of nodal domains for fermionic ground states is 2.^48–50 This implies that all even permutations form a simply connected domain in the space of electron coordinates. In this domain, the wave function sign is constant. Similarly, the odd permutations’ complementary spatial domain mirrors this with the opposite wave function sign. Consequently, in dimensions two and higher, the fermionic ground states generically exhibit a bisection of the particle configuration space with the boundary given by the node. (There are exceptions, but a further elaboration on this is out of the scope of this work.) For excitations, the situation is a bit more complicated; the nodal count can be two or higher depending on symmetries and other characteristics. It is therefore interesting to point out that although the obtained fixed-node biases are rather small, the nodes of most of our calculations with single-reference wave functions are topologically imperfect. In fact, for most of the states studied here, the count of nodal domains is 4. This is due to the single-reference trial wave functions with their anti-symmetric parts given as a product of spin-up times and spin-down determinants.

Interestingly, the following few cases happen to have the correct nodal topologies:

• ground state of the Si atom that has only two nodal domains since there is only one minority spin electron, i.e., the exchange in that channel is absent.

• SiH₂ and SiH₄ open-shell singlet excitations that are given as a linear combination of two-determinantal products. This form breaks the spin-up and spin-down artificial separation of spatial variables and leads to the fusion of the nodal domains into the minimal two, as is illustrated in Fig. 1, and can be checked numerically.

It is rather unexpected that there is no obvious significant difference between the fixed-node biases for correct (two domains) vs incorrect (four domains) topologies. Based on these observations, we suppose that in our systems, the spin-up and spin-down coordinate subspaces are energetically marginal, and the nodal shape in other locations is more important. Obviously, with enough variational freedom, the corresponding correct topologies will ultimately emerge as a result of appropriate optimizations. Unfortunately, our knowledge about the convergence of variational expansions in this respect is rather limited although some promising progress has been achieved very recently with CIPSI methods.^51–54 Therefore, our results can serve as a reference for future studies that would address this particular issue.

One of our goals was to shed light on how the fixed-node errors impact the excited states. Therefore, we have chosen several types of excitations: singlets vs triplets, low-lying states with excited energies on the order of 1 eV vs Rydberg-like states with excitations around 10 eV. Interestingly, we do not observe any clear, simple, or regular patterns. For SiH₂, which has a diradical ground state, we find mild decreases in the fixed-node biases. Presumably, the radical character of the state is still present, and therefore, the corresponding error persists into excitations with a minor decrease due to orbital restructuring. On the other hand, for strongly bonded systems with saturated single bonds, the increase is significant, between a factor of 2 and 3. This is clearly seen in silane that exhibits both the lowest percentage bias (in the ground state) and the highest percentage bias in its excited state. On an absolute scale, the increases are not substantial since they are of the order of 0.1 eV–0.2 eV. This is quite important for a qualitative assessment of the fixed-node biases in cases where it is, at present, very difficult to estimate their actual values. In particular, this is of significant interest for excitations in periodic (crystal) systems where QMC calculations involve hundreds of valence electrons.⁵⁵ 

We emphasize that it is straightforward to eliminate the fixed-node errors in the studied systems by employing multi-reference trial wave functions. With appropriate effort, a significant decrease in the biases can be achieved using CIPSI expansions, further boosted through trial wave function re-optimizations. However, our primary goal was to understand the behavior of single-reference fixed-node errors due to our interest in large systems where a conclusive, high-accuracy multi-reference study might not be feasible.

We have presented a high-accuracy study of Si_xH_y systems using both the QMC and many-body wave function methods based on expansions in basis sets. We have probed mainly for the level of fixed-node errors in ground states and especially in excited states. We have found that the single-reference trial wave functions in FN-DMC provide excellent descriptions of these systems with 1%–3.5% correlation variation in fixed-node biases. In general, for excitations, these errors were on the larger side but not always, and variations in fixed-node biases within the mentioned small range did not show any clear regular or discernible pattern. We have found excellent agreement with previous independent all-electron state-of-the-art studies, thus providing a clear support for the accuracy of ccECPs used in this study. In silane, the vertical singlet excitation appears to be larger in our calculations than in some experiments; however, we believe that new experiments with a more clear distinction between the vertical and adiabatic excitations would be very useful. Overall, we demonstrate remarkably small fixed-node biases that further confirm our previous finding. The small fixed-node errors are related to particular properties of s and p states in the main group elements of the second row. Similar results can be expected for the same group, such as Ge, Sn, Pb as well as P, As, Sb, and Bi especially in molecules and solids with similar single-bond patterns and closed-shell states.

The obtained data will also be useful for estimations of systematic errors of QMC results in similar systems with comparable chemical bonding and electronic state calculations.

The present work used the QWalk package and corresponding tools (55% of the effort) and was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences (BES), under Award No. de-sc0012314. The rest of this work has been using the QMCPACK package (45% of the effort) and was funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, as part of the Computational Materials Sciences Program and Center for Predictive Simulation of Functional Materials. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357. This research also used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC05-00OR22725. The project used also NERSC computational time allocations and resources. The authors would like to thank Anouar Benali and Anthony Scemama for their kind help with Quantum Package, comments, and discussions.

The input/output files and supporting data generated in this work are published in Materials Data Facility^56,57 and can be found in Ref. 58.