PHYSICAL REVIEW A 101, 062508 (2020)

Solving the Schrödinger equation of atoms and molecules using one- and two-electron integrals only
Hiroshi Nakatsuji ,* Hiroyuki Nakashima , and Yusaku I. Kurokawa
Quantum Chemistry Research Institute, Kyoto Technoscience Center 16 14 Yoshida Kawaramachi, Sakyo-ku, Kyoto 606-8305, Japan(Received 21 February 2020; revised manuscript received 13 April 2020; accepted 1 May 2020;
published 4 June 2020)
A variational theory called free-complement (FC) si j theory for solving the Schrödinger equation of atoms
and molecules using only one- and two-electron integrals over Slater or Gaussian functions is proposed. It isderived from the scaled Schrödinger equation [Phys Rev. Lett. 93, 030403 (2004)] by replacing the two-electron
part ri j of the scaling g operator with si j = ri2j , both of which solve the divergence difficulty inherent to theoriginal Schrödinger equation by avoiding the interelectron collisions. The si j function can be rewritten with only
one-electron functions so that when the initial wave function of the FC theory is composed of only one-electron
functions, the FC wave function including si j can be rewritten with only one-electron functions. Therefore,
the variational calculations of the FC si j theory can be performed with only one- and two-electron integrals.
However, in comparison with the ri j function, the si j function behaves less efficiently when two electrons comeclose to each other: the electron-electron cusp condition is not satisfied with the FC si j theory, though the wave
function has explicit ri j dependence. On the other hand, the electron-nuclear cusp condition is satisfied, even with
the Gaussian functions, for the presence of the electron-nuclear function riA in the g operator. Test applicationsof the FC si j theory were done to He, Li, and the 5 S o sp3 state of carbon atom and to hydrogen molecule using
the local-molecular orbital (MO)- and valence bond (VB)-type initial functions. We examined both Slater andGaussian functions. As the order of the FC theory increased, the wave functions were improved and the energies
approached from a few kcal/mol to less than 1 kcal/mol accuracies relative to the known exact values. Thus,
with the FC si j theory, the Schrödinger equation was solved to less than 1 kcal/mol accuracy for all the systems
examined. The costs for the si j calculations were much smaller than those of the ri j theory. However, for the 5 S o
sp3 state of the carbon atom, we observed that the convergence rate became slow when the calculation came closeto ∼ 1 kcal/mol accuracy. This suggests that for the FC si j theory, the calculations would become easier if therequired accuracy is within a few kcal/mol for the absolute energy. Actually, what we need in chemical studies
are mostly related to the difference energy, whose accuracy could be better than that of the absolute energy withthe robust variational method. Two lines of future studies were suggested.
DOI: 10.1103/PhysRevA.101.062508

the Hamiltonian of Eq. (2) in the variational calculations [3].
To overcome this difficulty, one of the authors introduced in
2004 the scaled SE (SSE) [3] given by

I. INTRODUCTION

Since the Schrödinger equation (SE) given by(H − E )ψ = 0,

(1)

 1
 ZA 
1
i −
A −
+
H =−
i 2
A 2
i,A riA
i< j ri j

ZA ZB
+
,
(2)
A<B RAB 1

is a governing principle of chemical science [1], a goal or astarting point of quantum chemistry is to establish a general
method of solving the SE. With this aim, starting from thestudy of the structure of the exact wave function [2], ourresearch has been focused on finding the method of solving
the SE as accurately and as easily as possible [3–9]. A reason
why the SE was left unsolved for about 80 years since its
discovery lay in the divergence difficulty of the integrals of*

Author to whomh.nakatsuji@qcri.or.jp

correspondence

2469-9926/2020/101(6)/062508(10)shouldbe

addressed:g(H − E )ψ = 0,
where the scaling operator g was defined by


g=
riA +
ri j .
i,A

(3)

(4)

i< j

This operator is always positive except for the colliding points
r0 , and even at r0 , this operator does not eliminate the information of the Hamiltonian since
lim gV = 0,

r→r0

(5)

where V is the electron-nuclear attraction and electronelectron repulsion potentials in the Hamiltonian given by
Eq. (2); we assume here Born-Oppenheimer (BO) approximation, though our theory is applicable also to the non-BO
case [10]. For these natures of the g operator, the SE and theSSE has the same set of solutions: we can divide SSE by gwithout affecting the solutions. So, we solve the SSE insteadof the SE. Then, we do not meet the divergence difficulty of

062508-1©2020 American Physical Society

NAKATSUJI, NAKASHIMA, AND KUROKAWA

PHYSICAL REVIEW A 101, 062508 (2020)

the integrals. The general method of solving the SE and SSE
was formulated in several variants [4]. Among them, a simple
way of solving the SSE is to use the simplest ICI (iterativecomplement interaction) formula [3,4,6] given byψn+1 = [1 + Cn g(H − En )]ψn ,(6)which leads to the exact solutions of the SE (or SSE) atconvergence.
A more efficient method, originally called free ICI [3] andrenamed later as the free-complement (FC) method [6], wasformulated from Eq. (6). Starting from some initial function
ψ0 , we apply the simplest ICI method as given by Eq. (6) tothe order n + 1. Then, from the right-hand side of Eq. (6),we select all the independent analytical functions, throw out
all the diverging functions, collect the rest that are referred
to as complement functions (cf’s) {φI }, give independent
coefficients {cI } to all of them, and we obtain the expressionof the FC wave function written asψn+1 =M
cI φI .

(7)

I=0

The cf φI is the function of riA and ri j and is written by
φI = fI (riA , ri j )φ0 ,(8)with φ0 = ψ0 and f0 = 1. This FC wave function becomes
potentially exact as n increases, and the structure of the FC
wave function is related to the theory published in 1985
[11]. Thus, when we calculate the coefficients {cI } using thevariational principle, we can get the essentially exact solutionsof the SE. For example, for the helium atom, we could obtain
the energy that was accurate to over 40 digits [12] and theessentially exact nature of the solutions was demonstrated
from the flatness of the local energy plane Hψ (r)/ψ (r)[13]. For the hydrogen molecule, the FC theory gave highly
accurate results [3,14]. Recently, we have reported essentiallyexact potential curves for the ground and many excited statesof different symmetries and reported the vibrational levelsof many different electronic states [15]. Recently, we could
obtain highly accurate wave functions for the ground andmany excited states of Li and Be atoms [16].
However, in the variational calculations of larger
molecules, the integral evaluations become difficult forthe cf’s including ri j functions: they cause three-, four-,and even more-electron integrals. When we use the local
Schrödinger equation (LSE) method [5], instead of thevariational principle, we become free from the integral
evaluations and we could obtain highly accurate solutions of
the SE using the sampling-type methodology for the first-row
atoms and small molecules [5–7,9]. It was also shown thatthe effect of the Pauli principle decays exponentially as the
distance of the electron pair increases and using this fact,
the interexchange (iExg) theory was proposed [7]. Combinedwith the iExg theory, the FC theory for solving the SE wasapplied to large systems like He fullerene [7], and further
simplifications were formulated for the FC theory to be
applied to large and even giant molecular systems [8].
In this paper, we consider a different method of solving
the SE, using the variational principle: the robust nature of
the variational method is very attractive. The above method

