# Spectroscopic accuracy directly from quantum chemistry: application to ground and excited states of beryllium dimer

###### Abstract

We combine explicit correlation via the canonical transcorrelation approach with the density matrix renormalization group and initiator full configuration interaction quantum Monte Carlo methods to compute a near-exact beryllium dimer curve, without the use of composite methods. In particular, our direct density matrix renormalization group calculations produce a well-depth of =931.2 cm which agrees very well with recent experimentally derived estimates =929.7 cm [Science, 324, 1548 (2009)] and =934.6 cm [Science, 326, 1382 (2009)]], as well the best composite theoretical estimates, =938 cm [J. Phys. Chem. A, 111, 12822 (2007)] and =935.1 cm [Phys. Chem. Chem. Phys., 13, 20311 (2011)]. Our results suggest possible inaccuracies in the functional form of the potential used at shorter bond lengths to fit the experimental data [Science, 324, 1548 (2009)]. With the density matrix renormalization group we also compute near-exact vertical excitation energies at the equilibrium geometry. These provide non-trivial benchmarks for quantum chemical methods for excited states, and illustrate the surprisingly large error that remains for 1 state with approximate multi-reference configuration interaction and equation-of-motion coupled cluster methods. Overall, we demonstrate that explicitly correlated density matrix renormalization group and initiator full configuration interaction quantum Monte Carlo methods allow us to fully converge to the basis set and correlation limit of the non-relativistic Schrödinger equation in small molecules.

## I Introduction

In basis-set quantum chemistry, we divide the challenge of solving the non-relativistic electronic Schrödinger equation into two parts: the treatment of -electron correlations, and the saturation of the one-particle basis. Recent years have seen significant advances in both these areas. In the first case, methods such as general order coupled cluster (CC) kallay , the density matrix renormalization group (DMRG) chan2011 , and initiator full configuration interaction quantum Monte Carlo (i-FCIQMC) cleland ; booth1 have been developed to achieve an efficient treatment of arbitrary -electron correlations in modestly sized molecules. In the second case, explicit correlation (F12) techniquesTen-no2004 augment the one-particle basis with geminal functions that represent the electron-electron cusp. Taken together, these advances provide the potential to converge to near-exact solutions of the non-relativistic electronic Schrödinger equation at the basis set limit. In this report, we describe the efficient combination of explicit correlation, via the canonical transcorrelation approach yanai2012canonical , with DMRG and with i-FCIQMC, and apply these combinations to determine the ground and excited state electronic structure of the beryllium dimer to very high accuracy.

## Ii Methods

The essence of explicit correlation (henceforth referred to as F12 theory) is to use a geminal correlation factor, to augment the doubles manifold of the virtual space Kong2011 ; tennoreview . The geminal can be thought of as including some excitations into a formally infinite basis of virtuals. Labelling the infinite virtual basis by , the geminal doubles excitation operator is written as

(1) |

where are the geminal doubles amplitudes. A significant practical advance was the realization that the geminal amplitudes are fixed to linear order by the electron-electron cusp condition Kutzelnigg1985 ; kut_klopp ; Ten-no2004 ,

(2) |

where is a projector, defined in terms of projectors and into the occupied and virtual space of the standard orbital basis,

(3) |

that ensures that excitations of the geminal factor are orthogonal to those of the standard orbital space sinanoglu .

Combining F12 methodology with the DMRG and i-FCIQMC methods involves practical hurdles not present in prior combinations of F12 theory with other correlation methods. For example, explicitly correlated coupled cluster theory formally starts from a non-Hermitian effective Hamiltonian, obtained by similarity transforming with the geminal excitation operator torheyden but the DMRG is most conveniently implemented with a Hermitian effective Hamiltonian. Similarly, the universal perturbative correction of Torheyden and Valeev, which has been used with i-FCIQMC booth ; torheyden2 does not introduce non-Hermiticity, but requires the one- and two-particle reduced density matrices which can be expensive to compute precisely in Monte Carlo methods. These practical complications are removed within the recently introduced canonical transcorrelation form of F12 theory of Yanai and Shiozaki yanai2012canonical . In this method, a Hermitian effective Hamiltonian is obtained from an anti-hermitian geminal doubles excitation operator,

(4) |

The canonical transcorrelated Hamiltonian is formally defined as

(5) |

The fully transformed Hamiltonian involves operators of high particle rank. To ameliorate the complexity, Yanai and Shiozaki invoke the same commutator approximations used in the canonical transformation theory Neuscamman2010 ; yanaict1 ; yanaict2 , and further simplify the quadratic commutator term by replacing the Hamiltonian with a generalized Fock operator ,

(6) |

