Explicitly correlated wavefunctions of the ground state and the lowest quintuplet state of the carbon atom

TL;DR

This paper calculates highly accurate nonrelativistic energies for the ground and lowest quintuplet states of the carbon atom using explicitly correlated Gaussian functions, achieving results close to experimental values.

Contribution

It introduces a new approach using symmetry-projected, explicitly correlated Gaussian functions to compute precise energies of carbon atom states, including relativistic corrections.

Findings

Estimated energy limits: -37.844906(4) and -37.691751(2) hartree.

Reproduced experimental excitation energy within about 7 cm$^{-1}$.

Achieved highly accurate nonrelativistic energy calculations.

Abstract

Variational, nonrelativisitic energies have been calculated for the ground state ($^3P_g$) and the lowest quintuplet state ($^5S_u$) of the carbon atom, with wavefunctions expressed in the basis of symmetry-projected, explicitly correlated Gaussian (ECG) lobe functions. New exact limits of these energies have been estimated, amounting to $-37.844906(4)$ and $-37.691751(2)$ hartree. With finite nuclear masses and leading, scalar relativistic corrections included, respective experimental excitation energy of $^{12}C$ has been reproduced with accuracy of about 7 cm$^{-1}$.

Tables

Table 1. Table 1: Nonrelativistic energies, squares of angular momentum and virial ratios for C ∞ superscript 𝐶 {}^{\infty}C . For extrapolated (E) results, standard deviations of the least significant digit are given in parentheses

	${}_{}{}^{3}P_g^{}$			${}_{}{}^{5}S_u^{}$
K	Energy	$⟨L^2⟩$	$-\frac{V}{T}$	Energy	$⟨L^2⟩$	$-\frac{V}{T}$
88	$-37.833555219$	2.00015841	2.000005668	$-37.688655780$	0.00018743	2.000004222
129	$-37.838026763$	2.00009725	1.999990991	$-37.690127107$	0.00008282	2.000009232
189	$-37.840615636$	2.00008198	2.000002689	$-37.690927579$	0.00006823	2.000001843
277	$-37.842526491$	2.00005762	1.999986861	$-37.691351361$	0.00004046	1.999999213
406	$-37.843672608$	2.00003558	1.999993489	$-37.691561153$	0.00002250	2.000000133
595	$-37.844247698$	2.00002376	1.999995273	$-37.691666242$	0.00000907	2.000000081
872	$-37.844533602$	2.00001554	1.999997023	$-37.691707096$	0.00000554	1.999999854
1278	$-37.844711824$	2.00000879	1.999999330	$-37.691727292$	0.00000336	1.999999806
1873	$-37.844794050$	2.00000523	1.999999866	$-37.691739801$	0.00000138	2.000000038
2745	$-37.844851675$	2.00000255	2.000000012	$-37.691745022$	0.00000069	1.999999997
4023	$-37.844877180$	2.00000125	2.000000010	$-37.691747780$	0.00000036	1.999999992
5896	$-37.844889402$	2.00000071	2.000000040
E	$-37.844906\left(4\right)$	2		$-37.691751\left(2\right)$	0

Table 2. Table 2: Comparison of calculated, nonrelativistic energies of C ∞ superscript 𝐶 {}^{\infty}C

method	${}_{}{}^{3}P_g^{}$	${}_{}{}^{5}S_u^{}$
CI-SDTQ ($l_{max}=6$)a	$-37.8393$	$-37.6893$
CI (selected configurations, $l_{max}=4$)b	$-37.8352652$
ECG, K=500c	$-37.84012879$
ECG, K=1000d	$-37.843333$
CCSD(T)-F12e	$-37.844334$
DMCf	$-37.84185\left(5\right)$	$-37.69026\left(3\right)$
DMCg	$-37.84446\left(6\right)$
FC-CFTh	$-37.845004\left(282\right)$	$-37.690998\left(99\right)$
estimated exacti	$-37.8450$
present work (ECG lobes):
variational	$-37.844889402$	$-37.691747780$
extrapolated	$-37.844906\left(4\right)$	$-37.691751\left(2\right)$
aRef. C-Sasaki , bRef. C-Ruiz , cRef. C-Sharkey , dRef. ChemRev , e Ref. Noga-F12 , f Ref. C-Maldonado , g Ref. atoms-Seth , h Ref. FC-Nakatsuji , i Ref. Chakravorty