includes the two-electron ri j terms in the wave function. Thisis an origin of the high accuracy of the theory, but at thesame time, it caused a difficulty in the variational method in
the evaluation and handling of the three-, four-, etc., electron
integrals. In contrast, popular quantum-chemical programs
handle usually only one- and two-electron integrals, though
they usually do not handle highly accurate wave functions
close to the Schrödinger accuracy. In this paper, we propose
a variational theory for solving the SE accurately by using
one- and two-electron integrals only. Test applications of the
proposed method will be performed for the He and Li atoms,
the 5 S (sp3 ) state of the carbon atom, and the hydrogenmolecule. Then, some concluding remarks will be given.
II. FREE-COMPLEMENT si j THEORY

When we use the variational method for solving the FC
theory explained above, we have to evaluate the integrals
originating from the g operator given by Eq. (4). The first riAterm caused no problem since it is a one-electron function,
but the second ri j term did, as explained above, because it is a
two-electron function. To relax this problem, we here consider
the possibility of using the g operator given by


g=
riA +
si j , si j = ri2j .
(9)
i,A

i< j

The first one-electron term is the same, but the second twoelectron term is different. Like the ri j function, the si j function
can eliminate the divergence difficulty caused by the electronelectron repulsion potential in the Hamiltonian. However, thesi j function would be less efficient than the ri j function, since
si j is smaller than ri j when two electrons i and j are in the“colliding” region less than ∼1 au. Further, the product of thesi j term and the potential 1/ri j becomes zero at the colliding
position ri j = 0, destroying the relation given by Eq. (5).
Nevertheless, the reason of introducing the si j function insteadof ri j lies in a valuable merit that si j is transformed into the
one-electron functions alone assi j = ri2j = xi2j + yi2j + zi2j
= (xi − x j )2 + (yi − y j )2 + (zi − z j )2
= ri2 + r 2j − 2xi x j − 2yi y j − 2zi z j ,(10)where the last two lines include only one-electron functions.
Our initial function ψ0 is composed of only one-electron
functions. So, when we use the si j function, instead of theri j function, all the cf’s {φI } of the FC theory are composedof only one-electron functions. Then, the Hamiltonian andoverlap integrals between the cf’s necessary for the variational
calculations are all written within the one- and two-electron
integrals [the two-electron integrals come from the interelectron potential 1/ri j in the Hamiltonian given by Eq. (2)]. We
refer to this theory as the FC si j theory and the original one
as the FC ri j theory. Like Eq. (8), the cf’s of the FC si j theory
are written as
φI = fI (riA , si j )φ0 .

(11)

Now, let us examine the FC si j theory from a general theoretical standpoint of solving the SE. As explained above,the si j operator is expected to be less efficient than the ri j

062508-2

SOLVING THE SCHRÖDINGER EQUATION OF ATOMS …

PHYSICAL REVIEW A 101, 062508 (2020)

operator. In particular, it does not satisfy Eq. (5) for theinterelectron repulsion potential Ve−e , so that the si j operator
erases the information of the SE at re−e = 0. This means thatthe electron-electron cusp condition given by∂ ψ̄ = λψ (ri j = 0),∂ri j ri j =0(12)where ψ̄ on the left-hand side is the spherically averaged wave
function around ri j = 0, is not satisfied for the wave function
calculated by the FC si j theory. The value of λ should be
1/2 and 1/4, respectively, for the singlet and triplet electronelectron pairs, but it is calculated to be 0 with the FC si j theory.
However, we note that the wave function of the FC si j theory
includes the explicit dependence of ri j in the form of si j , sothat the value of ∂ψ/∂ri j is not necessarily zero except at theelectron-electron collisional points. On the other hand, when
the wave function does not include the explicit ri j function,
the value of ∂ψ/∂ri j is zero everywhere.
The FC si j theory satisfies the electron-nuclear cusp condition which is important since the electron-nuclear attraction
potential becomes attractively infinite at the cusp region. The
riA function in the g operator of Eq. (9) helps to satisfy
the electron-nuclear cusp condition even when we use theGaussian functions.
On the other hand, the electron-electron cusp condition
should be less important than the electron-nuclear cusp condition since at the cusp region, the electron-electron repulsion
energy becomes infinitely large so that the probability of theoccurrence of the electron-electron collision must be small.
Particularly, for the same spin pair, the probability should bezero from the Pauli principle. Nevertheless, the e − e cusp
condition should be important for the theory close to the exact
level.
Different from the ri j function, the si j function is written
in some different expressions of the one-electron functions
like Eq. (10). In the variational si j calculations, we can utilize
this fact by relaxing the different elements. For example, we
may handle the ri2 + r 2j part and the xi x j + yi y j + zi z j part
independently, and in molecules, we may handle the x, y, and z
parts independently. With such relaxations, we can get a lower
energy than otherwise. In the applications below, we adopted
this relaxation method.
Since the computational merit of the FC si j theory that only
one- and two-electron integrals are necessary in the variational
calculations is so attractive over the demerit, we will examinein this paper how well this theory solves the SE for some small
systems.
III. FC si j CALCULATIONS WITH SLATERAND GAUSSIAN FUNCTIONS

We perform the FC si j calculations using both Slater andGaussian functions. Efficient integral evaluation methods forthese functions are available in the literature [17–19]. For
atoms, we used our own program for integral evaluations[20] and for the hydrogen molecules, we used the SMILES
program [21] for Slater functions and our own program [22]for Gaussian functions.