an approximation which is valid through second-order in perturbation theory. The subscript denotes that only one- and two-particle rank operators and density matrices are kept in the Mukherjee-Kutzelnigg normal-ordered form Mukherjee1997 ; kutzelnigg97 . In our calculations here, normal-ordering is carried out with respect to the Hartree-Fock reference, thus no density cumulants Mazziottiijqc ; kutzelniggcumulant ; Mazziotticpl ; MazziottiReview appear. The only error arises from the neglect of the three-particle normal-ordered operatorEricQuadraticCT generated by the excitations. As the orbital basis increases, tends to zero and the three-particle error also goes to zero, very different behaviour from density cumulant theories where full three-particle quantity reconstruction is performed EricQuadraticCT ; Mazziottijcp . In this sense, the three-particle error in this theory is part of the basis set error.

The practical advantages of the canonical transcorrelation formulation are that no correlated density matrices are required and the effective Hamiltonian is Hermitian and of two-particle form. It may thus be combined readily with any correlation treatment. Beyond the practical advantages, the canonical transcorrelation formulation uses a “perturb then diagonalize” approach, rather than the “diagonalize then perturb” approach of the a posteriori F12 treatment of Valeev previously combined with i-FCIQMC. This allows the geminal factors to automatically relax the parameters of the subsequent correlation treatment. This approach is similar in spirit to the similarity transformed F12 method of Ten-no which has been used with the PMC-SD method of Ohtsuka ohtsuka , but there the excitation operator is not anti-hermitian and alternative approximations are used in the simplification of the resulting equations. The use of the effective Hamiltonian (6) was denoted by Yanai and Shiozaki by the prefix F12-, thus in their nomenclature, the combinations with DMRG and i-FCIQMC in this work would be F12-DMRG and F12-i-FCIQMC respectively. However, as all our DMRG and i-FCIQMC calculations use this effective Hamiltonian here, we will usually omit the F12 prefix and simply refer to DMRG and i-FCIQMC.

We now briefly introduce the DMRG and i-FCIQMC correlation methods used in this work. The DMRG is a variational ansatz based on a matrix-product representation of the FCI amplitudes. Expanding the FCI wavefunction as

(7) |

where is the occupancy of orbital in the occupancy vector representation of the -particle determinant , , the one-site DMRG wavefunction approximates the FCI coefficient as the vector, matrix, …, matrix, vector product

(8) |

For each occupancy , the dimension of the corresponding matrix(vector) is (). is usually referred to as the number of renormalised states. The energy is determined by minimizing with respect to the matrix and vector coefficients in Eq. (8) Chan2002 ; Hachmann2006 ; chan2011 ; marti11 . As is increased, the DMRG energy converges towards the FCI limit. In practical DMRG calculations, the minimization is carried out with a slightly more flexible wavefunction form, where two , matrices on adjacent orbitals are fused into a single larger composite (two-site) matrix, . This introduces a larger variational space than in Eq. (8), which improves the numerical convergence. A measure of the error in a DMRG calculation is provided by the “discarded” weight, which is the (squared) difference in overlap between the two-site wavefunction and its best one-site approximation: this difference vanishes as . The discarded weight usually exhibits a linear relationship with the energy, and thus provides a convenient way to extrapolate the energy to the exact FCI result legeza ; Chan2002 .

The FCIQMC algorithm has been recently introduced by Alavi and co-workers Booth2013 ; booth1 ; cleland ; semistochastic . FCIQMC is a projector Monte Carlo method wherein the stochastic walk is done in determinant space BlaSug-PRD-83 ; trivedi , but instead of imposing a fixed node approximation Haaf1995 it uses computational power and cancellation algorithms to control the fermion sign problem. The exact ground state wavefunction is obtained by repeatedly applying a “projector” to an initial state,

(9) |

If is expanded in an orthogonal basis of determinants, the expansion coefficients evolve according to

(10) |

where labels the iterations, and is a time step, the maximum value of which is constrained by the inverse of the spectral range of the Hamiltonian. Since the number of basis states, is too large to permit storing all the coefficients, , a stochastic approach is used wherein “walkers” () sample the wavefunction. Although the distribution of walkers among the states at any time step is a crude approximation to the wavefunction, the infinite time average yields the ground state wavefunction exactly. The term in Eq. (10) leads to an increase or decrease in the weight of the walker on determinant while causes transitions of walkers from determinant to determinant . If walkers land on the same determinant, their weights are combined. However, because contributions to a given determinant can be of either sign for most systems, a fermionic sign problem results, where the signal becomes exponentially small compared to the noisefoulkes . As demonstrated by Alavi and coworkers, however, when cancellations are employed, for sufficiently large , the walk undergoes a transition into a regime where the sign problem is controlledbooth1 .

For the sufficiently large such that this cancellation is effective, FCIQMC is exact within a statistical error of order . However, the cost of this brute-force approach prevents application to realistic problems. A significant advance was the introduction of the initiator approximation (i-FCIQMC) booth1 ; cleland ; cleland11 . In the initiator approximation, only walkers beyond a certain initiator threshold are allowed to generate walkers on the unsampled determinants. The result is that low-weight determinants whose sign may not be sufficiently accurate, propagate according to a dynamically truncated hamiltonian, defined by the space of instantaneously occupied determinants. This concentrates the stochastic walk within a subspace of the full Hilbert space allowing for more effective cancellation, at the cost of introducing an initiator error that may be either positive or negative. However, as the total number of walkers is increased (for fixed ) the i-FCIQMC energy converges to the FCI limit. Additional large efficiency gains can be made by carrying out some of the walk non-stochastically and by using a multi-determinantal trial wave function when computing the energy estimator, giving rise to semistochastic quantum Monte Carlo semistochastic . This is not used in the results presented here, but, future studies will investigate the gain in efficiency and the possible reduction in initiator bias from doing so.