Table 3. Table 3: Nonrelativistic energies of carbon isotopes (E for extrapolated values); energy differences are in mhartree

K	$E{(}^{12}C)$	$E{(}^{13}C)$	$E{(}^{14}C)$	$E{(}^{13}C)-E{(}^{12}C)$	$E{(}^{14}C)-E{(}^{12}C)$
${}_{}{}^{3}P_g^{}$
4023	$-37.843165183$	$-37.843297308$	$-37.843410138$	$-0.132125$	$-0.244955$
5896	$-37.843177408$	$-37.843309534$	$-37.843422363$	$-0.132126$	$-0.244955$
E	$-37.843194\left(4\right)$	$-37.843326\left(4\right)$	$-37.843439\left(4\right)$
${}_{}{}^{5}S_u^{}$
2745	$-37.690072810$	$-37.690201865$	$-37.690312073$	$-0.129055$	$-0.239263$
4023	$-37.690075568$	$-37.690204623$	$-37.690314831$	$-0.129055$	$-0.239263$
E	-37.690079(2)	-37.690208(2)	-37.690318(2)
Isotopic shift for $D=E{(}^5{S}_u)-E{(}^3{P}_g)$
$D{(}^{13}C)-D{(}^{12}C)=0.003071$ (0.674cm-1), experiment (Ref. experiment ): 0.670(5)cm-1
$D{(}^{14}C)-D{(}^{12}C)=0.005692$ (1.249cm-1)

Table 4. Table 4: Convergence of the components of relativistic shift for C ∞ superscript 𝐶 {}^{\infty}C and comparison with literature data

K	$⟨{\hat{H}}_1⟩$	$⟨{\hat{H}}_2⟩$	$⟨{\hat{H}}_3⟩$	$⟨{\hat{H}}_4⟩$	$⟨{\hat{H}}_{RS}⟩$
${}_{}{}^{3}P_g^{}$
88	$-0.077592107$	$0.062485440$	$0.001154713$	$-0.000016296$	$-0.013968250$
129	$-0.077878128$	$0.062750198$	$0.001139876$	$-0.000015900$	$-0.014003954$
189	$-0.078232316$	$0.063087664$	$0.001130461$	$-0.000015467$	$-0.014029658$
277	$-0.078470993$	$0.063317266$	$0.001127264$	$-0.000015362$	$-0.014041825$
406	$-0.078622237$	$0.063461817$	$0.001121835$	$-0.000015335$	$-0.014053920$
595	$-0.078756303$	$0.063591427$	$0.001118551$	$-0.000015304$	$-0.014061628$
872	$-0.078842674$	$0.063674480$	$0.001115578$	$-0.000015264$	$-0.014067879$
1278	$-0.078927592$	$0.063758450$	$0.001114301$	$-0.000015230$	$-0.014070071$
1873	$-0.078982156$	$0.063811519$	$0.001113287$	$-0.000015222$	$-0.014072572$
2745	$-0.079013248$	$0.063842513$	$0.001112827$	$-0.000015209$	$-0.014073117$
4023	$-0.079034180$	$0.063862568$	$0.001112179$	$-0.000015203$	$-0.014074637$
5896	$-0.079045823$	$0.063874616$	$0.001111739$	$-0.000015200$	$-0.014074667$
HFa	$-0.07799$	$0.06300$	$0.00132$	—	$-0.01367$
MCHFb	$-0.078935$	$0.063934$	$0.001319$	—	$-0.013682$
DMCc	$-0.0790\left(2\right)$	$0.0639\left(2\right)$	$0.001115\left(8\right)$	$-0.000017\left(2\right)$	$-0.0140\left(3\right)$
${}_{}{}^{5}S_u^{}$
88	$-0.075975991$	$0.061299742$	$0.001088726$	$0.000056199$	$-0.013531324$
129	$-0.076395863$	$0.061707535$	$0.001085135$	$0.000056148$	$-0.013547045$
189	$-0.076664784$	$0.061964780$	$0.001082663$	$0.000056189$	$-0.013561152$
277	$-0.076929320$	$0.062222229$	$0.001080459$	$0.000056205$	$-0.013570427$
406	$-0.077035132$	$0.062326668$	$0.001077918$	$0.000056184$	$-0.013574362$
595	$-0.077161830$	$0.062450289$	$0.001077187$	$0.000056167$	$-0.013578169$
872	$-0.077223369$	$0.062511108$	$0.001075924$	$0.000056169$	$-0.013580168$
1278	$-0.077236771$	$0.062524877$	$0.001075463$	$0.000056165$	$-0.013580265$
1873	$-0.077289679$	$0.062575932$	$0.001074897$	$0.000056165$	$-0.013582685$
2745	$-0.077300396$	$0.062587016$	$0.001074547$	$0.000056165$	$-0.013582669$
4023	$-0.077320425$	$0.062606802$	$0.001074403$	$0.000056164$	$-0.013583055$
HFa	$-0.07623$	$0.06171$	$0.00127$	—	$-0.01325$
MCHFb	$-0.077202$	$0.062654$	$0.001278$	—	$-0.013270$
aHartree-Fock Davidson-CH2 , bmulticonfiguration Hartree-Fock C-MCHF   , cQuantum Monte Carlo atoms-Seth