The cf’s {φI } of the FC wave function are produced
from the simplest ICI wave function given by Eq. (6) by
the procedure written under Eq. (6). The change caused by
the replacement of ri j with si j is only that when we haveri j with this procedure, we have to replace it with si j . As
shown before in the applications of the FC theory to thehydrogen atom using both Slater and Gaussian functions [23],the Gaussian functions produce a larger number of cf’s thanthe Slater functions, mainly through the kinetic operator in
the Hamiltonian. For the hydrogen atom, the number of the
cf’s was three times larger with the Gaussians than with
the Slaters, but the results were similar when compared in thesame order; the results were slightly better with the Gaussians
than with the Slaters [23].
When we use the initial function composed of the Gaussianfunction, xil yim zin exp(−βri2 ), the FC theory introduces theunusual one-electron Gaussian functions ri xil yim zin exp(−βri2 ),which we call the r-Gauss function, including the ri function
originating from the g operator: we need to consider only
ri instead of rik since in the cases of k = 2n + 1, the usageof ri2 = xi2 + yi2 + zi2 transforms them into the above simple
r-Gauss formula. In the FC formulation, the presence ofthe r-Gauss function is essential to realize the results of
the Schrödinger accuracy: for example, the electron-nuclear
cusp condition is satisfied due to the presence of the r-Gauss
function as shown previously [23]. The one-center integrals
including r-Gauss functions are easily calculated, but mostof the multicenter integrals are not known. Therefore, forthe multicenter integrals including the r-Gauss function, weused the Gaussian-expansion method [24] due originally to Oohata, Taketa, and Huzinaga [25]. Since the electron-nuclear
cusp condition is mostly related to the one-center atomic integrals, the values of the electron-nuclear cusps in molecules are
scarcely affected by the accuracy of this Gaussian expansion
method.
From the expression of the si j operator given by Eq. (10),we understand that the FC si j calculations will require the cf’s
with higher angular momentum functions. For the calculations of atoms, this does not cause any problem. However,for molecules, this causes a problem since most molecularquantum-chemical programs available have limits on the highest angular quantum numbers. This is true for both Gaussianand Slater programs. Therefore, for the present FC si j calculations, we have to improve the codes so as to adapt to higher
angular momentum integrals [22]. The theoretical efforts on
this subject are found in the literature for both Slater [26] andGaussian [27] functions.
Now, let us apply the FC si j theory to the helium atom. We
used both Slater and Gaussian functions. As the initial functions, we used two doubly occupied 1s orbitals of different
exponents: 1s(k) = exp(−α1s,k r) and 1s(k) = exp(−α1s,k r 2 )(k = 1,2) for Slaters and Gaussians, respectively. We generated the cf’s of the FC theory using the procedure described
below of Eq. (6). At each order of the FC theory, we optimized only the two exponents α1s,1 and α1s,2 included in
φ0 commonly to all the cf’s involved, which we referred toas “double-zeta” optimization. The results are summarized inTable I. Initially at n = 0, the Gaussian functions were muchworse than the Slater functions, since it was composed of only
two primitive Gaussian functions. As the order n increased,

062508-3

NAKATSUJI, NAKASHIMA, AND KUROKAWA

PHYSICAL REVIEW A 101, 062508 (2020)

TABLE I. Free-complement si j theory applied to He atom using optimized double-zeta Slater and Gaussian functions.a
Slater functionnM

Energy (a.u.)

0
1
2
3
4

2
6
14
26
44

4

8
9
15
17
18
21
23
28

−2.847 656−2.886 147−2.901 524−2.902 698−2.903 167(cf selection)−2.899 141−2.900 571−2.902 342−2.902 568−2.902 626−2.902 802−2.902 950−2.903 088−2.903 724b

Exact

Gaussian function

E = EFC − Eexact (kcal/mol)
35.18311.030
1.381
0.644
0.349

n

M

0
1
2
3

2
14
64
164

3
2.876 (<3)
1.979 (<2)
0.867 (<1)
0.725 (<0.8)
0.689 (<0.7)
0.578 (<0.6)
0.486 (<0.5)
0.399 (<0.4)

Energy (a.u.)−2.473 516−2.886 477−2.902 432−2.903 234(cf selection)
23−2.899 100
27−2.900 712
38−2.902 148
42−2.902 311
45−2.902 451
47−2.902 655
52−2.902 795
60−2.902 951
73−2.903 096−2.903 724b

E = EFC − Eexact (kcal/mol)
269.959
10.8230.811
0.307
2.901 (<3)
1.890 (<2)
0.989 (<1)
0.887 (<0.9)
0.799 (<0.8)
0.671 (<0.7)
0.583 (<0.6)
0.485 (<0.5)
0.394 (<0.4)

Initial functions with the double-zeta 1s orbitals: φ0(k) = (1s(k) )2 αβ (k = 1, 2). The Slater orbitals used were 1s(k) = exp(−α1s,k r) withthe optimized sets of (α1s,1 , α1s,2 ) = (1.688, 1.688), (1.947, 3.365), (1.818, 3.546), (2.203, 4.262), (2.508, 5.012) for n = 0 to 4: at n = 0,
the optimized values of α1s,1 and α1s,2 became equal to the value of α1s of the single-zeta case. The Gaussian orbitals used were 1s(k) =exp(−α1s,k r 2 ) using the optimized sets of (α1s,1 , α1s,2 ) = (0.552, 2.087), (0.417, 0.904), (0.328, 2.713), (0.326, 3.779) for n = 0 − 3.
bExact energy from Ref. [12].
athe number of the cf’s, M increases, and the energies for both
Slater and Gaussian functions were improved and quickly approached the exact value [12]. The differences from the exact
energy E became smaller than 1 kcal/mol at n = 3 andn = 2 with the Slater and Gaussian functions, respectively.
Note that at this accurate level, the Gaussian functions worked
better than the Slater functions when compared in the order of
the FC calculations.
After obtaining the chemical accuracy, we performed the
cf-selection method: we examined the energetic contributionof each cf included in the wave function, and ordered themin the energetic contribution. With this method, we can findthe smallest set of the cf’s necessary for the desired accuracy. Then, we can reproduce the chemical accuracy withthe smaller number of the cf’s. This cf-selection method is
explained in the Appendix.
We applied the cf-selection method to the n = 4 Slater
results and to the n = 3 Gaussian results. It was shown thatfor obtaining the accuracy of 3 kcal/mol relative to the exact
energy, 8 Slater and 23 Gaussian cf’s were necessary, forthe 2 kcal/mol accuracy, 9 Slater, and 27 Gaussian cf’s, forthe 1 kcal/mol accuracy (chemical accuracy), 15 Slater, and
38 Gaussian cf’s, for the 0.8 kcal/mol accuracy, 17 Slater,and 45 Gaussian cf’s, and so forth. Generally, the Gaussians
requires two to three times more cf’s than the Slaters. We cancompare the present results of the si j theory to those of the ri jtheory [12] and we find that the si j theory is less efficient than
the ri j theory, though the calculations are easier with the si jtheory.
We next applied the FC si j theory to the Li atom using both
Slater and Gaussian functions and the results were summarized in Table II. We used two sets of (1s)2 (2s)1 functions with
1s = exp(−α1s r) and 2s = r exp(−α2s r) − d exp(−α2s r) for
the Slater case having different sets of exponents and coeffi-

