Drachmann’s regularization approach is implemented for floating explicitly correlated Gaussians (fECGs) and molecular systems. Earlier applications of drachmannized relativistic corrections for molecular systems were hindered due to the unknown analytic matrix elements of 1/rix1/rjy-type operators with fECGs. In the present work, one of the 1/r factors is approximated by a linear combination of Gaussians, which results in calculable integrals. The numerical approach is found to be precise and robust over a range of molecular systems and nuclear configurations, and thus, it opens the route toward an automated evaluation of high-precision relativistic corrections over potential energy surfaces of polyatomic systems. Furthermore, the newly developed integration approach makes it possible to construct the matrix representation of the square of the electronic Hamiltonian relevant for energy lower-bound as well as time-dependent computations of molecular systems with a flexible and high-precision fECG basis representation.
The leading-order (α2Eh) relativistic corrections1 have been computed to high precision for two-, three, and four-electron atoms2–8 and two- and three-electron diatomic molecules.6,9–14 We are aware of a single computation of a polyatomic system () for which the regularized leading-order relativistic corrections were computed (at the equilibrium structure).12
The terms contributing to the leading-order relativistic correction are large in absolute value and have opposite signs, so the correction values must be converged to several digits to arrive at precise final values. Alternative routes to the mainstream perturbation theory approach have been explored and are promising, potentially also for spectroscopic applications,11,15–22 but currently available only for two-electron (or two-spin-1/2-fermion19) systems.
Hence, it is relevant to pursue also the traditional approach1 and develop automatization schemes for the precise evaluation of the involved singular operators over generally applicable (explicitly correlated) Gaussian basis sets. Due to the incorrect description of the (electron–electron and electron–nucleus) coalescence points of the non-relativistic wave function represented as a linear combination of explicitly correlated Gaussians (ECGs), the highly localized, short-range operator expectation values must be regularized to speed up convergence with respect to the basis size.2,6,12 There have been two major directions to correct for this deficiency of general basis sets: (a) using operator identities, which replace local operators with global expressions, often called “drachmannization,” named after the pioneer of this technique,2 and (b) using integral transformation techniques,6,12 which map the short-range behavior to the long-range, for which analytic asymptotic expressions can be derived, and can be used to replace the long-tail part of the integral.
Despite the engaging formal theoretical background12 of the integral transformation technique,6 it is critical to be able to identify a lower threshold value for the long-range part, which is to be represented by the analytic asymptotic expression. Numerical experience shows that a considerable amount of manual fiddling is required to find the optimal threshold value (separating the important physical features represented by numerical integration and the asymptotic limiting function).12,13 This arbitrariness introduces large uncertainties.
For a precise and automated computation of the leading-order relativistic correction value over a multi-dimensional potential energy surface, potentially at hundreds or thousands of nuclear geometries, a robust methodology is required, which can be used in an automated fashion without much manual adjustment.
The drachmannization approach appears to be more robust than the integral transformation technique since it does not contain any parameters to adjust, but the necessary ECG integrals for molecular systems have been thought to be incalculable (in a closed, analytic form).23 In particular, molecular computations with clamped nuclei can be efficiently carried out using the so-called floating explicitly correlated Gaussians (fECGs), in which the shifted (“floating”) centers (treated as variational parameters) are necessary to efficiently account for the significant wave function amplitude off-centered from the origin. The ⟨1/rix1/rjy⟩-type integral expressions (i, j label electron coordinates; x, y label electron or nuclear coordinates), necessary for the evaluation of the drachmannized relativistic corrections, are not known for fECGs, unless the ECG centers are fixed at the origin (zero shift).
Interestingly, the same type of incalculable integrals appear in lower-bound theory, and most importantly its Pollak–Martinazzo24,25 version, which has been demonstrated to be powerful for atomic systems (zero-shift ECGs with analytic integrals),26,27 but no molecular computations have been performed due to the missing integral expressions. Potentially, similar integrals will be required in time-dependent simulations recently proposed with explicitly correlated Gaussians.28
In this work, we develop a numerical approach for the computation of the ⟨1/rix1/rjy⟩-type integrals with fECGs by writing one of the 1/r factors as a linear combination of Gaussians. The representation of 1/r as a linear combination of Gaussians has already been successfully used for deformed explicitly correlated Gaussians by Varga and co-workers,29 using the mathematical scheme of Beylkin and Monzón.30
As a result, a robust numerical drachmannization scheme is developed in this work, and regularized relativistic corrections (computable in a black-box fashion) are reported for two-, three-, and four-electron di- and triatomic molecular systems.
First, we have tested the Gaussian approximation for three special cases, for which analytic integral expressions are also available. Figure 1 highlights the robustness and excellent performance of the numerical approach. As a result of extensive numerical testing, we conclude that high-precision results can be obtained robustly without manual fiddling of the parameters of Eq. (15). For further computations, we have selected the following set of values: M = 200, a = −31, and b = 31, which is a computationally well-affordable choice and allows for sufficiently precise results (close to machine precision in double-precision arithmetic) without further adjustment. An alternative, similarly useful and practical parameterization was M = 400, a = −45, and b = 45, for which the computational cost remains low, and was used to test the (M, a, b) = (200, −31, 31) parameterization. Further refinement of these parameters is possible in the future, but these choices have already provided excellent results with a low computational cost and thereby open the route for the computation of relativistic corrections for polyelectronic and polyatomic molecular systems with unprecedented precision in an automated manner.
![Relative error on a logarithmic scale, log |δ|, of the Gaussian approximation, Eqs. (14) and (15) (with −a = b), tested for special cases for which analytic expressions are also available: (a) Ukl=φk1rijφl; (b) Ukl=φk1rij1rpqφls=0; and (c) Ukl=φk1rij2φl. δ=[Ukl(analytic)−Ukl]/Ukl. Randomly selected pairs of basis functions were used to prepare the figure (deposited in the supplementary material), and they exemplify the general behavior observed for all basis sets used in this study.](/view-large/figure/89101472/211102_1_5.0213079.figures.online.f1.jpg)
Relative error on a logarithmic scale, log |δ|, of the Gaussian approximation, Eqs. (14) and (15) (with −a = b), tested for special cases for which analytic expressions are also available: (a) ; (b) ; and (c) . . Randomly selected pairs of basis functions were used to prepare the figure (deposited in the supplementary material), and they exemplify the general behavior observed for all basis sets used in this study.
Relative error on a logarithmic scale, log |δ|, of the Gaussian approximation, Eqs. (14) and (15) (with −a = b), tested for special cases for which analytic expressions are also available: (a) ; (b) ; and (c) . . Randomly selected pairs of basis functions were used to prepare the figure (deposited in the supplementary material), and they exemplify the general behavior observed for all basis sets used in this study.
Tables I and II present the test results of this numerical drachmannization approach for two-, three-, and four-electron atoms and the two-electron hydrogen molecule, for which benchmark-quality data were already available in the literature. In these examples, we computed the ground state for each system, but the approach is similarly applicable to excited states, too.
Demonstration of the numerical drachmannization approach for the ground state of atoms, He: 11S, Li: 22S, and Be: 21S, and comparison with high-precision reference data from the literature. δEnr = Eref − Enr in Eh. All computations have been carried out with the same (M, a, b) = (200, −31, 31) parameterization of the Gaussian expansion, Eqs. (14) and (15). Convergence tables are provided in the supplementary material.
| −δEnr | ||||
|---|---|---|---|---|
| He | 6 × 10−11 | −13.522 054 | 7.241 717 15 | 0.106 345 364 |
| He3 | 1 × 10−15 | −13.522 016 81 | 7.241 717 274 | 0.106 345 371 2 |
| Li | 3 × 10−10 | −78.557 5 | 41.527 76 | 0.544 322 |
| Li4 a | 2 × 10−13 | −78.556 143 062 | 41.527 828 927 | 0.544 329 79 |
| Be | 5 × 10−7 | −270.706 5 | 141.475 80 | 1.605 301 5 |
| Be7 a | 3 × 10−10 | −270.703 76 | 141.476 010 4 | 1.605 305 33 |
| −δEnr | ||||
|---|---|---|---|---|
| He | 6 × 10−11 | −13.522 054 | 7.241 717 15 | 0.106 345 364 |
| He3 | 1 × 10−15 | −13.522 016 81 | 7.241 717 274 | 0.106 345 371 2 |
| Li | 3 × 10−10 | −78.557 5 | 41.527 76 | 0.544 322 |
| Li4 a | 2 × 10−13 | −78.556 143 062 | 41.527 828 927 | 0.544 329 79 |
| Be | 5 × 10−7 | −270.706 5 | 141.475 80 | 1.605 301 5 |
| Be7 a | 3 × 10−10 | −270.703 76 | 141.476 010 4 | 1.605 305 33 |
Illustration of the enhanced convergence rate by the numerical drachmannization approach developed in this work for the example of the H2 molecule (, Rp–p = 1.40 bohr). δEnr = Eref − Enr in Eh. All computations have been carried out with the same (M, a, b) = (200, −31, 31) parameterization of the Gaussian expansion, Eqs. (14) and (15). Reference 10: Enr = −1.174 475 714 220 443 4(5) Eh.
| −δEnr | ||||
|---|---|---|---|---|
| Direct | ||||
| 1.7 × 10−7 | −1.65 | 0.918 | 0.0169 | −0.205 |
| 1.6 × 10−8 | −1.65 | 0.919 | 0.0168 | −0.205 |
| 1.4 × 10−9 | −1.65 | 0.919 | 0.016 78 | −0.2055 |
| 3.8 × 10−10 | −1.654 | 0.919 | 0.016 77 | −0.2055 |
| Numerical drachmannization | ||||
| 1.7 × 10−7 | −1.655 | 0.919 33 | 0.016 742 | −0.206 |
| 1.6 × 10−8 | −1.6548 | 0.919 335 5 | 0.016 743 1 | −0.2058 |
| 1.4 × 10−9 | −1.654 80 | 0.919 335 7 | 0.016 743 21 | −0.205 75 |
| 3.8 × 10−10 | −1.654 79 | 0.919 335 8 | 0.016 743 24 | −0.205 733 |
| Reference 10 | ||||
| “Direct” | −1.654 837 5(1) | 0.919 37(7) | 0.016 743 5(4) | −0.205 682 (8) |
| “St. reg.” | −1.654 745 (2) | 0.919 336 211 2(18) | 0.016 743 278 1(7) | −0.205 689 (2) |
| “rECG & m.reg.” | −1.654 744 522 5(1) | 0.919 336 211(2) | 0.016 743 278 3(3) | −0.205 688 526(7) |
| −δEnr | ||||
|---|---|---|---|---|
| Direct | ||||
| 1.7 × 10−7 | −1.65 | 0.918 | 0.0169 | −0.205 |
| 1.6 × 10−8 | −1.65 | 0.919 | 0.0168 | −0.205 |
| 1.4 × 10−9 | −1.65 | 0.919 | 0.016 78 | −0.2055 |
| 3.8 × 10−10 | −1.654 | 0.919 | 0.016 77 | −0.2055 |
| Numerical drachmannization | ||||
| 1.7 × 10−7 | −1.655 | 0.919 33 | 0.016 742 | −0.206 |
| 1.6 × 10−8 | −1.6548 | 0.919 335 5 | 0.016 743 1 | −0.2058 |
| 1.4 × 10−9 | −1.654 80 | 0.919 335 7 | 0.016 743 21 | −0.205 75 |
| 3.8 × 10−10 | −1.654 79 | 0.919 335 8 | 0.016 743 24 | −0.205 733 |
| Reference 10 | ||||
| “Direct” | −1.654 837 5(1) | 0.919 37(7) | 0.016 743 5(4) | −0.205 682 (8) |
| “St. reg.” | −1.654 745 (2) | 0.919 336 211 2(18) | 0.016 743 278 1(7) | −0.205 689 (2) |
| “rECG & m.reg.” | −1.654 744 522 5(1) | 0.919 336 211(2) | 0.016 743 278 3(3) | −0.205 688 526(7) |
New numerical results for two-, three-, and four-electron di- and triatomic molecular systems are reported in Table III. These results are the first computation of these leading-order relativistic corrections or improve upon the already available values in the literature. Furthermore, we report the curve of relativistic corrections to the electronic state (Fig. 2) to demonstrate the general applicability of the numerical drachmannization approach along potential energy curves.
Molecular systemsa: numerically drachmannized expectation values and the leading-order relativistic correction. −δEnr is the estimated convergence error, in Eh of the (variational) non-relativistic energy. All converged and 1–2 additional digits are shown. All computations have been carried out with the same (M, a, b) = (200, −31, 31) parameterization of the Gaussian expansion. Convergence tables are provided in the supplementary material.
| −δEnr | ⟨HOO⟩ | |||||
|---|---|---|---|---|---|---|
| 10−12b | −0.624 952 638 6 | 0.318 293 258 6 | ⋯ | ⋯ | −0.124 978 757 0 | |
| b34 c | 10−9 | −0.624 952 7 | 0.318 293 2 | ⋯ | ⋯ | −0.124 978 9 |
| c | 10−9 | −1.933 449 7 | 1.089 654 23 | 0.018 334 666 | −0.057 217 520 | −0.221 442 3 |
| i12 c | 10−9 | −1.933 424 | 1.089 655 | 0.018 335 | −0.057 218 | −0.221 422 |
| HeH+ | 10−9 | −13.419 497 | 7.216 244 8 | 0.101 121 71 | −0.141 241 61 | −1.907 804 |
| i15 c | 10−9 | −13.419 29 | 7.216 25 | 0.101 122 | −0.141 242 | −1.907 58 |
| 10−8 | −23.867 58 | 12.537 54 | 0.117 596 4 | −0.086 577 | −3.890 79 | |
| i13 c | 10−8 | −23.866 14 | 12.537 51 | 0.117 63 | −0.086 58 | −3.889 30 |
| H3 | 10−8 | −2.277 870 | 1.236 566 | 0.016 712 | −0.046 780 1 | −0.329 754 2 |
| He2 | 10−6 | −24.157 271 | 12.660 991 | 0.119 967 0 | −0.073 154 | −3.965 700 |
| −δEnr | ⟨HOO⟩ | |||||
|---|---|---|---|---|---|---|
| 10−12b | −0.624 952 638 6 | 0.318 293 258 6 | ⋯ | ⋯ | −0.124 978 757 0 | |
| b34 c | 10−9 | −0.624 952 7 | 0.318 293 2 | ⋯ | ⋯ | −0.124 978 9 |
| c | 10−9 | −1.933 449 7 | 1.089 654 23 | 0.018 334 666 | −0.057 217 520 | −0.221 442 3 |
| i12 c | 10−9 | −1.933 424 | 1.089 655 | 0.018 335 | −0.057 218 | −0.221 422 |
| HeH+ | 10−9 | −13.419 497 | 7.216 244 8 | 0.101 121 71 | −0.141 241 61 | −1.907 804 |
| i15 c | 10−9 | −13.419 29 | 7.216 25 | 0.101 122 | −0.141 242 | −1.907 58 |
| 10−8 | −23.867 58 | 12.537 54 | 0.117 596 4 | −0.086 577 | −3.890 79 | |
| i13 c | 10−8 | −23.866 14 | 12.537 51 | 0.117 63 | −0.086 58 | −3.889 30 |
| H3 | 10−8 | −2.277 870 | 1.236 566 | 0.016 712 | −0.046 780 1 | −0.329 754 2 |
| He2 | 10−6 | −24.157 271 | 12.660 991 | 0.119 967 0 | −0.073 154 | −3.965 700 |
The electronic state and molecular geometry was , : Rp–p = 12.00 bohr; H2, : Rp–p = 1.40 bohr; , 1: Rp–p = 1.65 bohr, ; H3, 12Σ+: bohr, bohr, ; HeH+, 11Σ+: Rp–α = 1.46 bohr; , : Rα–α = 2.042 bohr; He2, : Rα–α = 2.00 bohr.
The non-relativistic energy has been reported with higher precision in Ref. 36.
Regularized relativistic corrections to the electronic state near the local minimum (near Rpp = 12.5 bohr) computed by the automated numerical drachmannization approach proposed in this work. The computed values confirm (and are more precise) than the ones reported by Kennedy and Howells using special basis functions.34
Regularized relativistic corrections to the electronic state near the local minimum (near Rpp = 12.5 bohr) computed by the automated numerical drachmannization approach proposed in this work. The computed values confirm (and are more precise) than the ones reported by Kennedy and Howells using special basis functions.34
This work reports a robust and generally applicable numerical drachmannization approach for the computation of the leading-order relativistic corrections with floating explicitly correlated Gaussian basis functions. Formerly incalculable ⟨1/rix1/rjy⟩-type floating ECG integrals are made calculable using a simple and robust Gaussian approximation to one of the 1/r factors. The computational approach is demonstrated to be numerically robust and generally applicable for two-, three-, and four-electron atomic, diatomic, and polyatomic systems and along potential energy curves (and potentially, surfaces). The numerical drachmannization approach developed in this work opens the route to the automated computation of high-precision relativistic corrections to potential energy curves and surfaces in the future.
SUPPLEMENTARY MATERIAL
The supplementary material includes (1) details of the integral calculation and implementation; (2) convergence tables; and (3) basis set parameterizations.
The ACKNOWLEDGMENTS
Financial support of the European Research Council through a Starting Grant (No. 851421) is gratefully acknowledged. We thank Péter Jeszenszki and Robbie Ireland for discussions and joint work over the past years.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Balázs Rácsai: Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Dávid Ferenc: Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Ádám Margócsy: Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Edit Mátyus: Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available within the article and its supplementary material.
