We present high-accuracy correlated calculations of small SixHy molecular systems in both the ground and excited states. We employ quantum Monte Carlo (QMC) together with a variety of many-body wave function approaches based on basis set expansions. The calculations are carried out in a valence-only framework using recently derived correlation consistent effective core potentials. Our primary goal is to understand the fixed-node diffusion QMC errors in both the ground and excited states with single-reference trial wave functions. Using a combination of methods, we demonstrate the very high accuracy of the QMC atomization energies being within ≈0.07 eV or better when compared with essentially exact results. By employing proper choices for trial wave functions, we have found that the fixed-node QMC biases for total energies are remarkably uniform ranging between 1% and 3.5% with absolute values at most ≈0.2 eV across the systems and several types of excitations such as singlets and triplets as well as low-lying and Rydberg-like states. Our results further corroborate that Si systems, and presumably also related main group IV and V elements of the periodic table (Ge, Sn, etc), exhibit some of the lowest fixed-node biases found in valence-only electronic structure QMC calculations.
I. INTRODUCTION
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 few4,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.6 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 ccECPs8–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 SixHy 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.14 Surprisingly, our previous results suggested that this is true even though only single-reference trial wave functions have been employed.14 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.
II. METHODS
We employ our recently generated ccECPs for Si9 and H.10 Previously, extensive ccECP transferability tests on atomic spectra and on SiO, Si2, H2 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.15 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 as16
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.17 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 (Ndet). 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 SiH4, 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−1 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 results1.
All molecular geometries in this work are obtained from experimental values.20–22 Table I provides these geometries. Software packages MOLPRO23 and MRCC24 were used for CC calculations. PySCF,25 GAMESS,26 and QUANTUM packages27 were employed to generate DFT orbitals and CI expansions. Subsequent DMC calculations were performed using QWalk28 and QMCPACK29,30 packages.
SixHy molecular geometries used in this work. Experimental values are used.20–22 All values are in Å.
System . | Atom . | x . | y . | z . |
---|---|---|---|---|
SiH2 | Si | 0.000 000 | 0.000 000 | 0.000 000 |
(GS and vertical EX) | H | 0.000 000 | 0.000 000 | 1.514 020 |
H | −1.513 113 | 0.000 000 | −0.052 390 | |
SiH2 | Si | 0.000 000 | 0.000 000 | 0.000 000 |
(adiabatic EX 1B1) | H | 0.000 000 | 0.000 000 | 1.485 320 |
H | −1.253 519 | 0.000 000 | −0.796 785 | |
SiH2 | Si | 0.000 000 | 0.000 000 | 0.000 000 |
(adiabatic EX 3B1) | H | 0.000 000 | 0.000 000 | 1.476 800 |
H | −1.300 289 | 0.000 000 | −0.700 133 | |
SiH4 | Si | 0.000 0 | 0.000 0 | 0.000 0 |
H | 0.854 4 | 0.8544 | 0.8544 | |
H | −0.854 4 | −0.8544 | 0.8544 | |
H | −0.854 4 | 0.8544 | −0.8544 | |
H | 0.854 4 | −0.8544 | −0.8544 | |
Si2H6 | Si | 0.000 0 | 0.0000 | 1.1600 |
Si | 0.000 0 | 0.0000 | −1.1600 | |
H | 0.000 0 | 1.3865 | 1.6483 | |
H | −1.200 8 | −0.6933 | 1.6483 | |
H | 1.200 8 | −0.6933 | 1.6483 | |
H | 0.000 0 | −1.3865 | −1.6483 | |
H | −1.200 8 | 0.6933 | −1.6483 | |
H | 1.200 8 | 0.6933 | −1.6483 |
System . | Atom . | x . | y . | z . |
---|---|---|---|---|
SiH2 | Si | 0.000 000 | 0.000 000 | 0.000 000 |
(GS and vertical EX) | H | 0.000 000 | 0.000 000 | 1.514 020 |
H | −1.513 113 | 0.000 000 | −0.052 390 | |
SiH2 | Si | 0.000 000 | 0.000 000 | 0.000 000 |
(adiabatic EX 1B1) | H | 0.000 000 | 0.000 000 | 1.485 320 |
H | −1.253 519 | 0.000 000 | −0.796 785 | |
SiH2 | Si | 0.000 000 | 0.000 000 | 0.000 000 |
(adiabatic EX 3B1) | H | 0.000 000 | 0.000 000 | 1.476 800 |
H | −1.300 289 | 0.000 000 | −0.700 133 | |
SiH4 | Si | 0.000 0 | 0.000 0 | 0.000 0 |
H | 0.854 4 | 0.8544 | 0.8544 | |
H | −0.854 4 | −0.8544 | 0.8544 | |
H | −0.854 4 | 0.8544 | −0.8544 | |
H | 0.854 4 | −0.8544 | −0.8544 | |
Si2H6 | Si | 0.000 0 | 0.0000 | 1.1600 |
Si | 0.000 0 | 0.0000 | −1.1600 | |
H | 0.000 0 | 1.3865 | 1.6483 | |
H | −1.200 8 | −0.6933 | 1.6483 | |
H | 1.200 8 | −0.6933 | 1.6483 | |
H | 0.000 0 | −1.3865 | −1.6483 | |
H | −1.200 8 | 0.6933 | −1.6483 | |
H | 1.200 8 | 0.6933 | −1.6483 |
III. RESULTS AND DATA
A. Total energies for ground states and excitations
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 SiH2 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.20 The total energies and vertical gaps for SiH4 and Si2H6 are listed in Tables III and IV, respectively.
SiH2 total energies (Ha) and excitation gaps (eV) using various methods. All geometries were adopted from experiments.20
. | GS . | Vertical excitation . | Adiabatic excitation . | ||||||
---|---|---|---|---|---|---|---|---|---|
. | (1A1) . | (1B1) . | (3B1) . | (1B1) . | (3B1) . | ||||
. | Total . | Total . | Gap . | Total . | Gap . | Total . | Gap . | Total . | Gap . |
Method . | (hartree) . | (hartree) . | (eV) . | (hartree) . | (eV) . | (hartree) . | (eV) . | (hartree) . | (eV) . |
CISD/RHFa | –4.994 38 | –4.83224 | 4.4121 | –4.95106 | 1.1788 | –4.84958 | 3.9402 | –4.96533 | 0.7905 |
CISDT/RHFb | –4.998 62 | –4.90946 | 2.4262 | –4.95583 | 1.1644 | –4.92695 | 1.9502 | –4.96979 | 0.7845 |
exFCI/NatOrbc | –5.006 76 | –4.91838 | 2.4051 | –4.96030 | 1.2642 | –4.93527 | 1.9453 | –4.97397 | 0.8924 |
RCCSD(T)/RHFd | –5.006 27(9) | –4.95958(8) | 1.270(3) | –4.97339(8) | 0.895(3) | ||||
CCSDT(Q)/RHFe | –5.007 26(9) | –4.96088(8) | 1.262(3) | –4.97456(8) | 0.890(3) | ||||
DMC/RHF | –5.001 7(1) | –4.9091(1) | 2.520(4) | –4.9553(1) | 1.263(4) | –4.9279(1) | 2.008(4) | –4.9697(1) | 0.871(4) |
DMC/PBEf | –5.002 6(1) | –4.9136(1) | 2.422(4) | –4.9579(1) | 1.216(4) | –4.9317(1) | 1.929(4) | –4.9718(1) | 0.838(4) |
Yao et al.g | 0.897 | ||||||||
Experiment | 2.36(10)h | 1.27(10)g | 1.92768(1)i | 0.91(4)j |
. | GS . | Vertical excitation . | Adiabatic excitation . | ||||||
---|---|---|---|---|---|---|---|---|---|
. | (1A1) . | (1B1) . | (3B1) . | (1B1) . | (3B1) . | ||||
. | Total . | Total . | Gap . | Total . | Gap . | Total . | Gap . | Total . | Gap . |
Method . | (hartree) . | (hartree) . | (eV) . | (hartree) . | (eV) . | (hartree) . | (eV) . | (hartree) . | (eV) . |
CISD/RHFa | –4.994 38 | –4.83224 | 4.4121 | –4.95106 | 1.1788 | –4.84958 | 3.9402 | –4.96533 | 0.7905 |
CISDT/RHFb | –4.998 62 | –4.90946 | 2.4262 | –4.95583 | 1.1644 | –4.92695 | 1.9502 | –4.96979 | 0.7845 |
exFCI/NatOrbc | –5.006 76 | –4.91838 | 2.4051 | –4.96030 | 1.2642 | –4.93527 | 1.9453 | –4.97397 | 0.8924 |
RCCSD(T)/RHFd | –5.006 27(9) | –4.95958(8) | 1.270(3) | –4.97339(8) | 0.895(3) | ||||
CCSDT(Q)/RHFe | –5.007 26(9) | –4.96088(8) | 1.262(3) | –4.97456(8) | 0.890(3) | ||||
DMC/RHF | –5.001 7(1) | –4.9091(1) | 2.520(4) | –4.9553(1) | 1.263(4) | –4.9279(1) | 2.008(4) | –4.9697(1) | 0.871(4) |
DMC/PBEf | –5.002 6(1) | –4.9136(1) | 2.422(4) | –4.9579(1) | 1.216(4) | –4.9317(1) | 1.929(4) | –4.9718(1) | 0.838(4) |
Yao et al.g | 0.897 | ||||||||
Experiment | 2.36(10)h | 1.27(10)g | 1.92768(1)i | 0.91(4)j |
CBS values using TZ-5Z extrapolation.
CBS values using DZ-QZ extrapolation.
Ndet ≈ 200 K. Natural orbitals (NatOrbs) are obtained from smaller runs with Ndet ≈ 10 K. CBS values using DZ-QZ extrapolation.
CBS values using TZ-6Z extrapolation.
Estimated from an energy decrease compared to RCCSD(T) at the TZ basis set level.
The lowest GS energy corresponds to DMC/PBE0.
CBS values using DZ-5Z extrapolation with the SHCI method from Ref. 31.
Estimated from Fig. 9 in Ref. 32.
Reference 20.
Reference 33.
SiH4 total energies (Ha) and vertical excitation gaps (eV) using various methods. The geometries were obtained from experimental values.21 The number of determinants is given in thousands (k). Note that, counterintuitively, DMC/CIS/PBE (Ndet = 2.1 k) gives lower energy than DMC/CIPSI/NatOrb (Ndet = 360 k).
. | . | GS (1A1) . | EX (1T2) . | EX (3T2) . | ||
---|---|---|---|---|---|---|
Method . | Ndet . | Total (hartree) . | Total (hartree) . | Gap (eV) . | Total (hartree) . | Gap (eV) . |
CISD/RHFa | −6.2642 | −5.8231 | 12.003 | −5.9358 | 8.936 | |
CISD/CAS(8e−, 14o)b | −6.2770 | −5.9230 | 9.633 | |||
RCCSD(T)/RHFc | −6.2792(1) | −5.9563(1) | 8.787(4) | |||
CCSDT(Q)/RHFd | −6.2799(1) | −5.9599(1) | 8.708(4) | |||
DMC/RHF | −6.2762(1) | −5.9246(2) | 9.568(6) | −5.9465(1) | 8.972(4) | |
DMC/PBEe | −6.2777(1) | −5.9301(2) | 9.459(6) | −5.9532(1) | 8.830(4) | |
VMC/CIS/PBE(aQZ) | 2.1 k | −6.2691(3) | −5.9210(3) | 9.47(1) | ||
DMC/CIS/PBE(aQZ) | 2.1 k | −6.2767(2) | −5.9319(2) | 9.382(8) | ||
CIPSI+PT2/PBE(aQZ) | 350 k | −6.2746 | −5.9224 | 9.584 | ||
VMC/CIPSI/PBE(aQZ) | 350 k | −6.2644(2) | −5.9123(2) | 9.581(8) | ||
CIPSI+PT2/NatOrb(Aqz)f | 360 k | −6.2774 | −5.9332 | 9.366 | ||
VMC/CIPSI/NatOrb(Aqz)f | 360 k | −6.2695(1) | −5.9215(1) | 9.470(4) | ||
DMC/CIPSI/NatOrb(aQZ)f | 360 k | −6.2774(4) | −5.9302(3) | 9.45(1) | ||
Lehtonen et al.g | 9.53 | |||||
Porter et al.h | 9.44(5) | 8.89(4) | ||||
Grossman et al.i | 9.1(1) | 8.7(1) | ||||
Experiment37,38 | 8.9, 9.7 | 8.7 | ||||
Experiment39 | 9.39(15) | |||||
Experiment40 | 9.43(4) |
. | . | GS (1A1) . | EX (1T2) . | EX (3T2) . | ||
---|---|---|---|---|---|---|
Method . | Ndet . | Total (hartree) . | Total (hartree) . | Gap (eV) . | Total (hartree) . | Gap (eV) . |
CISD/RHFa | −6.2642 | −5.8231 | 12.003 | −5.9358 | 8.936 | |
CISD/CAS(8e−, 14o)b | −6.2770 | −5.9230 | 9.633 | |||
RCCSD(T)/RHFc | −6.2792(1) | −5.9563(1) | 8.787(4) | |||
CCSDT(Q)/RHFd | −6.2799(1) | −5.9599(1) | 8.708(4) | |||
DMC/RHF | −6.2762(1) | −5.9246(2) | 9.568(6) | −5.9465(1) | 8.972(4) | |
DMC/PBEe | −6.2777(1) | −5.9301(2) | 9.459(6) | −5.9532(1) | 8.830(4) | |
VMC/CIS/PBE(aQZ) | 2.1 k | −6.2691(3) | −5.9210(3) | 9.47(1) | ||
DMC/CIS/PBE(aQZ) | 2.1 k | −6.2767(2) | −5.9319(2) | 9.382(8) | ||
CIPSI+PT2/PBE(aQZ) | 350 k | −6.2746 | −5.9224 | 9.584 | ||
VMC/CIPSI/PBE(aQZ) | 350 k | −6.2644(2) | −5.9123(2) | 9.581(8) | ||
CIPSI+PT2/NatOrb(Aqz)f | 360 k | −6.2774 | −5.9332 | 9.366 | ||
VMC/CIPSI/NatOrb(Aqz)f | 360 k | −6.2695(1) | −5.9215(1) | 9.470(4) | ||
DMC/CIPSI/NatOrb(aQZ)f | 360 k | −6.2774(4) | −5.9302(3) | 9.45(1) | ||
Lehtonen et al.g | 9.53 | |||||
Porter et al.h | 9.44(5) | 8.89(4) | ||||
Grossman et al.i | 9.1(1) | 8.7(1) | ||||
Experiment37,38 | 8.9, 9.7 | 8.7 | ||||
Experiment39 | 9.39(15) | |||||
Experiment40 | 9.43(4) |
CBS values using DZ-QZ extrapolation.
CBS values using TZ-5Z extrapolation.
CBS values using TZ-6Z extrapolation.
Estimated from an energy decrease compared to RCCSD(T) at the DZ basis set level.
The EX(1T2) state uses EX(3T2) UKS/PBE orbitals, which resulted in lower energy than using GS orbitals. Spin contamination in UKS was observed to be negligible (⟨2S + 1⟩ = 3.002 083 1).
Natural orbitals (NatOrbs) are obtained from a run with Ndet ≈ 350 K.
AE CC2 calculation from Ref. 34.
DMC/CASSCF calculation from Ref. 35.
DMC/CASSCF calculation from Ref. 36.
Si2H6 total energies (Ha) and vertical excitation gaps (eV) using various methods. Experimental geometry was used for this molecule.22
. | GS (1A1g) . | EX (1E1u) . | EX (3A1g) . | ||
---|---|---|---|---|---|
Method . | Total (hartree) . | Total (hartree) . | Gap (eV) . | Total (hartree) . | Gap (eV) . |
CISD/RHFa | −11.3287 | −10.8183 | 13.889 | −11.0825 | 6.699 |
CISD/CAS(14e−, 13o)b | −11.3400 | −11.0623 | 7.557 | ||
RCCSD(T)/RHFc | −11.3766(3) | −11.1308(5) | 6.69(2) | ||
CCSDT(Q)/RHFd | −11.3782(3) | −11.1336(5) | 6.66(2) | ||
DMC/RHF | −11.3708(2) | −11.0855(2) | 7.763(8) | −11.1215(2) | 6.784(8) |
DMC/PBEe | −11.3725(2) | −11.0934(2) | 7.595(8) | −11.1248(2) | 6.740(8) |
Lehtonen et al.f | 7.61 | ||||
Experiment41 | 7.6 | ≈6.7g | |||
Experiment39 | 7.56 |
. | GS (1A1g) . | EX (1E1u) . | EX (3A1g) . | ||
---|---|---|---|---|---|
Method . | Total (hartree) . | Total (hartree) . | Gap (eV) . | Total (hartree) . | Gap (eV) . |
CISD/RHFa | −11.3287 | −10.8183 | 13.889 | −11.0825 | 6.699 |
CISD/CAS(14e−, 13o)b | −11.3400 | −11.0623 | 7.557 | ||
RCCSD(T)/RHFc | −11.3766(3) | −11.1308(5) | 6.69(2) | ||
CCSDT(Q)/RHFd | −11.3782(3) | −11.1336(5) | 6.66(2) | ||
DMC/RHF | −11.3708(2) | −11.0855(2) | 7.763(8) | −11.1215(2) | 6.784(8) |
DMC/PBEe | −11.3725(2) | −11.0934(2) | 7.595(8) | −11.1248(2) | 6.740(8) |
Lehtonen et al.f | 7.61 | ||||
Experiment41 | 7.6 | ≈6.7g | |||
Experiment39 | 7.56 |
CBS values using DZ-QZ extrapolation.
CBS values using DZ-QZ extrapolation.
CBS values using TZ-6Z extrapolation.
Estimated from an energy decrease compared to RCCSD(T) at the DZ basis set level.
The lowest EX(3A1g) energy corresponds to DMC/PBE0.
AE CC2 values from Ref. 34.
Experimental value in Ref. 36 of 6.3 eV is probably a mix of vertical and adiabatic excitations.
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 SiH2 and Si2H6, 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 SiH2, 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 SiH4, we observe higher DMC/PBE excitation gap discrepancies with some experiments for the vertical states. Inspection of the experiments37 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.
B. Binding energies
Table V presents the SixHy 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.
SixHy atomization energies (eV) using various methods calculated using RCCSD(T), CCSDT(Q), and DMC/DFT methods compared with previous correlated calculations35,42–45 and with experiments. The experimental values are corrected to correspond to the bottom of the well (static ions) atomization energies. Experimental values of zero-point energy have been used (SiH2: 0.3092 eV, SiH4: 0.8339 eV, and Si2H6: 1.3042 eV).46
System . | Method . | De (eV) . |
---|---|---|
SiH2 | RCCSD(T) | 6.671(3) |
CCSDT(Q) | 6.673(3) | |
DMC/PBE0 | 6.599(4) | |
Feller et al.a | 6.66506 | |
Haunschild et al.b | 6.665 | |
Greef et al.c | 6.57(2) | |
Yao et al.d | 6.685 | |
Experimentse | 6.40(9), 6.7(1) | |
SiH4 | RCCSD(T) | 14.098(3) |
CCSDT(Q) | 14.091(3) | |
DMC/PBE | 14.085(4) | |
Feller et al.a | 14.0847 | |
Haunschild et al.b | 14.074 | |
Greef et al.c | 14.05(2) | |
Martin et al.f | 14.0235 | |
Porter et al.g | 14.19(2) | |
Yao et al.d | 14.107 | |
Experimentse | 13.96(2), 14.00(2) | |
Si2H6 | RCCSD(T) | 23.248(9) |
CCSDT(Q) | 23.241(9) | |
DMC/PBE | 23.192(8) | |
Feller et al.a | 23.2345 | |
Haunschild et al.b | 23.222 | |
Greef et al.c | 23.20(3) | |
Yao et al.d | 23.255 | |
Experimentse | 22.991(4), 23.08(1) |
System . | Method . | De (eV) . |
---|---|---|
SiH2 | RCCSD(T) | 6.671(3) |
CCSDT(Q) | 6.673(3) | |
DMC/PBE0 | 6.599(4) | |
Feller et al.a | 6.66506 | |
Haunschild et al.b | 6.665 | |
Greef et al.c | 6.57(2) | |
Yao et al.d | 6.685 | |
Experimentse | 6.40(9), 6.7(1) | |
SiH4 | RCCSD(T) | 14.098(3) |
CCSDT(Q) | 14.091(3) | |
DMC/PBE | 14.085(4) | |
Feller et al.a | 14.0847 | |
Haunschild et al.b | 14.074 | |
Greef et al.c | 14.05(2) | |
Martin et al.f | 14.0235 | |
Porter et al.g | 14.19(2) | |
Yao et al.d | 14.107 | |
Experimentse | 13.96(2), 14.00(2) | |
Si2H6 | RCCSD(T) | 23.248(9) |
CCSDT(Q) | 23.241(9) | |
DMC/PBE | 23.192(8) | |
Feller et al.a | 23.2345 | |
Haunschild et al.b | 23.222 | |
Greef et al.c | 23.20(3) | |
Yao et al.d | 23.255 | |
Experimentse | 22.991(4), 23.08(1) |
AE UCCSD(T) calculations from Ref. 42.
AE frozen core CCSD(T) method calculations from Ref. 43.
DMC calculations from Ref. 45.
CBS values using DZ-5Z extrapolation with the SHCI method.31
Experimental values summarized in Ref. 42.
AE CCSD(T) calculation from Ref. 44.
DMC/CASSCF calculation from Ref. 35.
Included also are results from published high-accuracy studies42–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.
C. Fixed-node errors
1. Fixed-node biases for single-reference trial wave functions
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.47). 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.
Impact of orbitals and basis sets on fixed-node DMC energies (Ha) for SiH4 GS and vertical excitation in the aug-cc-pVnZ basis set. Total energies, kinetic energies, and total energy variances are shown for the Slater–Jastrow wave function (single-reference for GS, two-reference for EX).
Basis . | GS . | EX . | ||||
---|---|---|---|---|---|---|
nZ . | Total . | Kinetic . | Variance . | Total . | Kinetic . | Variance . |
DZ | −6.2740(2) | 3.9521(28) | 0.0589(3) | −5.9229(3) | 3.7568(16) | 0.0559(2) |
TZ | −6.2765(1) | 3.9658(20) | 0.0507(1) | −5.9252(3) | 3.7651(27) | 0.0468(2) |
QZ | −6.2774(1) | 3.9690(22) | 0.0428(2) | −5.9277(2) | 3.7650(36) | 0.0392(2) |
5Z | −6.2776(1) | 3.9771(15) | 0.0365(1) | −5.9279(2) | 3.7683(26) | 0.0346(2) |
6Z | −6.2777(1) | 3.9779(13) | 0.0355(1) | −5.9272(2) | 3.7625(23) | 0.0330(1) |
Basis . | GS . | EX . | ||||
---|---|---|---|---|---|---|
nZ . | Total . | Kinetic . | Variance . | Total . | Kinetic . | Variance . |
DZ | −6.2740(2) | 3.9521(28) | 0.0589(3) | −5.9229(3) | 3.7568(16) | 0.0559(2) |
TZ | −6.2765(1) | 3.9658(20) | 0.0507(1) | −5.9252(3) | 3.7651(27) | 0.0468(2) |
QZ | −6.2774(1) | 3.9690(22) | 0.0428(2) | −5.9277(2) | 3.7650(36) | 0.0392(2) |
5Z | −6.2776(1) | 3.9771(15) | 0.0365(1) | −5.9279(2) | 3.7683(26) | 0.0346(2) |
6Z | −6.2777(1) | 3.9779(13) | 0.0355(1) | −5.9272(2) | 3.7625(23) | 0.0330(1) |
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 SiH2 and the vertical excited triplet state of silane.
Fixed-node error analysis from single-reference DMC and CC energies (Ha). For comparison, we have also included preliminary results for Si crystal, to be published elsewhere. MAD denotes mean average deviation.
. | . | DMC/DFTa . | SCF/RHF . | Estim. Exact . | FN err. . | FN err./Corr. . | FN err./Nelec . |
---|---|---|---|---|---|---|---|
System . | State . | (Ha) . | (Ha) . | Corr. (Ha) . | (eV) . | (%) . | (eV/electron) . |
Si | 3P | −3.7601(1) | −3.672 4778(1) | −0.08957(6) | 0.053(3) | 2.2(1) | 0.0133(8) |
SiH2 | 1A1 | −5.0026(1) | −4.853 64(8) | −0.1536(1) | 0.127(4) | 3.03(9) | 0.0211(6) |
SiH2 | 1B1(vert) | −4.9136(1) | −4.770 51(7)b | −0.14787(7) | 0.130(3) | 3.23(7) | 0.0217(5) |
SiH2 | 1B1(adia) | −4.9317(1) | −4.792 82(7)b | −0.14245(7) | 0.097(3) | 2.51(7) | 0.0162(5) |
SiH2 | 3B1(vert) | −4.9579(1) | −4.831 00(7) | −0.1299(1) | 0.081(3) | 2.3(1) | 0.0135(6) |
SiH2 | 3B1(adia) | −4.9718(1) | −4.846 33(7) | −0.1282(1) | 0.075(3) | 2.2(1) | 0.0125(6) |
SiH4 | 1A1 | −6.2777(1) | −6.088 0(1) | −0.1919(1) | 0.060(4) | 1.15(7) | 0.0075(5) |
SiH4 | 3T2 | −5.9532(1) | −5.773 5(1) | −0.1864(1) | 0.182(4) | 3.59(7) | 0.0228(5) |
Si2H6 | 1A1g | −11.3725(2) | −11.0232 (3) | −0.3550(4) | 0.16(1) | 1.6(1) | 0.0111(7) |
Si2H6 | 3B1g | −11.1248(2) | −10.787 8(4) | −0.3458(7) | 0.24(1) | 2.5(2) | 0.017(1) |
Si cryst/atomc | −0.14342(2) | 0.05(1) | 1.3(3) | 0.013(3) | |||
MAD ground states | 0.099(3) | 1.99(5) | 0.0132(3) | ||||
MAD excited states | 0.134(3) | 2.72(4) | 0.0173(3) |
. | . | DMC/DFTa . | SCF/RHF . | Estim. Exact . | FN err. . | FN err./Corr. . | FN err./Nelec . |
---|---|---|---|---|---|---|---|
System . | State . | (Ha) . | (Ha) . | Corr. (Ha) . | (eV) . | (%) . | (eV/electron) . |
Si | 3P | −3.7601(1) | −3.672 4778(1) | −0.08957(6) | 0.053(3) | 2.2(1) | 0.0133(8) |
SiH2 | 1A1 | −5.0026(1) | −4.853 64(8) | −0.1536(1) | 0.127(4) | 3.03(9) | 0.0211(6) |
SiH2 | 1B1(vert) | −4.9136(1) | −4.770 51(7)b | −0.14787(7) | 0.130(3) | 3.23(7) | 0.0217(5) |
SiH2 | 1B1(adia) | −4.9317(1) | −4.792 82(7)b | −0.14245(7) | 0.097(3) | 2.51(7) | 0.0162(5) |
SiH2 | 3B1(vert) | −4.9579(1) | −4.831 00(7) | −0.1299(1) | 0.081(3) | 2.3(1) | 0.0135(6) |
SiH2 | 3B1(adia) | −4.9718(1) | −4.846 33(7) | −0.1282(1) | 0.075(3) | 2.2(1) | 0.0125(6) |
SiH4 | 1A1 | −6.2777(1) | −6.088 0(1) | −0.1919(1) | 0.060(4) | 1.15(7) | 0.0075(5) |
SiH4 | 3T2 | −5.9532(1) | −5.773 5(1) | −0.1864(1) | 0.182(4) | 3.59(7) | 0.0228(5) |
Si2H6 | 1A1g | −11.3725(2) | −11.0232 (3) | −0.3550(4) | 0.16(1) | 1.6(1) | 0.0111(7) |
Si2H6 | 3B1g | −11.1248(2) | −10.787 8(4) | −0.3458(7) | 0.24(1) | 2.5(2) | 0.017(1) |
Si cryst/atomc | −0.14342(2) | 0.05(1) | 1.3(3) | 0.013(3) | |||
MAD ground states | 0.099(3) | 1.99(5) | 0.0132(3) | ||||
MAD excited states | 0.134(3) | 2.72(4) | 0.0173(3) |
DMC/DFT corresponds to the lowest single-reference DMC energies obtained in this work (two-reference for EX singlet).
SiH2 vertical and adiabatic singlet SCF/RHF energies correspond to CASSCF(2e−, 2°) energy with TZ-6Z extrapolations.
Assumed experimental cohesion of 4.68(1) eV/at. and the extrapolated DMC/PBE0 value for the Si crystal ground state.
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.
2. Impact of nodal domain topologies vs nodal shapes
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.
SiH2 and SiH4 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.
Nodes of the SiH4 molecule using 2e− (up+down electrons at the same point) scan for (a) ground state (single-reference) and (b) excited state singlet (two-reference). Quadrisections of the configurations in case (a) and bisection in case (b) are clearly visible.
Nodes of the SiH4 molecule using 2e− (up+down electrons at the same point) scan for (a) ground state (single-reference) and (b) excited state singlet (two-reference). Quadrisections of the configurations in case (a) and bisection in case (b) are clearly visible.
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.
3. Fixed-node biases in excitations
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 SiH2, 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.55
4. Multi-reference trial wave functions
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.
IV. CONCLUSIONS
We have presented a high-accuracy study of SixHy 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.
ACKNOWLEDGMENTS
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.
DATA AVAILABILITY
REFERENCES
Future work to be published on Si crystal using QMC.