cient d as given in the footnote of the table. For Gaussians,
the functions were of the type exp(−αr 2 ) but others were
essentially the same. As seen from Table II, the difference
from the exact energy E approaches rapidly from a few
kcal/mol to less than 1 kcal/mol; the chemical accuracy was
obtained at n = 4 with the Slater functions and at n = 3 with
the Gaussian functions. Again, the order n necessary to obtainthe chemical accuracy was smaller with the Gaussians thanwith the Slaters. When we apply the cf-selection method after
obtaining the chemical accuracy, we see that the chemical
accuracy was realized with the selected 37 cf’s for the Slatercase and with the selected 164 cf’s for the Gaussian case,which are much smaller than the total numbers of the cf’s at
the orders 4 and 3, respectively, for the Slater and Gaussian
cases, respectively.
Thus, with the FC si j theory, we could obtain chemical
accuracy for both Slater and Gaussian cases. From the ordern, Gaussian was more efficient than Slater as we observed
earlier for the hydrogen atom [23]. However, the number ofthe cf’s is larger with the Gaussians than with the Slaters,
since with the former, the cf’s are generated not only fromthe potential operators in the Hamiltonian, but also from thekinetic operators, in the procedures described from Eqs. (6)and (7). When we compare the results of the FC si j theory
with those of the FC ri j theory, we see that the former is lessefficient than the latter: the result of the ri j method can be seen
for Li from our recent paper [16] from the ground-state line of
Table VIII. The FC ri j theory gave more compact descriptions
than the FC si j theory.
We next apply the FC si j theory to solving the SE of the
5 o
S (sp3 ) state of the carbon atom. This state is not only themost important state of the carbon atom, which is the mostimportant element in chemistry and biology, but also a hot
subject in the variational computational quantum chemistry.

062508-4

SOLVING THE SCHRÖDINGER EQUATION OF ATOMS …

PHYSICAL REVIEW A 101, 062508 (2020)

TABLE II. Free-complement si j theory applied to Li atom using optimized double-zeta Slater and Gaussian functions.a
Slater function
nb

M

0
1
2
3
4

2
10
36
100
222

4

Exact

Energy (a.u.)−7.419 190−7.458 890−7.474 094−7.476 270−7.476 932(cf selection)
17−7.473 381
22−7.475 025
37−7.476 489
43−7.476 628
64−7.476 791

Gaussian function

E = EFC − Eexact (kcal/mol)
36.94112.029
2.489
1.123
0.707

n

M

0
1
2
3

2
34
322
1666

3
2.936 (<3)
1.904 (<2)
0.986 (<1)
0.898 (<0.9)
0.797 (<0.8)−7.478 060b

Energy (a.u.)−6.586 026−7.444 839−7.475 776−7.477 270(cf selection)
73−7.473 416
102−7.474 909
164−7.476 472
181−7.476 627
202−7.476 786
231−7.476 949
312−7.477 105872−7.477 264−7.478 060b

E = EFC − Eexact (kcal/mol)
559.760
20.8461.433
0.496
2.914 (<3)
1.977 (<2)
0.997 (<1)
0.899 (<0.9)
0.799 (<0.8)
0.697 (<0.7)
0.599 (<0.6)
0.499 (<0.5)

Initial functions using double-zeta 1s and 2s orbitals: φ0(1) = (1s(k) )2 αβ(2s(k) )α (k = 1, 2). The Slater orbitals usedwere 1s(k) = exp(−α1s,k r), 2s(k) = r exp(−α2s,k r) − dk exp(−α2s,k r) with the optimized sets (α1s,1 , α2s,1 , d1 , α1s,2 , α2s,2 , d2 ) =(2.689, 0.677, 0.4605, 2.700, 0.865, 0.2730), (3.038, 0.696, 0.5011, 5.179, 0.690, 0.4211), (2.897, 0.746, 0.4747, 5.541, 0.671, 0.2266),(3.023, 0.843, −6.8742, 6.330, 0.681, 0.2842), (3.390, 0.895, −5.4772, 7.656, 0.686, 0.4515) for n = 0 − 4. The Gaussian orbitals used were1s(1) = exp(−α1s,1 r 2 ), 2s(1) = r exp(−α2s,1 r 2 ) − d1 exp(−α2s,1 r 2 ), 1s(2) = exp(−α1s,2 r 2 ), 2s(2) = r exp(−α2s,2 r 2 ) − d2 exp(−α2s,2 r 2 ) with the
optimized sets (α1s,1 , α2s,1 , d1 , α1s,2 , α2s,2 , d2 ) = (1.445, 0.0628, −4.9078, 5.457, 0.0654, −0.3614), (1.072, 0.0643, −3.7342, 2.536, 0.0677,−0.4632), (0.826, 0.0633, −10.768, 3.967, 0.0798, 0.04690), (0.802, 0.0635, −13.2310, 5.123, 0.0590, −1.1715) for n = 0 − 3.
bExact energy from Ref. [16].
aTo the best of our knowledge, Noga et al. [28] used the
CCSD(T)-F12 method and Brown et al. [29] used the quantum
Monte Carlo method for the 3 P ground state and we used the
FC-LSE method for the 3 P and 5 S o states of the carbon atom[9], all methods being nonvariational. Recently, Strasburger[30] performed the variational calculations for the 3 P and5 o
S states with the symmetry-projected, explicitly correlated
Gaussian lobe method and reported for the 5 S o state theenergy lower by 0.751 mH than the value estimated from theexperimental value [9].
We performed the FC si j calculations for the carbon 5 S o sp3
configuration (1s)2 (2s)1 (2px )1 (2py )1 (2pz )1 using the Slater
functions only and the results were summarized in Table III.
The calculations here are basically similar to those shownin Tables I and II: we optimized two sets of the parameters(α1s , α2s , d, α2p ) commonly to all the cf’s involved. After step5, since the numbers of the cf’s to be handled are larger thanthose of the above two cases, we used the cf-selection step
before increasing the order n of the calculations to keep thesize of the calculations feasible. We note that we have two
reference exact energies, one estimated from the NIST experimental values [31,32] and the other the lowest variationalenergy reported by Strasburger [30], which are referred to as
eE-a and eE-b, respectively, where “eE” denotes the exact
energy.
From step 0 to step 5, we performed the FC calculations
increasing both ri and si j operators or only ri operators,
optimizing the two sets of the parameters (α1s , α2s , d, α2p ).
We obtained the result of 2.60 kcal/mol accuracy relative to
eE-a or of 3.07 kcal/mol relative to eE-b. Then, we applied
the cf-selection method to make the size of the wave function