Equations

$$\hat{H}=-\frac{{\nabla }_{nuc}^2}{2m_{nuc}}+\sum _{i=1}^n\left(-\frac{{\nabla }_i^2}{2}-\frac{Z}{r_i}\right)+\sum _{i>j=1}^n\frac{1}{r_{ij}}$$

$$\Psi ({\mathbf{r}}_1,s_1,\dots ,{\mathbf{r}}_n,s_n)=\sum _{I=1}^K{C}_I\hat{A}{\Theta }_I(s_1,\dots ,s_n)\hat{P}{\chi }_I({\mathbf{r}}_1,\dots ,{\mathbf{r}}_n)$$

$${\chi }_I({\mathbf{r}}_1,\dots ,{\mathbf{r}}_n)=\exp \left[-\sum _{i=1}^n{\alpha }_{I,i}({\mathbf{r}}_i-{\mathbf{R}}_{I,i}{)}^2-\sum _{i>j=1}^n{\beta }_{I,ij}{r}_{ij}^2\right]$$

$$\hat{P}{(}^3{P}_g)=\hat{E}+{\hat{C}}_4^1+{\hat{C}}_2+{\hat{C}}_4^3-{\hat{\sigma }}_{v1}-{\hat{\sigma }}_{v2}-{\hat{\sigma }}_{d1}-{\hat{\sigma }}_{d2}$$

$$\Theta (s_1,\dots ,s_6)=[\alpha (1)\beta (2)-\beta (1)\alpha (2)][\alpha (3)\beta (4)-\beta (3)\alpha (4)]\alpha (5)\alpha (6)$$

$$\Theta (s_1,\dots ,s_6)=[\alpha (1)\beta (2)-\beta (1)\alpha (2)]\alpha (3)\alpha (4)\alpha (5)\alpha (6)$$

$$E=E^{(0)}+E^{(2)}+E^{(3)}+\cdots$$

$${\hat{H}}_{RS}={\hat{H}}_1+{\hat{H}}_2+{\hat{H}}_3+{\hat{H}}_4$$

$${\hat{H}}_1=-\frac{1}{8c^2}\sum _{i=1}^n{\nabla }_i^4$$

$${\hat{H}}_2=\frac{Z\pi }{2c^2}\sum _{i=1}^n\delta ({\mathbf{r}}_i)$$

$${\hat{H}}_3=\frac{\pi }{c^2}\sum _{i>j=1}^n\delta ({\mathbf{r}}_{ij})$$

$${\hat{H}}_4=\frac{1}{2c^2}\sum _{i>j=1}^n\left(\frac{{\nabla }_i\cdot {\nabla }_j}{r_{ij}}+\frac{{\mathbf{r}}_{ij}\cdot [({\mathbf{r}}_{ij}\cdot {\nabla }_i){\nabla }_j]}{r_{ij}^3}\right)$$