## Iii Results and Discussion

We now describe the application of the DMRG and i-FCIQMC methods to the beryllium dimer. The beryllium dimer has been of long-standing interest to theory and experiment. (See Refs. Patkowski1 ; Merritt for an overview of earlier theoretical and experimental work). Simple molecular orbital arguments would say that the molecule is unbound, however, Be can in fact be observed in the gas phase. The observed bond is significantly stronger than that of other van der Waal’s closed shell diatomics such as He and Ne Patkowski1 ; Patkowski2 . The unusual bonding arises from electron correlation effects that are enhanced by the near degeneracy of the Be atom. This near-degeneracy, coupled with the need for very large basis sets to describe the long bond-lengths, presents a challenge for modern electronic structure methods, while the weak bond makes accurate experimental measurement challenging. The lack of accurate theoretical data has also hindered the intepretation of experiment, as the a priori assumed functional form of the potential energy curve biases the extraction of parameters from the spectral lines. Thus, for many years, there had been significant disagreement between theory and experiment.

The earliest experimental estimate of the well-depth () was 79030 cm (Ref. Bondybey, ; Bondybey2, ), but this used a Morse potential in the fitting that has the wrong shape at large distances, where van der Waal’s forces dominate. Theoretical calculations generally yielded much deeper wells. Composite coupled cluster/full-configuration interaction schemes that sum over core/valence (CV), complete basis set (CBS), high-order correlation effects, and relativistic corrections, gave as 944 cm (Ref. Martin, ), 93815 cm (Ref. Patkowski1, ) and 93510 cm (Ref. koput, ). We believe the latter calculation to be the most accurate to date. Variants of multireference configuration interaction gave similar, but slightly shallower wells: 9038 cm (Ref. Gdanitz1999, , -MR-ACPF with relativistic corrections), 912 cm (Ref. Schmidt2010, , MRCI with CV, CBS, and relativistic corrections), and 923 cm (Ref. koput, , MRCI+Q, no error bar). Only recently, remeasurements by Merritt et al.Merritt , together with an improved fitting of the experimental spectrum, yielded an experimentally derived consistent with theory: 929.72.0 cm, which lies within the error bars of the calculations. A further refit of Merritt et al.’s measurements to a “fine-tuned” version of the potential of Ref. Patkowski2, gave a slightly modified well-depth of 934.6 cm, presumably with similar error bars to Ref. Merritt, . This can be regarded as the most accurate “experimental” estimate of to date.

With the recent resolution of the disagreement between theory and experiment, bonding in the beryllium dimer can now be considered to be satisfactorily understood, at least from a computational perspective. Nonetheless, the theoretical efforts so far have required careful composite schemes to separately saturate basis set effects, high-order correlation, and core contributions. While such additive schemes perform quite well, the need to assume additivity between large contributions is theoretically unsatisfactory and can potentially introduce some uncertainty into the final predicted result. For example, the all-electron FCI calculation in Ref. Patkowski1, could only be carried out in an aug-cc-pVDZ basis, and gave a well-depth of only 181 cm, while the CCSD(T) calculations in the largest aug-cc-pV7Z basis koput gave a well-depth of only 696 cm. Thus, in reaching the value of 935 cm a large degree of transferability amongst incremental contributions was assumed. The only non-composite method, the -MR-ACPF calculation of Gdanitz Gdanitz1999 gave a non-relativistic cm, which remains quite far from the best experimental or theoretical results.

M | Energy | Discarded weight |
---|---|---|

500 | 1.57 | -29.338592 |

1000 | 2.28 | -29.338647 |

1500 | 5.81 | -29.338655 |

2000 | 1.56 | -29.338657 |

(l.) | - | -29.338657 |

(q.) | - | -29.338658 |