smaller and then increased the order n. With the cf-selection
method, we reduced the size from 3930 cf’s of step 5 tothe important 615 cf’s at the energy of 3.00 kcal/mol above
the eE-a. Then, we applied ri and si j operators and obtained
M = 4452 cf’s whose E was 1.45 kcal/mol above the eE-a.
Next, we reduced the cf size by the cf-selection method to
M = 1133 and then relaxed the radial and angular elementsof the si j operators obtaining 5059 cf’s whose E was 1.248
kcal/mol above the eE-a. Then, we selected 1490 important
cf’s and applied the ri and relaxed s1s,1s operators since theintra-1s correlations seemed to be important and obtained
9387 cf’s whose E was 0.740 kcal/mol above the eE-a and1.211 kcal/mol above eE-b. Then, we reduced the number ofthe cf’s to 2102 by the cf-selection method and applied ri andrelaxed s1s,1s operators and obtained the energy −37.690 224
au which is 0.487 kcal/mol above eE-a and 0.958 kcal/mol
above the eE-b. We repeated the same selection and the application of the ri and relaxed s1s,1s operators and obtained theenergy −37.690 418 a.u. which is 0.365 kcal/mol above eE-aand 0.837 kcal/mol above the eE-b. The energies obtained in
the eighth and ninth steps were within the chemical accuracy
from both estimates of the exact energy: this result supports
strongly the variational result obtained by Strasburger as
being the best variational energy of the 5 S o (sp3 ) state of the
carbon atom so far obtained. The difference of the two energies, −37.691 au from experimental values and −37.691 751
from the nonrelativistic variational calculation by Strasburger,
would be attributed to the relativistic and other effects [16].
We note that the costs for the FC si j calculations were muchsmaller than that of the ri j calculations, even though thecalculational process of Table III was far from smart.

062508-5

NAKATSUJI, NAKASHIMA, AND KUROKAWA

PHYSICAL REVIEW A 101, 062508 (2020)

TABLE III. Free-complement si j theory applied to the 5 S o sp3 state of carbon atom using optimized double-zeta Slater functions.a
Step

Operation

M

9ri , si jri , si jririri
cf selectionr i , si j
cf selection
si j relax
cf selectionri , relaxed s1s,1s
cf selectionri , relaxed s1s,1s
cf selectionri , relaxed s1s,1s

2
18
124
500
1540
3930
615
4452
1133
5059
1490
9387
2102
13019
2154
13814

−37.531 477−37.625 807−37.672 690−37.681 544−37.685 353−37.686 851−37.686 220−37.688 684−37.688 207−37.689 012−37.688 845−37.689 821−37.689 599−37.690 224−37.689 864−37.690 418

9

cf selection

2942
3091
3442

−37.690 162−37.690189−37.690220−37.691c (-37.691 751)d

0
1
2
3
4
5
6

7
8

Energy (a.u.)ExactE = EFC − Eexact (kcal/mol)b
100.102 (100.573)
40.909 (41.380)
11.489 (11.961)
5.934 (6.405)
3.543 (4.015)
2.603 (3.075)
3.000 (3.471)
1.453 (1.925)
1.753 (2.224)
1.248 (1.719)
1.352 (1.824)
0.740 (1.211)
0.879 (1.350)
0.487 (0.958)
0.713 (1.184)
0.365 (0.837)
0.526 (0.997)
0.509 (0.980)
0.489 (0.960)

Initial functions with double-zeta 1s, 2s, and 2p orbitals were used: φ0(k) = (1s(k) )2 αβ(2s(k) 2px (k) 2py (k) 2pz (k) )αααα (k = 1, 2). The
Slater orbitals used were 1s(k) = exp(−α1s,1 r), 2s(k) = r exp(−α2s,k r) − dk exp(−α2s,k r), 2px,y,z (k) = (x, y, z) exp(−α2p,k r) with the optimized
sets (α1s,1 , α2s,1 , α2p,1 , d1 , α1s,2 , α2s,2 , α2p,2 , d2 ) = (5.684, 1.748, 1.586, 0.03104, 6.218, 4.868, 1.844, 0.1861), (5.482, 1.670, 2.029, 0.3312,
6.561, 2.123, 1.480, 0.2145), (5.556, 1.728, 1.654, 0.2217, 10.759, 1.7634, 1.531, 0.1298), (5.610, 1.953, 1.888, 0.1248, 10.983, 1.940, 1.495,
0.2677), (5.523, 2.185, 2.001, −0.4849, 10.943, 1.986, 1.626, 0.2623), (5.905, 2.337, 2.141, −0.8030, 11.709, 2.125, 1.740, 0.1630) at step= 0 − 5. For step = 6 − 9, the same set as step = 5 was used.
bE without parentheses is from footnote c and E with parentheses is from footnote d.
cExact energy estimated from the absolute ground-state energy given in Ref. [31] and the NIST experimental excitation energies of Ref. [32].
dLowest variational energy so far obtained from Ref. [30].
aWe applied the cf-selection method to the result of theninth step and the results were shown at the bottom of
Table III. We obtained the chemical accuracy result against
eE-b as E = 0.997 kcal/mol with 2941 cf’s, which wasmuch smaller than the total number of the cf’s, 13 814 of step
9. Continuing further the cf-selection method, we obtained
E = 0.980 kcal/mol with 3091 cf’s and E = 0.960
kcal/mol with 3442 cf’s. These results of the cf-selection
method are useful since with the FC si j method, the results
can be handled only with the one- and two-electron integrals.
Thus, we could obtain a highly accurate result close to
the Schrödinger limit by using only one- and two-electron
integrals with the FC si j theory. However, in comparison withthe results obtained for the smaller atoms shown in Tables I
and II, the convergence rate for the carbon atom became slow
after step 5 reaching about 3 kcal/mol accuracy. The methodof calculations after step 5 was far from smart, though it
was partly due to the limitations of the computers we used.
Also, the wave functions obtained at the chemical accuracy
at steps 8 and 9 were not compact enough as the reference
accurate atomic wave function, though the cf-selected wave
functions obtained therefrom would have some value, since
they are relatively compact and can be handled with one- andtwo-electron integrals only. However, if we can permit therequired accuracy to be within a few kcal/mol for the absoluteenergy, the calculations of the FC si j theory would becomemuch easier: for the present case, we can use the result of
step 5. In chemistry, what we need in most studies is theaccuracy of the difference energy, not the absolute energy. For
example, the heat of formation of a molecule is the differenceof the absolute energy of a molecule relative to the sum of
the absolute energies of the constituent atoms. So, if eachabsolute energy is calculated in some constant accuracy withthe variational method, the accuracy of the difference energy
would become better. The robust nature of the variational
method is important for this type of calculations.
Now, as a last test calculation, we apply the present FC
si j theory to the hydrogen molecule in the equilibrium geometry of R = 1.4011 au using both Slater and Gaussian functions. For the H2 molecule, some two-center integrals including r-Gauss functions must be calculated by the Gaussianexpansion method [24] with the six-term expansion in thepresent case. First, we show the calculations starting from thelocal-molecular orbital (MO) initial wave function composedof the three local-MO functions given by
φ01 = (1sA + 1sB )2 αβ,
φ02 = (1sA − 1sB )2 αβ,(13)φ03 = (1sA + 1sB )(xA 1sA − xB 1sB )(αβ − βα),where 1sA or B = exp(−αrA or B ) for the Slater function and
1sA or B = exp(−αrA2 or B ) for the Gaussian function. The first