$$E_{exp}{(}^3{P}_g)=\frac{1}{9}\left[E{(}^3{P}_0)+3E{(}^3{P}_1)+5E{(}^3{P}_2)\right]=29.591c{m}^{-1}.$$

Full text

I Introduction

Carbon has the richest chemistry among all elements. Literally thousands of scientific papers, reporting various calculations on carbon compounds, are written every year. On the other hand, determination of the carbon atom properties on the grounds of theory, with accuracy comparable to that offered by spectroscopic experiments, remains a challenge for computational chemistry. The present work does not provide an ultimate solution of this problem. The aim is more modest – to demonstrate that an explicitly correlated wavefunction is able to yield about 10 $\mu$hartree energy accuracy for six-electron atom and that an ansatz with inaccurate angular dependency of basis functions may prove to be more efficient from that which is seemingly better or at least more elegant.

For a long time, the best variational nonrelativistic energy of the carbon atom ground state was that obtained with the configuration interaction method in the year 1974 C-Sasaki , using Slater orbitals corresponding to the angular momentum quantum number $l$ up to 6, and up to quadruply excited configurations. That result still represents the most accurate published CI energy of this state. Recent CI calculation in which configurations were selected carefully, considering their energy contributions, but built from orbitals the with $l$ limited to 4, gave a little higher energy C-Ruiz . Comparison with the work devoted to boron atom and anion (the latter being isoelectronic with carbon atom’s ground state) B-Almora-Diaz , makes it clear that orbitals with much higher $l$ are needed for building a many-electron basis, capable of yielding the accuracy of about 1 mhartree. CI variational energies have been surpassed by these obtained in calculations with explicitly correlated Gaussian functions C-Sharkey ; ChemRev , but the best result reported is still about 1.5 mhartree above the estimate of exact nonrelativistic energy Chakravorty . This estimate is approached well by nonvariational (or at least not strictly variational) methods – coupled clusters with exponential correlation factor (CC-F12) Noga-F12 , diffusion quantum Monte Carlo simulations C-Maldonado ; atoms-Seth and the “free complement” method FC-Nakatsuji . In the latter, regularized Krylov sequences of functions generated by the system’s Hamiltonian are used as the basis. The results, reported in cited references, differ however even by few milihartree and those, for which standard deviations are given atoms-Seth ; C-Maldonado ; FC-Nakatsuji , do not overlap (table 2 in further text). Calculations on the ${}^{5}S_{u}$ state (the one spin quintuplet below the ionization limit of the carbon atom, known widely for the $sp^{3}$ orbital hybridization model) C-Sasaki ; C-Maldonado ; FC-Nakatsuji yielded somewhat more consistent results, with the discrepancies reaching several hundreds microhartrees.

A comparison of theoretical results with spectroscopic data does not require absolute energies of states. Good agreement with experimental excitation energy of the carbon atom from the ground state to the ${}^{5}S_{u}$ state (and other states too, but they are not the subject of present work) has been achieved in multiconfiguration Hartree-Fock calculations C-MCHF , with omitted correlation of $1s$ core electrons. According to the same reference publication, the leading relativistic energy corrections contribute about 90 cm*-1*. They are partially taken into consideration in the present work.

For few-electron systems, for which high accuracy is desired, an explicitly correlated ansatz, used in the variational framework, is most effective. Unfortunately, the associated computational cost grows rapidly with the number of electrons. Despite of technological progress, both on computer hardware and software side, almost two decades passed between first publications of the ground state energies, computed with nearly-microhartree accuracy, for beryllium Be-KCR and boron B-BA ; B-PKP atoms. Concerning the analytical forms of explicitly correlated basis functions, many of them are tractable for two- and three-electron systems, while only the Hylleraas-CI method and ECGs are competitive in practice for four-electron atoms and atomic ions. Both methods provided energies accurate to few nanohartrees Be-Sims ; Li-Sims ; Be-PPK . The latter ansatz is at present the one applicable for systems containing five and more electrons, due to relatively simple form of integrals appearing in Hamiltonian matrix elements. Other types of basis functions are used too, but occurence of complicated many-electron integrals, without known analytical solutions, forces resorting to stochastic techniques atoms-Seth ; C-Maldonado ; FC-Nakatsuji or using the resolution of identity Noga-F12 for reduction of their complexity.