We can now carry out a direct calculation, with saturated large basis sets and explicit correlation as well as a full account of the -electron correlations, using the canonically transcorrelated DMRG and i-FCIQMC methods, thus eliminating the need for composite approaches. We have computed several points along the ground-state Be potential energy curve using a series of cc-pCVnZ-F12 basis sets basis with =D, T, Q (henceforth referred to as DZ, TZ, and QZ, for short) and cc-pCVnZ-F12_OPTRI basis basis sets with =D, T, Q respectively for the resolution of the identity (RI) basis sets. These basis sets contain 68, 124, and 192 basis functions respectively, with up to functions in the QZ basis, and the RI basis sets contain 164, 190 and 188 basis functions respectively. The DMRG calculations were carried out using the Block code sharma2012spin . This DMRG implementation incorporates two symmetries not commonly found in other implementations: spin-adaptation (SU(2)) and symmetries. Spin-adapted DMRG implementations for quantum chemistry were described by Wouters et al. wouters and our group sharma2012spin , based on earlier work by McCulloch mcculloch3 . Compared to non-spin-adapted DMRG with only symmetry, we find that calculations with spin-adapted states correspond in accuracy to approximately renormalized non-spin-adapted states in the calculation sharma2012spin . Our implementation of symmetry resembles that for spin-symmetry, where the Wigner-Eckart theorem is used to simplify the evaluation of matrix elements as well as to reduce storage. We find that symmetry brings an additional factor of 2 in the effective over the use of only symmetry. Consequently, with both spin and adaptation, our reported energies here with renormalized states are roughly comparable in accuracy to similar calculations with renormalized states in a conventional DMRG code with only and symmetries. Our calculation at the bond length of took a wall clock time of 150 hours running in parallel on 72 Intel Xeon E5-2670 cores, totalling 10,800 core hours.

The i-FCIQMC calculations were carried out using the Neci code booth1 ; booth2011 ; Booth-neci . These calculations used the Abelian rotational subgroup of , as described in Ref. booth2011, . This symmetrized determinant space is smaller than that for the group, especially for large angular momentum basis sets, but is larger than that for the full group by less than a factor of 2 because it retains only one-dimensional irreducible representations.

The F12 integrals and transcorrelated Hamiltonian were generated using the Orz code, using the F12 exponent . The well-depth was calculated from the energy at Å. All 8 electrons were correlated, thus the largest calculation formally involved more than determinants. In the DMRG calculations, we also computed the lowest 4 excited states in the class of irreps (, , , ) at the ground-state equilibrium geometry of . For comparison, we also present results of CCSD(T), CCSD(T)-F12, and F12-CCSD(T) adler2007 ; knizia ; yanai2012canonical calculations for the ground-state curve, and MRCI-F12 shiozaki2011 , MRCI, and EOM-CCSD calculations for the excited states. These computations were performed using the Molpro packagewerner2012 ; the F12-CCSD(T) calculations used the MRCC program with the transcorrelated Hamiltonian as inputmrcc .

CCSD(T) | CCSD(T)/CBS | CCSD(T)-F12b | F12-CCSD(T) | ||||
---|---|---|---|---|---|---|---|

/Å | QZ | 5Z | 6Z | QZ | QZ | ||

2.20 | 0.46 | 0.73 | 0.86 | 0.97 | 1.05 | 0.88 | 0.86 |

2.40 | 2.73 | 2.94 | 3.04 | 3.12 | 3.18 | 3.07 | 3.05 |

2.45 | 2.83 | 3.03 | 3.12 | 3.20 | 3.25 | 3.15 | 3.14 |

2.50 | 2.83 | 3.02 | 3.11 | 3.18 | 3.23 | 3.14 | 3.12 |

3.00 | 1.45 | 1.55 | 1.60 | 1.64 | 1.67 | 1.60 | 1.61 |

5.00 | 0.34 | 0.35 | 0.36 | 0.36 | 0.36 | 0.35 | 0.35 |

DMRG | i-FCIQMC | experiment | ||||||
---|---|---|---|---|---|---|---|---|

r/Å | DZ | TZ | QZ | CBS/BSSE/rel. | QZ (50) | QZ(200) | Merritt | Patkowski |

2.20 | 0.68 | 1.76 | 2.11 | 2.23(0.08) | 2.06 (0.02) | 2.41 | 2.21 | |

2.30 | 2.33 | 3.41 | – | – | – | 3.67 | 3.71 | |

2.40 | 3.02 | 3.99 | 4.20 | 4.26(0.05) | 4.22 (0.02) | 4.17 | 4.22 | |

2.45 | 3.13 | 4.09 | 4.24 | 4.30(0.04) | 4.27 (0.05) | 4.21(0.04) | 4.24 | 4.26 |

2.50 | 3.13 | 4.00 | 4.18 | 4.24(0.04) | 4.32 (0.03) | 4.11(0.04) | 4.20 | 4.19 |

2.60 | 2.95 | 3.70 | – | – | – | 3.89 | 3.86 | |

2.70 | 2.64 | 3.27 | – | – | – | 3.44 | 3.42 | |

3.00 | 1.72 | 2.12 | 2.23 | 2.26(0.02) | 2.32 (0.05) | 2.18 | 2.22 | |

5.00 | 0.31 | 0.37 | 0.39 | 0.39(0.01) | 0.55 (0.03) | 0.39 | 0.40 |

Method | |
---|---|

CCSD(T)-F12b/BSSE | 699.3 |

DMRG | 931.2 |

i-FCIQMC | 924(9) |

DMRG/CBS/BSSE/rel. | 944(10) |

Author | |

Merritt(E/T)Merritt | 929.7(2) |

Patkowski(E/T)Patkowski2 | 934.6 |

Patkowski(T)Patkowski1 | 938.0(15) |

Schmidt(T)Schmidt2010 | 915.5 |

Koput(T)koput | 935.1(10) |