062508-6

SOLVING THE SCHRÖDINGER EQUATION OF ATOMS …

PHYSICAL REVIEW A 101, 062508 (2020)

TABLE IV. Free-complement si j theory applied to H2 molecule starting from the MO-type initial functions using optimized double-zeta
Slater and Gaussian functions.
Slater functionanM

Energy (a.u.)

0
1
2
3

6
28
94
234

3

14
19
33
39
64

−1.151 301−1.169 851−1.172 934−1.1732 932(cf selection)−1.169 928−1.171 424−1.172 889−1.173 045−1.173 223−1.174 476c

Gaussian functionb

E = EFC − Eexact (kcal/mol)nM

Energy (a.u.)

14.543
2.902
0.968
0.742

0
1
2

6
116
710

−1.045 981−1.169 399−1.173 771

2.854 (<3)
1.915 (<2)
0.996 (<1)
0.898 (<0.9)
0.786 (<0.8)

2

53
67
114
124
139

(cf selection)−1.169 852−1.171 347−1.172 894−1.173 046−1.173 203−1.174 476c

E = EFC − Eexact (kcal/mol)
80.6323.186
0.442

2.902 (<3)
1.963 (<2)
0.993 (<1)
0.897 (<0.9)
0.799 (<0.8)

a
The Slater orbitals used are 1s(k) = exp(−αk r) (k = 1, 2) with the optimized exponents of (α1 , α2 ) = (1.272, 1.654), (1.470, 2.270), (1.344,
2.330), (1.650, 2.521) at n = 0 − 3.
b
The Gaussian orbitals used are 1s(k) = exp(−αk r 2 ) (k = 1, 2) with the optimized exponents of (α1 , α2 ) = (0.3014, 1.1453), (0.1922, 0.6970),(0.2117, 1.4347) at n = 0 − 2.
cExact energy from Refs. [14,15].

one of Eq. (13) is the bonding local MO, the second oneis the two-electron excited configuration to the antibonding
local MO, and the third one is the polarization function pair.
We performed double-zeta calculations for the 1sA and 1sBorbitals: at each level of the calculations, two exponents α1
and α2 were optimized commonly to all the cf’s. Since themolecule is not spherically symmetric, we relaxed the xi x jpart from the yi y j + zi z j part of the expression given by thethird line of Eq. (10) with x being parallel to the molecularaxis. The results of the FC si j calculations starting from thelocal-MO initial functions are summarized in Table IV: theleft-hand side is for the Slater function and the right-hand side
is for the Gaussian function.
We discuss first the results for the Slater functions. Starting from the initial set of functions that was higher by
14.54 kcal/mol from the exact energy, we got at first order(n = 1) the accuracy of 2.90 kcal/mol, at second order 0.97
kcal/mol, and at n = 3, 0.74 kcal/mol accuracy, the last two
being less than 1 kcal/mol accuracy (chemical accuracy). By
performing the cf selection for the n = 3 wave function, we
found that the chemical accuracy was achieved with 33 cf’s incontrast to 94 cf’s of the second-order wave function.
We next discuss the results of the Gaussian functions. The
first impression of the results was that the energy of the initial
function was very poor: the E value was 80.6 kcal/mol
above the exact value. This is because they are the results of
two primitive Gaussian functions. However, even from such acrude initial wave function, when we apply the FC si j theory,
the accuracy was improved at n = 1 to E = 3.19 kcal/mol,which was close to the result of the Slater function, and at n =2, the result reached already the chemical accuracy of E =0.44 kcal/mol. In comparison with the Slater functions, theresults of the Gaussian functions were good if we compare at
the same order, however the number of the cf’s were muchlarger with the Gaussians than with the Slaters. When weapplied the cf-selection method to the second-order result, thechemical accuracy was obtained with the 114 cf’s; though thisnumber is much smaller than 710 of n = 2, it is still about
three times larger than the Slater-function case. Thus, as far as
the compactness is concerned, the Slater functions are muchbetter than the Gaussian functions. This is the general trendthat could be seen also from Tables I and II for He and Li
atoms, respectively.
Next, we applied the present FC si j theory to the valencebond (VB)-type initial functions of the H2 molecule given by
φ01 = 1sA 1sB (αβ − βα),φ02 = [(1sA )2 + (1sB )2 ]αβ,φ03 = (xA 1sA 1sB − 1sA xB 1sB )(αβ − βα),φ04 = xA 1sA xB 1sB (αβ − βα),(14)where 1sA or B = exp(−αrA or B ) for the Slater function and
1sA or B = exp(−αrA2 or B ) for the Gaussian function. φ01 represents the covalent function, φ02 the ionic function, φ03 thepolarization function, and φ04 the double polarization function.
Table V shows the results. In this calculation, we handled
φ01 as being the zeroth-order function and the other threefunctions as first-order functions: the cf’s at order n of Table V
is the sum of the order n cf’s produced from φ01 and the order(n − 1) cf’s produced from φ02 , φ03 , and φ04 .
When we use the Slater-type functions, we could get
accuracy of less than 2 kcal/mol at order 2 with 34 cf’s andthe chemical accuracy at n = 3 as E = 0.74 kcal/mol. By
applying the cf-selection method to the order 3 result, wefound that actually only 30 cf’s were necessary for getting thechemical accuracy. We also showed the results for the higher
accuracies. The right-hand side of Table V gives the results
for the Gaussian functions. Though the energy of the covalent
double-zeta Gaussian functions were very crude (E = 89.7
kcal/mol), the FC theory improved the energy to E = 1.62
kcal/mol at n = 1 and we could get the chemical accuracy
already at n = 2. Probably, this is due to the fact that we had amuch larger number of functions with the Gaussians than with
the Slaters as we have seen similar trends in Tables I and II.