In the present calculation, the ansatz of explicitly correlated Gaussian lobe functions (called also “Gaussians with shifted centers”, see Ref. MVD for an example) is employed. This ansatz was applied, with a success, in studies of small molecules, molecular ions and van der Waals complexes H2H3+KCR ; HeHeKR1 ; HeHeKR2 ; H3+PA ; LiH-Tung . Free atoms have spherical symmetry, therefore their exact wavefunctions are eigenfunctions of not only Hamiltonian, but also square of angular momentum ($\hat{L}^{2}$) and $z$-component of angular momentum ($\hat{L}_{z}$) operators. Basis functions $\chi$ for atomic states are constructed, as a rule, so that the relations $\hat{L}^{2}\chi=L(L+1)\chi$ and $\hat{L}_{z}\chi=M\chi$ are fulfilled a priori for particular values of $L$ and $M$ quantum numbers B-BA ; B-PKP ; C-Sharkey ; Be-P-Bubin ; Li-D-Bubin ; L3-Sharkey ; N-Sharkey . On the contrary, a lobe function centered off the nucleus is not an eigenfunction of these operators. Convergence towards desired state may however be enforced by variational optimization of trial wavefunction, with proper symmetry constraint. This method, introduced in earlier papers devoted to high $L$ states of the lithium atom Li-KS and various states of many-electron harmonium harm3-JC ; harm4-JC ; harm4W-JC ; harm56-JC ; harm6-KS , will be shortly described in next section. Atomic units are used unless stated otherwise.

II Methods

II.1 Nonrelativistic wavefunction

The stationary Schrödinger equation is solved with the nonrelativistic Hamiltonian of n-electron atom

[TABLE]

The wavefunction, depending on spatial (${\bf r}_{i}$) and spin ($s_{i}$) coordinates,

[TABLE]

with proper permutational symmetry ensured by $\hat{A}$ and primitives $\chi_{I}$ being explicitly correlated Gaussian lobe functions

[TABLE]

is not an eigenfunction of $\hat{L^{2}}$ for non-zero ${\bf R}_{I,i}$ vectors. Deviation of $\langle L^{2}\rangle$ from exact $L(L+1)$ eigenvalue is effectively diminished by the procedure of variational energy minimization, in which nonlinear parameters $\alpha_{I,i}$, $\beta_{I,ij}$ and ${\bf R}_{I,i}$) are optimized. Linear coefficients are determined by solution of the eigenvalue problem, for given set of nonlinear parameters. The convergence towards desired state is ensured by the spatial symmetry projector $\hat{P}$, proper for an irreducible representation of selected finite point group. Action of $\hat{P}$ upon basis functions annihilates a finite subset of their unwanted components whose symmetry properties are specific to others, than the desired one, representations of the $K_{h}$ infinite point group. Particularly, $A_{2}$ representation of the $C_{4v}$ point group was used for the symmetry projector of the ${}^{3}P_{g}$ state

[TABLE]

Confinement of all ${\bf R}_{I,i}$ vectors to the $xy$ plane ensures proper parity (even) of the wavefunction. Lifting this constraint while using the projector proper to the $A_{2g}$ representation of the $D_{4h}$ group offered only negligible energy lowering at substantial increase of computation time even for small basis sets, therefore this alternative path has been abandoned at early stage of the work. $A_{1u}$ representation of the $O_{h}$ point group was employed for the ${}^{5}S_{u}$ state. The symmetry projector is simply too long (48 operations) to be written here explicitly. Identity and all rotation operators, that form the $O$ group, enter this projector with coefficients equal to $1$, and remaining operators (products of the former with the inversion operator) – with coefficients equal to $-1$.

Single spin functions:

[TABLE]

for the triplet, and