Tables 3 and 4 present our accumulated data for the DMRG and i-FCIQMC ground-state Be calculations, as well as selected computed and reference data for the well-depths. All DMRG energies correspond to (see below) while all i-FCIQMC calculations were carried out with and (see below). Figure 1 shows the convergence of the DMRG energy as a function of the discarded weight and for the QZ basis at ; energies as a function of are given in Table 1. We note that the DMRG energies presented in Table 1 and Figure 1 were obtained by first carrying out standard DMRG calculations up to =2500, and then backtracking (by decreasing in subsequent sweeps) down to =500 in steps of 500, to obtain the tabulated energies at =500, 1000, 1500, 2000. This ensures that the energy at each is well converged and free from any initialization bias, leading to more accurate extrapolation. We calculate the DMRG extrapolated energy by fitting to linear and quadratic functions of the discarded weight. Due to the high cost of calculation, insufficient sweeps were performed at =2500 to attain full convergence, hence the DMRG energies at =2500 were not themselves used in the extrapolation. The maximum difference between the linear and quadratic extrapolations is 6 , and we use this as an upper estimate of the remaining error in the DMRG energy. Examining Fig. 1, we find that the DMRG energy converges extremely rapidly with : even by , the total DMRG energy in the QZ basis appears within 10 (2 cm) of the extrapolated result!

The i-FCIQMC energies contain two sources of error: statistical error (due to the finite simulation time), and initiator error (due to the finite walker population). The statistical errors are listed in the Table 3 and are on the order of 20-50 . The remaining discrepancy between the i-FCIQMC energies and the DMRG energies is due to initiator error. Note that the initiator error can be of either sign. Because of the small energy scales of this system, the initiator error is significant at some bond-lengths. For example, at , the initiator error with is 0.14 m, or about , causing the i-FCIQMC curve to have an unphysical shape (the energy at is below that at the equilibrium distance ). The initiator error can be removed by carrying out simulations with larger number of walkers. At and we recomputed the i-FCIQMC using walkers. These i-FCIQMC are now in better agreement with the converged DMRG energies and restore the physical shape of the potential. However, such calculations were 3-4 times more expensive than the corresponding DMRG calculations.

We now discuss the possible remaining sources of error and non-optimality in our calculations. These include basis set superposition error (BSSE), relativistic effects, non-optimality of the F12 exponent, geometry effects, errors associated with the F12 approximations in the canonical transcorrelation approach and basis set incompleteness error. BSSE error can be estimated from the counterpoise correction boys . We find the counterpoise contribution to the F12-DMRG well-depth to be -11 (-2.4 cm) at the QZ level. Our relativistic correction using the CCSD(T)/aug-cc-pCVQZ method with the second-order Douglas-Kroll-Hess (DKH) one-electron Hamiltonian is -4.2 cm, which is in good agreement with previous studiesPatkowski1 ; Gdanitz1999 . We have checked the optimality of the F12 exponent and the bond-length effects through CCSD(T)-F12 calculations adler2007 ; knizia . At the QZ level, yielded the same CCSD(T)-F12 =3.2 m to within 2 (0.4 cm) and thus we conclude that our exponent of is near-optimal. The difference in energy between the CCSD(T)-F12/QZ equilibrium bond-length energy (at 2.46Å), and the energy at our assumed is only 3 (0.6 cm).

The F12 canonical transcorrelation approach contains two kinds of error. The first is the auxiliary basis integral approximations used to compute the F12 integrals, and the second is the neglect of normal-ordered three-particle operators in the canonical transcorrelated Hamiltonian as described above. (We recall that in this work all three-particle cumulants are zero in our definition of , since we normal order with respect to a Hartree-Fock reference). Both the above errors are non-variational, which can be seen from the DMRG atomic energies as we increase the basis cardinal number; these are -14.66674 (DZ), -14.66669 (TZ), -14.66721 (QZ). For comparison, the best variational calculation for the beryllium atom that we are aware of, using exponentially correlated Gaussian expansions, is -14.66736 beatom . However, both errors also go identically to zero as the orbital basis is increased, because the F12 factor (and the amplitude) is only used to represent the correlation not captured within the basis set.

To obtain more insight into the error from the F12 canonical transcorrelated Hamiltonian, we have computed in Table 2 the F12-CCSD(T)/cc-pCVQZ-F12 binding energies (i.e. CCSD(T) using the canonical transcorrelated Hamiltonian) using the MRCC program of Kállay mrcc , and the conventional CCSD(T)-F12b/cc-pCVQZ-F12 binding energies using the Molpro program package werner2012 . (As pointed out by Knizia et al.ccsdf12b , the CCSD(T)-F12b variant is to be preferred with the large basis sets used here). We observe that the CCSD(T)-F12b/cc-pCVQZ-F12 and F12-CCSD(T)/cc-pCVQZ-F12 binding energies agree very well (to within 5 cm along the entire binding curve). Perfect agreement between the methods is not expected as they correspond to different F12 theories, but these results show that the neglect of three-particle operators in the canonical transcorrelated Hamiltonian produces a description with no significant differences from a standard F12 approach.