062508-7

NAKATSUJI, NAKASHIMA, AND KUROKAWA

PHYSICAL REVIEW A 101, 062508 (2020)

TABLE V. Free-complement si j theory applied to H2 molecule starting from the VB-type initial functions using optimized double-zeta
Slater and Gaussian functions.
Slater functionanM

Energy (a.u.)

0
1
2
3

2
12
34
72

3

14
15
30
39
52

−1.141 658−1.161 817−1.171 340−1.173 289(cf selection)−1.170 109−1.171 647−1.172 890−1.173 043−1.173 203−1.174 476c

Gaussian functionb

E = EFC − Eexact (kcal/mol)nM

Energy (a.u.)

20.593
7.944
1.968
0.745

0
1
2
3

2
32
286
1597

2.740 (<3)
1.775 (<2)
0.995 (<1)
0.899 (<0.9)
0.799 (<0.8)

3

43
61
97
103
114

−1.031 520−1.171 887−1.173 088−1.173 939(cf selection)−1.169 774−1.171 325−1.172 884−1.173 047−1.173 205−1.174 476c

E = EFC − Eexact (kcal/mol)
89.706
1.624
0.871
0.337
2.951 (<3)
1.977 (<2)
0.999 (<1)
0.897 (<0.9)
0.797 (<0.8)

a
The Slater orbitals used are 1s(k) = exp(−αk r) (k = 1, 2) with the optimized exponents of (α1 , α2 ) = (1.3056, 1.5452), (1.4465, 1.7124),(1.3282, 2.0011), (1.6589, 2.7933) at n = 0 − 3.
b
The Gaussian orbitals used are 1s(k) = exp(−αk r 2 ) (k = 1, 2) with the optimized exponents of (α1 , α2 ) = (0.2981, 1.0334), (0.2201, 0.8791),(0.2337, 1.3185) at n = 0 − 2. At n = 3, (α1 , α2 ) = (0.2337, 1.3185) which is the same as that of n = 2 is used.
c
Exact energy, Refs. [14,15].

At n = 3, we could get a highly accurate result of E =0.34 kcal/mol. When we apply the cf-selection method tothis result, we see that only 97 functions were necessary for
obtaining the chemical accuracy.
In comparison with the MO-type calculations shown in
Table IV, the present VB-type results seem to be better thanthe MO-type ones for both Slater and Gaussian functions,even though the order-counting of the initial functions wasfavorable for the local-MO calculations.
Finally, we compare the results of the H2 molecule withthose of the He atom, both two-electron systems. From theorder that is necessary to get the chemical accuracy, both
calculations were roughly similar, but when we compare fromthe necessary numbers of the cf’s for getting the chemical
accuracy, the H2 calculations were less efficient than the He
calculations: for getting the similar accuracy, about a twice
larger number of cf’s was necessary for the H2 calculations
than for the He calculations, for both the MO-type and VBtype calculations.
IV. CONCLUSION

Based on the scaled Schrödinger equation [3], we haveproposed the FC si j theory with si j being ri2j as a variational
theory that requires only one- and two-electron integrals forsolving the SE in some accuracy, in contrast to the original
theory referred to as the FC ri j theory that requires us toevaluate three- and more-electron integrals. The si j operator
can prevent the divergence difficulty of the integrals, likewise
the ri j operator, though si j is less efficient than ri j because
si j behaves less satisfactorily than ri j in particular near theelectron-electron cusp region. An important merit associatedto si j is that it can be rewritten only with the one-electron
functions as shown by Eq. (10), so that the wave function
produced by the FC si j theory is composed of only oneelectron functions. Therefore, the variational calculations of
the FC si j theory can be performed within one- and two-electron integrals: the two-electron integrals come from theintegrals of the 1/ri j operator included in the Hamiltonian.
Because this merit is attractive, we have examined in thispaper how accurately the FC si j theory can solve the SE. By
applying to the He, Li, 5 S o (sp3 ) state of the C atom and H2
molecule using Slater and Gaussian functions, we found thatas the order of the theory was increased, the variational FC si jtheory improved the solutions from a few kcal/mol accuracy
to the chemical accuracy of less than 1 kcal/mol error from thereference exact energy by using only one- and two-electron
integrals. The convergence rates to the exact limits were slow
in comparison with the FC ri j calculations performed before,
but the computational costs were much smaller. Judging fromthe order n of the FC theory, the Gaussian functions converged
faster than the Slater functions, but from the number of the
cf’s necessary, the Slater functions were more compact thanthe Gaussian functions. For the 5 S o (sp3 ) state of the C atom,
we could get the result that supports the recent calculations
performed by Strasburger [30]: the nonrelativistic variationalenergy of the 5 S o (sp3 ) state of the C atom is lower than
−37.691 a.u. but above −37.691 751 a.u. Though the present
calculation was performed as a test calculation of the FC si jtheory, the result was already at a state-of-the-art level. The
local-MO- and VB-type calculations of the H2 molecule with
the Slater and Gaussian functions showed similar behaviors to
those performed for atoms. However, the convergence rate to
the chemical accuracy was about twice slower than that for
the He atom, the same two-electron system, when counted
from the numbers of the necessary cf’s. Thus, with the FC
si j theory, we could obtain the chemical accuracy of less than
1 kcal/mol for all the systems examined here.
However, we have also observed for the 5 S o (sp3 ) state of
the C atom that the convergence rate became slow when thecalculation came close to ∼1 kcal/mol accuracy: up to a few
kcal/mol accuracy, the calculations were straightforward, butthereafter, the calculations were not smart. However, if we
can permit the required accuracy of the FC si j calculations

062508-8

SOLVING THE SCHRÖDINGER EQUATION OF ATOMS …

PHYSICAL REVIEW A 101, 062508 (2020)

to be within a few kcal/mol for the absolute energy, the calculations would become much easier. Actually, most chemical
phenomena are governed by the difference energy, not theabsolute energy, involved in the phenomena. Therefore, it isnecessary to examine the accuracy of the difference energy inthe variational calculations of the FC si j theory. If the above
relaxation in the necessary accuracy of the theory is permitted,
the practical calculations would become easier. The robust
nature of the variational method would be helpful in doing
such computational research.
Thus, we would like to suggest two lines of future studies.
First, it is still very important to investigate how to improve
the theory itself to obtain easily the chemical accuracy solutions of the SE for the absolute energy. Second, we must
investigate how much is the accuracy of the difference energy
obtained with the variational calculations of the FC si j theory
applied to real chemical phenomena.