[TABLE]

for the quintuplet, are sufficient to ensure convergence to correct variational limits, as the spatial functions are nonorthogonal.

The optimizations of basis set parameters were carried out for infinite nuclear mass. The eigenvalue problem was then solved, in the same basis, for various isotopes of carbon. The center of mass (CM) motion was not separated explicitly from the Hamiltonian, as the wavefunction depends on relative coordinates only (${\bf r}_{i}$ in all equations is the vector of coordinates of $i^{th}$ electron relative to the nucleus). In such case, total kinetic energy operator in the laboratory coordinate frame, acting upon the wavefunction (or a basis function) gives the same result as action of the kinetic energy operator of relative motion, because $\hat{T}=\hat{T}_{CM}+\hat{T}_{rel}$, and $\hat{T}_{CM}\Psi_{rel}=0$ for any definition of internal coordinates BHA-Kutz .

II.2 Relativistic corrections

The relativistic corrections to the energy are obtained in the perturbative series in the fine structure constant $\alpha=\frac{1}{4\pi\epsilon_{0}}\frac{e^{2}}{\hbar c}$. In atomic units, the value of $\alpha$ is equal to the reciprocal of the speed of light in vacuum, $c=137.036$. With the rest mass energy omitted, successive terms in the expansion

[TABLE]

are calculated as expectation values of respective operators, with known nonrelativistic wavefunction. $E^{(0)}$ is the nonrelativistic energy, $E^{(2)}$ contains the Breit-Pauli corrections and higher order terms are the QED (radiative) corrections. The Breit-Pauli Hamiltonian contains the $\hat{H}_{RS}$ operator, which is responsible for the scalar relativistic correction, shifting the energies of whole terms, and the fine and hyperfine structure operators. Only the former is being considered in this work. For fixed nucleus, the relativistic shift Hamiltonian

[TABLE]

consists of following components:

[TABLE]

is the mass-velocity correction,

[TABLE]

is the electron-nucleus Darwin term,

[TABLE]

represents the sum of electron-electron Darwin term and spin-spin Fermi contact interaction (the latter after integration over spin variables Davidson-CH2 ), and

[TABLE]

describes the interaction of magnetic dipoles arising from orbital motion of the electrons. Expectation values of $\hat{H}_{1}$, $\hat{H}_{2}$ and $\hat{H}_{3}$ are known to converge slowly, because these operators sample the wavefunction for short interparticle distances, where ECG functions have an incorrect analytical behaviour (do not describe the wavefunction’s cusps). This deficiency may be overcome by regularization of the problem Drachman , respective technique is however not implemented yet in author’s program, so direct formulas were used in the calculations.

III Results and discussion

The most time-consuming part of the calculations was the optimization of nonlinear parameters of basis functions. It was commenced with very small sets, with sizes of 1, 2 and 3 (stages 1, 2 and 3) respectively. Basis enlarging was performed, with reuse of previously optimized functions in mind. Therefore functions from stage $k-2$ were appended to the set obtained at stage $k$. Consequently, the successive basis sizes formed the Narayana’s cows sequence integers . Each new basis was optimized, function by function, in cycles. Then the expectation value of the $\hat{L}^{2}$ operator and the virial ratio of potential and kinetic energies ($-\frac{V}{T}$) were computed. The results are given in table 1, beginning with 88 ECGs. The convergence of the energy and $\langle L^{2}\rangle$ was substantially better for the ${}^{5}S_{u}$ state, therefore the calculations for this state were finished with 4023 basis functions, while basis of 5896 functions was additionally built for the ground state.

There is no regularity to be found among virial ratios. Their values, close to 2, say only that all parameters were optimized reasonably well. No parameter scaling based on virial ratios was attempted. On the other hand, the squares of angular momentum converge to known exact limits, making it possible to try to extrapolate the energies as functions of $x=\langle L^{2}\rangle-L(L+1)$. There is no theoretical foundation for this extrapolation, other than an observation, that the deviation of $\langle L^{2}\rangle$ is linearly proportional to the rotation energy error, and assumption that the latter is a slowly varying fraction of the total energy error. Extrapolation to $x=0$ with least squares linear regression (fig. 1), using 5 points for each state, yields estimates of exact nonrelativistic energies of the ${}^{\infty}C$ atom. For the ground state, this estimate differs by about 0.1 mhartree from the previous one Chakravorty and is certainly more precise.

