The Dirac–Coulomb equation with positive-energy projection is solved using explicitly correlated Gaussian functions. The algorithm and computational procedure aims for a parts-per-billion convergence of the energy to provide a starting point for further comparison and further developments in relation with high-resolution atomic and molecular spectroscopy. Besides a detailed discussion of the implementation of the fundamental spinor structure, permutation, and point-group symmetries, various options for the positive-energy projection procedure are presented. The no-pair Dirac–Coulomb energy converged to a parts-per-billion precision is compared with perturbative results for atomic and molecular systems with small nuclear charge numbers. Paper II [D. Ferenc, P. Jeszenszki, and E. Mátyus, J. Chem. Phys. 156, 084110 (2022).] describes the implementation of the Breit interaction in this framework.

A variational solution procedure is reported for the many-particle no-pair Dirac–Coulomb and Dirac–Coulomb–Breit Hamiltonians aiming at a parts-per-billion (ppb) convergence of the atomic and molecular energies, described within the fixed nuclei approximation. The procedure is tested for nuclear charge numbers from Z = 1 (hydrogen) to 28 (iron). Already for the lowest Z values, a significant difference is observed from leading-order Foldy–Woythusen perturbation theory, but the observed deviations are smaller than the estimated self-energy and vacuum polarization corrections.

The Breit interaction is implemented in the no-pair variational Dirac–Coulomb (DC) framework using an explicitly correlated Gaussian basis reported in the previous paper [P. Jeszenszki, D. Ferenc, and E. Mátyus, J. Chem. Phys. 156, 084111 (2022)]. Both a perturbative and a fully variational inclusion of the Breit term are considered. The no-pair DC plus perturbative Breit and the no-pair DC–Breit energies are compared with perturbation theory results including the Breit–Pauli Hamiltonian and leading-order non-radiative quantum electrodynamics corrections for low Z values. Possible reasons for the observed deviations are discussed.

This paper elaborates the integral transformation technique and uses it for the case of the non‐relativistic kinetic and Coulomb potential energy operators, as well as for the relativistic mass‐velocity and Darwin terms. The techniques are tested for the ground electronic state of the helium atom and perturbative relativistic energies are reported for the ground electronic state of the H3+ molecular ion near its equilibrium structure.

Variational and perturbative relativistic energies are computed and compared for two-electron atoms and molecules with low nuclear charge numbers. In general, good agreement of the two approaches is observed. Remaining deviations can be attributed to higher-order relativistic, also called non-radiative quantum electrodynamics (QED), corrections of the perturbative approach that are automatically included in the variational solution of the no-pair Dirac--Coulomb--Breit (DCB) equation to all orders of the $\alpha$ fine-structure constant. The analysis of the polynomial $\alpha$ dependence of the DCB energy makes it possible to determine the leading-order relativistic correction to the non-relativistic energy to high precision without regularization. Contributions from the Breit--Pauli Hamiltonian, for which expectation values converge slowly due the singular terms, are implicitly included in the variational procedure. The $\alpha$ dependence of the no-pair DCB energy shows that the higher-order ($\alpha^4 E_\mathrm{h}$) non-radiative QED correction is 5~\% of the leading-order ($\alpha^3 E_\mathrm{h}$) non-radiative QED correction for $Z=2$ (He), but it is 40~\% already for $Z=4$ (Be$^{2+}$), which indicates that resummation provided by the variational procedure is important already for intermediate nuclear charge numbers.