To extrapolate the remaining F12 and basis set errors to zero, we carry out a further basis-set completeness (CBS) study. In Table 2, we give the CCSD(T)/aug-cc-pCVnZ binding energies for =4, 5, 6. Following Koputkoput we use the following two basis extrapolation formulae to provide error bars on the complete basis result:

(11) | ||||

(12) |

From Table 2 we observe that the F12-CCSD(T)/cc-pCVQZ-F12 binding energies correspond closely to those of CCSD(T)/aug-cc-pCV6Z. Using Koput’s prescription, we obtain the extrapolated energy as the average of Eqs. (11), (12). At the equilibrium bond length we obtain a basis set limit correction to the DMRG calculation of 87 (19 cm) and an uncertainty of 43 (10 cm). (We estimate the uncertainty as half the extrapolation correction). Thus, the basis-set error remains the largest source of uncertainty in our calculations.

Compared to the experimentally derived well-depths, we find that our directly calculated DMRG (and i-FCIQMC) well-depths, 931.2 cm (924 9 cm), are in excellent agreement with the “experimental” of 929.7 cm (Merritt et al Merritt ) and 934.6 cm (Patkowski et al Patkowski2 ) (Table 4). Including the estimated CBS correction (19 cm), the counterpoise correction (-2 cm), and the relativistic correction (-4 cm), yields a corrected well-depth of 944 cm (DMRG) with an error estimate of 10 cm, which is slightly larger, but still in good agreement with the experimental well-depths. Thus, corrected or otherwise, our calculations compare favorably to the very best experimentally derived well-depths to date. Compared to CCSD(T)-F12, we find that quadruples and higher correlations contribute 25% of the binding energy, indicating significant correlation effects in the ground-state.

The largest absolute discrepancy between our calculations and the experimentally derived curve appears at the shorter bond-length of , where we find the energy (CBS/BSSE/rel. corrected) to be 2.07 m above the equilibrium point as compared to 1.83 m and 2.05 m respectively, in the experimental numbers of Merritt et al. Merritt and Patkowski et al. Patkowski2 . Given the close agreement between our computations and experiment at all other points on the curve (the agreement between the corrected DMRG curve with Patkowski’s curve is better than 0.05 m at all points) the discrepancy with Merritt’s experimental number is quite large. When measured as a multiple of the theoretical uncertainty, we also find that the largest errors are at (2.4) and at (4.0). We note that the inadequacies of Merritt’s fit at longer distances have already been discussed in Ref. Patkowski2 . Our results further suggest that there are inaccuracies in Merritt’s experimental fit at shorter distances as well.

State | DMRG/TZ | DMRG/QZ | DMRG/CBS | MRCI-F12 | MRCI+Q-F12 | EOM-CCSD |
---|---|---|---|---|---|---|

3.61 | 3.59 | 3.57(0.02) | 3.60 | 3.54 | 3.97 | |

3.58 | 3.56 | 3.55(0.01) | 3.70 | 3.55 | 3.48 | |

7.69 | 7.66 | 7.64(0.03) | 8.27 | 8.13 | 7.33 | |

4.81 | 4.78 | 4.77(0.02) | 4.80 | 4.75 | 5.96 |

M | Discarded weight | Energy | Discarded weight | Energy | Discarded weight | Energy | Discarded weight | Energy |
---|---|---|---|---|---|---|---|---|

500 | 7.40 | -29.206566 | 2.43 | -29.207784 | 1.83 | -29.056883 | 6.18 | -29.162740 |

1000 | 1.27 | -29.206794 | 3.40 | -29.207877 | 1.01 | -29.057249 | 1.02 | -29.162905 |

1500 | 4.71 | -29.206827 | 9.48 | -29.207890 | 2.92 | -29.057277 | 2.87 | -29.162931 |

2000 | 2.18 | -29.206836 | 2.67 | -29.207894 | 9.74 | -29.057282 | 1.11 | -29.162939 |

(l.) | -29.206844 | -29.207894 | -29.057280 | -29.162940 | ||||

(q.) | -29.206845 | -29.207895 | -29.057287 | -29.162943 |

We now turn to the excited state DMRG calculations. While accurate ground-state energies can be obtained through composite techniques, this is much more difficult for excited states, due to significantly larger correlation effects. Near exact excited states, however, can be accessed through a state-averaged DMRG calculation Dorando2007 . Combined with the saturated basis set treatment here, the DMRG excitation energies now allow us to present very accurate excitation energies for large basis sets, against which other methods may be compared. The (F12-)DMRG excitation energies, with comparison MRCI-F12, MRCI+Q-F12, and EOM-CCSD energies, are shown in Table 5. The active space used in the MRCI calculations was a 4 electron, 8 orbital complete active space.

The convergence of the DMRG excitation energies with is shown in Table 6 and is plotted in Figure 2. These show that the DMRG energies are converged to within 10 of the formal exact result, and are thus negligible on the eV scale (on the order of tenths of meV’s). The basis set errors for the excitation energies are larger than for the ground-state, because we use the F12 canonical transcorrelated Hamiltonian derived for the ground state to compute all the excitation energies, thus the correlation factor is biased towards the ground-state. To estimate the complete basis set limit of the excitation energy we use Eq. (13) (derived from a fit to CCSD-F12b energies across a large data sethill-extrap ).