The energy convergence with projected ECG lobe functions appears better than with ECGs mutiplied by a proper polynomial of electrons coordinates, which are eigenfunctions of $\hat{L}^{2}$ (table 2). The energies obtained in the present work, for 189 and 406 functions, are lower than those published in refs. C-Sharkey ; ChemRev , computed with 500 and 1000 basis functions respectively. The best published result of Diffusion Quantum Monte Carlo simulation atoms-Seth is surpassed with 872 functions, while final variational energy is lower by about 400 $\mu$hartree. This difference exceeds by far the numerical uncertainty of that simulation. The present result is matched only by the FC-CFT method by Nakatsuji et al. FC-Nakatsuji , at least for the ground state, but with significantly wider error margins.

Nonadiabatic calculations were carried out with nuclear masses $m_{{}^{12}C}=21868.6618$, $m_{{}^{13}C}=23697.6661$, and $m_{{}^{14}C}=25520.3468$, calculated from known molar masses of carbon isotopes (12$u$, 13.003355$u$ and 14.003241$u$ respectively) – dividing them by the Avogadro number and subtracting 6 electron masses. The same basis sets were used as for fixed nucleus, only the linear coefficients in the wavefunctions (Eq. 2) were obtained independently for each isotope. Extrapolations were also based on an assumption that the gap from the best variational energy to the limit, and standard deviation of extrapolated energy, do not change with the nuclear mass. For a given state, the differences of nonrelativistic energies of isotopes converge very quickly and remain stable, therefore only the results obtained with two largest basis sets are given in table 3. The effect of finite nuclear mass contributes $-39.78\mu$hartree (or $-8.731$cm*-1*) to the energy of ${}^{3}P_{g}\rightarrow^{5}S_{u}$ excitation of ${}^{12}C$. The isotopic shift between ${}^{13}C$ and ${}^{12}C$ agrees perfectly with experimental data experiment so a reliable prediction for the ${}^{14}C$ isotope is possible with nonrelativistic wavefunctions (bottom of table 3).

Even with finite nuclear mass taken into account, the nonrelativistic theory is not sufficient to calculate accurate energy differences between atomic states. The ground state term has a fine structure. According to spectroscopic data experiment , the terms of ${}^{12}C$, with $J=1$ and $J=2$, appear respectively at 16.4167(13)cm*-1* and 43.4135(13)cm*-1* above that with $J=0$. The calculation of this split could not be completed in this work, because of lacking implementation of expectation values of spin-orbit and spin-spin coupling operators. At this stage, it is only possible to refer theoretical results to weighted average energy of the ${}^{3}P_{g}$ term:

[TABLE]

The ${}^{5}S_{u}$ term appears at 33735.121(18)cm*-1*, so the reference “excitation energy” amounts to 33705.530cm*-1*. Subtraction of extrapolated, nonrelativistic energies (table 3) gives 33605(1)cm*-1*, which misses the experimental result by about 100cm*-1*.

The convergence of relativistic shifts and their components, for fixed nucleus, is illustrated with the data presented in table 4. As expected, the mass-velocity and electron-nucleus Darwin terms are dominant and their convergence with increasing number of basis functions is unsatisfactory. The differences between the results obtained in two successive largest basis sets still exceed 10$\mu$hartree, for both states. The expectation value of $\hat{H}_{3}$, which contains two-electron Dirac delta operator, converges also slowly but the differences fall below 1 $\mu$hartree. The orbit-orbit interaction energies are converged within few nanohartree. Fortunately the errors of individual components cancel to some extent, owing to the optimization of nonlinear parameters of the wavefunction Cencek-relat , so that five decimal digits of total relativistic scalar corrections for both states seem to be stable and converged even better than their nonrelativistic energies. There is however no perspective to extrapolate these corrections to infinite basis set limit. The results agree fairly well with published relativistic corrections, obtained with less accurate wave functions Davidson-CH2 ; C-MCHF ($\langle\hat{H}_{1}\rangle$, $\langle\hat{H}_{2}\rangle$ and $\langle\hat{H}_{3}\rangle$). The orbit-orbit term has been calculated only for the ground state, in the Quantum Monte Carlo simulation atoms-Seth .