APPENDIX: SELECTION OF IMPORTANT
COMPLEMENT FUNCTIONS

We thank N. Kawakami for his continuous support of theactivities of QCRI. This research was also supported by the
JSPS KAKENHI, Grants No. 16H02257, No. 17H06233, No.
17H04867, and No. 17KT0103. The authors thank the Research Center for Computational Science, Institute for Molecular Science, Okazaki, for letting us use their computers.

After a set of the FC si j calculations, it is useful to havean efficient method of selecting the important complement
functions (cf’s) in the orders of their energy contributions.
This is necessary for obtaining a compact set of cf’s that gives
some energy accuracy and for continuing the series of the FC
si j calculations without making the dimensions too large forhandling.
The cf-selection method is designed to produce a compact
set from a diffuse one by choosing the important cf’s, one
by one, from the diffuse set without much loss of accuracy.
As we did in the present calculations, we obtain the result of
the desired accuracy, chemical accuracy for example, by increasing the order n of the FC theory with a diffuse dimension
M. Then, we try to produce a compact set by the systematic
selections of the cf’s involved in the diffuse set. The algorithm
is as follows. (i) First, we choose the initial function(s) φ0
as the most important cf (set). (ii) Add ith cf (i = 0) chosenfrom the diffuse set to φ0 , solve the eigenvalue problem of thedimension 2, and obtain the energy Ei (i < M ). Repeat it for
all the M − 1 candidates of the diffuse set. (iii) Select the cf
φselect that gives the largest energy lowering using the energy
set {Em } and replace φ0 by φ0 + φselect . (iv) Repeat (ii) and (iii)
until we get the desired accuracy for the system. (v) When we
get the desired accuracy at the kth iteration, the first k + 1 cf’s
are the compact set we want to obtain.

[1] P. A. M. Dirac, Proc. R. Soc. Ser. A (London) 123, 714(1929).
[2] H. Nakatsuji, J. Chem. Phys. 113, 2949 (2000).
[3] H. Nakatsuji, Phys Rev. Lett. 93, 030403 (2004).
[4] H. Nakatsuji, Phys. Rev. A 72, 062110 (2005).
[5] H. Nakatsuji, H. Nakashima, Y. I. Kurokawa, and A. Ishikawa,Phys. Rev. Lett. 99, 240402 (2007).
[6] H. Nakatsuji, Acc. Chem. Res. 45, 1480 (2012).
[7] H. Nakatsuji and H. Nakashima, J. Chem. Phys. 142, 194101(2015).
[8] H. Nakatsuji, H. Nakashima, and Y. Kurokawa, J. Chem. Phys.
149, 114105 (2018).
[9] H. Nakatsuji, H. Nakashima, and Y. Kurokawa, J. Chem. Phys.
149, 114106 (2018).
[10] Y. Hijikata, H. Nakashima, and H. Nakatsuji, J. Chem. Phys.
130, 024102 (2009).
[11] H. Nakatsuji, J. Chem. Phys. 83, 5743 (1985).
[12] H. Nakashima and H. Nakatsuji, J. Chem. Phys. 127, 224104(2007); Y. Kurokawa, H. Nakashima, and H. Nakatsuji, Phys.
Chem. Chem. Phys. 10, 4486 (2008).
[13] H. Nakashima and H. Nakatsuji, Phys. Rev. Lett. 101, 240406(2008).
[14] Y. Kurokawa, H. Nakashima, and H. Nakatsuji, Phys. Rev. A
72, 062502 (2005).
[15] Y. Kurokawa, H. Nakashima, and H. Nakatsuji, Phys. Chem.
Chem. Phys. 21, 6327 (2019).
[16] H. Nakatsuji and H. Nakashima, J. Chem. Phys. 150, 044105(2019).

[17] V. Saunders, in Handbook of Molecular Physics and Quantum
Chemistry (J. Wiley & Sons, New York, 2003), Vol. 2, Chap.
29, p. 641 and Chap. 30, p. 682.
[18] Recent Advances in Computational Chemistry: MolecularIntegrals over Slater Orbitals, edited by T. Ozdogan and
M. B. Ruiz (Transworld Research Network, Kerala, India,
2008).
[19] S. Huzinaga, Suppl. Prog. Theor. Phys. 40, 52 (1967).
[20] H. Nakashima and H. Nakatsuji, QCRI atomic integral program
package (unreleased).
[21] H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schutz,
P. Celani, W. Gyorffy, D. Kats, T. Korona, R. Lindh et al.,MOLPRO, version 2019.2; J. F. Rico, R. Lopez, G. Ramirez,
I. Ema, D. Zorrilla, and K. Ishida, SMILES a package formolecular calculations with Slater functions.
[22] QCRI Gaussian molecular program (unreleased).
[23] H. Nakatsuji and H. Nakashima, Intern. J. Quantum Chem. 109,
2248 (2009).
[24] Y. Kurokawa, H. Nakashima, and H. Nakatsuji (unpublished).
[25] K. O-ohata, H. Taketa and S. Huzinaga, J. Phys. Soc. Jpn. 21,2306 (1966); H. Taketa, S. Huzinaga, and K. O-ohata, ibid. 21,2313 (1966).
[26] M. Lesiuk and R. Moszynski, Phys. Rev. E 90, 063318 (2014);
90, 063319 (2014).
[27] D. W. Schwenke and H. F. King, J. Phys. Chem. A 120, 9348
(2016).
[28] J. Noga, S. Kedzuch, J. Simunek, and S. Ten-no, J. Chem. Phys.
128, 174103 (2008).

ACKNOWLEDGMENTS

062508-9

NAKATSUJI, NAKASHIMA, AND KUROKAWA

PHYSICAL REVIEW A 101, 062508 (2020)[29] M. D. Brown, J. R. Trail, P. Lopez Rios, and R. J. Needs,
J. Chem. Phys. 126, 224110 (2007).
[30] K. Strasburger, Phys. Rev. A 99, 052512 (2019).
[31] S. J. Chakravorty, S. R. Gwaltney, E. R. Davidson, F. A. Parpia,
and C. F. Fischer, Phys. Rev. A 47, 3649 (1993).

[32] C. E. Moore, Atomic Energy Levels: As Derived fromthe Analyses of Optical Spectra (National Bureau ofStandards, Washington, DC, 1952), GPO, Vol. I, p. 227;
see https://physics.nist.gov/ for NIST Atomic Spectra
Database.

062508-10