(13) |

Since excited state complete basis set extrapolation is less well studied, we estimate the uncertainty conservatively as twice the difference between the estimated complete basis value and the QZ value. As for the ground-states, the basis set error remains the largest uncertainty in the calculations, but even with our conservative estimate ranges only from 0.01 to 0.03 eV.

Overall, MRCI+Q-F12 gives the best agreement with DMRG, with errors of less than 0.05 eV for 3 out of the 4 states. The effect of the Q size-consistency correction is significant, contributing as much as 0.15 eV to the excitation energy. The EOM-CCSD excitation energy errors are large for all states, which is unsurprising given the multireference nature of the ground-state. However, what is most surprising is that for the state, the error of the MRCI+Q excitation energy is as large as 0.4 eV! This indicates extremely strong correlation effects in this state. The state of the beryllium dimer is thus a good benchmark state for the development of excited state methods.

To summarize, in this work we have used explicit correlation via the canonical transcorrelation approach, in conjunction with the density matrix renormalization group and initiator full configuration interaction quantum Monte Carlo methods, to compute the binding curve of the beryllium dimer without the use of composite methods. Our calculations correlate all 8 electrons in basis sets with an orbital basis set of up to 192 basis functions (cc-pCVQZ-F12). Our direct DMRG calculations produce a well-depth of =931.2 cm which agrees very well with the best experimental and theoretical estimates. The remaining basis set effects, BSSE, and relativistic effects, contribute to a final well-depth of =944 cm. We find a significant discrepancy between our computed binding energies and the experimentally derived energies of Merritt et al. at shorter bond-lengths () that suggest inaccuracies in the experimental fits. Finally, using DMRG, we have also computed the excited states at the equilibrium geometry to unprecedented accuracy, highlighting surprisingly strong correlation in the excited states. Overall, we have demonstrated that, by combining explicit correlation with the DMRG or i-FCIQMC methods, it is now possible to directly solve the non-relativistic Schrödinger equation without significant basis set or correlation error for small molecules.

Acknowledgements

This work was supported by National Science Foundation (NSF) through Grant No.
NSF-CHE-1265277. TY was supported in part by the Core Research for Grant-in-Aid for Scientific Research (C) (Grant No. 21550027) from Ministry of Education, Culture, Sports, Science and Technology-Japan (MEXT). CJU would like to acknowledge the NSF grant NSF-CHE-1112097 and the Department of Energy (DOE) grant DE-SC0006650.

## References