Assuming the same scalar relativistic corrections for the ${}^{12}C$ isotope as for ${}^{\infty}C$ (which is expected to be correct within a fraction of $\mu$hartree for total correction Stanke-FNM ) and adding their values obtained in the largest basis sets to extrapolated nonrelativistic energies from table 3, corrected energies are obtained: $E(^{3}P_{g})=-37.857269$ and $E(^{5}S_{u})=-37.703662$. Their difference amounts to $0.153607$ hartree, or 33713 cm*-1*. It is not possible to calculate its standard deviation, because of lacking error estimation for relativistic corrections. Assuming arbitrarily that the error range is doubled, it would amount to 2 cm*-1*. The missing contribution of at least 6 cm*-1* to the excitation energy might stem from radiative corrections.

Notice should be taken, that the present result is almost equal to that published in Ref. C-MCHF (33711 cm*-1*), but the latter is accurate owing to a fortunate cancellation of errors. The finite nuclear mass effect (about $-9$cm*-1*) was not taken into account there. The orbit-orbit magnetic interaction energy was also omitted, resulting in underestimated relativistic correction to the excitation energy in older calculations Davidson-CH2 ; C-MCHF . This term has the smallest absolute value among all components of scalar Breit-Pauli corrections, but is the one with opposite signs for ${}^{3}P_{g}$ and ${}^{5}S_{u}$ states. It contributes nearly 16 cm*-1* to the excitation energy – almost 15% of the total contribution of relativistic corrections, amounting to 108 cm*-1*.

IV Conclusions

The optimized ECG lobe functions, projected onto proper representations of finite point groups, appear to be a powerful tool for studying the properties of atomic states. Quite surprisingly, they form a more efficient basis, giving lower variational energies at noticeably shorter expansions, than the ECGs with preexponential factors, which are eigenfunctions of $\hat{L}^{2}$. For the first time, the nonrelativistic energies of an 6-electron atom were calculated with accuracy better than 20 $\mu$hartree. Apparent weakness, manifesting oneself in $\langle L^{2}\rangle$ deviating from exact eigenvalue of this operator, may be utilized for energy extrapolation, leading to new estimations of nonrelativistic energies of the lowest triplet and quintuplet states of the carbon atom.

There are ${}^{1}D_{g}$ and ${}^{1}S_{g}$ states of the carbon atom, stemming from the same orbital electron configuration as the ground state ($1s^{2}2s^{2}2p^{2}$), and with energies lower than that of the ${}^{5}S_{u}$ term. They were not studied here, due to excessive computational cost of the optimization of thousands of nonlinear parameters of explicitly correlated basis functions, and they are postponed to a future work. Realization of the present project lasted for about two years and engaged varying number of modern CPU cores, reaching several hundreds working in parallel on calculation of matrix elements.

Concerning the goal to achieve spectroscopic accuracy of quantum-chemical calculations for the carbon atom, there is still a long way to go. The experimental accuracy of the energy difference of the two states considered here, amounting to 0.018cm*-1*, i.e. 82 nanohartree, is by two orders of magnitude better than what the present calculations may offer. The fine structure as well as the contribution of radiative corrections have to be addressed by future work.

V Acknowledgments

This work was financed by the statutory activity subsidy from the Polish Ministry of Science and Higher Education for the Faculty of Chemistry of Wrocław University of Science and Technology (contract number 0401/0121/18). Completion of the project was possible owing to extraordinarily friendly policy of Wrocław Center for Networking and Supercomputing (WCSS, http://wcss.pl), where most calculations have been carried out. Author is grateful for allowing flexible resource allocation, tailored to the needs of the researcher. Thanks are also due to dr Paweł Kędzierski for maintenance of the computer laboratory, where small basis sets were optimized before moving to WCSS.