- (1) M. Kállay, P. R. Surján, The Journal of Chemical Physics 115, 2945 (2001).
- (2) G. K.-L. Chan, S. Sharma, Annual Review of Physical Chemistry 62, 465 (2011).
- (3) D. Cleland, G. H. Booth, A. Alavi, The Journal of Chemical Physics 132, 41103 (2010).
- (4) G. H. Booth, A. J. W. Thom, A. Alavi, The Journal of Chemical Physics 131, 54106 (2009).
- (5) S. Ten-no, Chemical Physics Letters 398, 56 (2004).
- (6) T. Yanai, T. Shiozaki, The Journal of Chemical Physics 136, 84107 (2012).
- (7) L. Kong, F. A. Bischoff, E. F. Valeev, Chemical Reviews 112, 75 (2011).
- (8) S. Ten-no, Theoretical Chemistry Accounts 131, 1 (2012).
- (9) W. Kutzelnigg, Theoretica Chimica Acta 68, 445 (1985).
- (10) W. Kutzelnigg, W. Klopper, The Journal of Chemical Physics 94, 1985 (1991).
- (11) O. Sinanoglu, The Journal of Chemical Physics 36, 706 (1962).
- (12) M. Torheyden, E. F. Valeev, Physical Chemistry Chemical Physics 10, 3410 (2008).
- (13) G. H. Booth, D. Cleland, A. Alavi, D. P. Tew, The Journal of Chemical Physics 137, 164112 (2012).
- (14) M. Torheyden, E. F. Valeev, The Journal of Chemical Physics 131, 171103 (2009).
- (15) E. Neuscamman, T. Yanai, G. K.-L. Chan, International Review of Physical Chemistry 29, 231 (2010).
- (16) E. Neuscamman, T. Yanai, G. K.-L. Chan, The Journal of Chemical Physics 130, 124102 (2009).
- (17) T. Yanai, G. K.-L. Chan, The Journal of Chemical Physics 124, 194106 (2006).
- (18) T. Yanai, G. K.-L. Chan, The Journal of Chemical Physics 127, 104107 (2007).
- (19) D. Mukherjee, Chemical Physics Letters 274, 561 (1997).
- (20) W. Kutzelnigg, D. Mukherjee, The Journal of Chemical Physics 107, 432 (1997).
- (21) W. Kutzelnigg, D. Mukherjee, The Journal of Chemical Physics 110, 2800 (1999).
- (22) Y. Ohtsuka, S. Ten-no, AIP Conference Proceedings 1456, 97 (2012).
- (23) G. K. L. Chan, M. Head-Gordon, The Journal of Chemical Physics 116, 4462 (2002).
- (24) J. Hachmann, W. Cardoen, G. K. L. Chan, The Journal of Chemical Physics 125, 144101 (2006).
- (25) K. H. Marti, M. Reiher, Physical Chemistry Chemical Physics 13, 6750 (2011).
- (26) O. Legeza, G. Fáth, Physical Review B 53, 14349 (1996).
- (27) G. H. Booth, A. Gruneis, G. Kresse, A. Alavi, Nature 493, 365 (2013).
- (28) F. R. Petruzielo, A. A. Holmes, H. J. Changlani, M. P. Nightingale, C. J. Umrigar, Physical Review Letters 109, 230201 (2012).
- (29) R. Blankenbecler, R. L. Sugar, Physical Review D 27, 1304 (1983).
- (30) N. Trivedi, D. M. Ceperley, Physical Review B 41, 4552 (1990).
- (31) D. F. B. ten Haaf, H. J. M. van Bemmel, J. M. J. van Leeuwen, W. van Saarloos, D. Ceperley, Phys. Rev. B 51, 13039 (1995).
- (32) M. H. Kolodrubetz, J. S. Spencer, B. K. Clark, W. M. C. Foulkes, The Journal of Chemical Physics 138, 024110 (2013).
- (33) D. M. Cleland, G. H. Booth, A. Alavi, The Journal of Chemical Physics 134, 24112 (2011).
- (34) K. Patkowski, R. Podeszwa, K. Szalewicz, The Journal of Physical Chemistry A 111, 12822 (2007).
- (35) J. M. Merritt, V. E. Bondybey, M. C. Heaven, Science 324, 1548 (2009).
- (36) K. Patkowski, V. Špirko, K. Szalewicz, Science 326, 1382 (2009).
- (37) V. E. Bondybey, Science 227, 125 (1985).
- (38) V. E. Bondybey, Chemical Physics Letters 109, 436 (1984).
- (39) J. M. L. Martin, Chemical Physics Letters 303, 399 (1999).
- (40) J. Koput, Physical Chemistry Chemical Physics 13, 20311 (2011).
- (41) R. J. Gdanitz, Chemical Physics Letters 312, 578 (1999).
- (42) M. W. Schmidt, J. Ivanic, K. Ruedenberg, The Journal of Physical Chemistry A 114, 8687 (2010).
- (43) J. G. Hill, K. A. Peterson, Physical Chemistry Chemical Physics 12, 10460 (2010).
- (44) S. Sharma, G. K.-L. Chan, The Journal of Chemical Physics 136, 124121 (2012).
- (45) S. Wouters, P. A. Limacher, D. V. Neck, P. W. Ayers, The Journal of Chemical Physics 136, 134110 (2012).
- (46) I. P. McCulloch, M. Gulacsi, Europhysics Letters 57, 852 (2002).
- (47) G. H. Booth, D. Cleland, A. J. W. Thom, A. Alavi, The Journal of Chemical Physics 135, 84104 (2011).
- (48) G. H. Booth, S. D. Smart, A. Alavi, http://arxiv.org/abs/1305.6981 (2013).
- (49) T. B. Adler, G. Knizia, H.-J. Werner, The Journal of Chemical Physics 127, 221106 (2007).
- (50) J. G. Hill, K. A. Peterson, G. Knizia, H.-J. Werner, The Journal of Chemical Physics 131, 194195 (2009).
- (51) G. Knizia, T. B. Adler, H.-J. Werner, The Journal of Chemical Physics 130, 54104 (2009).
- (52) D. A. Mazziotti, Chemical Reviews 112, 244 (2012).
- (53) D. A. Mazziotti, Chemical Physics Letters 289, 419 (1998).
- (54) D. A. Mazziotti, International Journal of Quantum Chemistry 70, 557 (1998).
- (55) A. E. DePrince, D. A. Mazziotti, The Journal of Chemical Physics 127, 104104 (2007).
- (56) T. Shiozaki, G. Knizia, H.-J. Werner, The Journal of Chemical Physics 134, 34113 (2011).
- (57) H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, Wiley Interdisciplinary Reviews: Computational Molecular Science 2, 242 (2012).
- (58) S. F. Boys, F. Bernardi, Molecular Physics 19, 553 (1970).
- (59) J. Komasa, J. Rychlewski, K. Jankowski, Physical Review A 65, 42507 (2002).
- (60) M. Kállay, P. R. Surján, The Journal of Chemical Physics 115, 2945 (2001).
- (61) J. J. Dorando, J. Hachmann, G. K.-L. Chan, The Journal of Chemical Physics 127, 84109 (2007).
- (62) G. Knizia, T. B. Adler, H.-J. Werner, The Journal of Chemical Physics 130, 054104 (2009).