Variational Methods for Atoms and the Virial Theorem

: In the case of the one-electron Dirac equation with a point nucleus, the virial theorem (VT) states that the ratio of the kinetic energy to potential energy is exactly − 1, a ratio that can be an independent test of the accuracy of a computed solution. This paper studies the virial theorem for subshells of equivalent electrons and their interactions in many-electron atoms. This shows that the linear scaling of the dilation is achieved through the balancing of the contributions to the potential of an electron from inner and outer regions that some Slater integrals impose conditions on a single subshell, but others impose conditions between subshells. The latter slows the rate of convergence of the self-consistent ﬁeld process in which radial functions are updated one at a time. Several cases are considered. Results are also extended to the nonrelativistic case.


Introduction
In classical mechanics, the virial theorem (VT) is well-known [1]. It relates the total kinetic energy of a stable system of discrete particles bound by potential forces with that of the total potential energy of the system. The validity of this theorem for atoms was proven by Finkelstein [2] for the Schrödinger equation in 1928. Fock [3,4], who used the variational method to derive the Hartree-Fock method for atoms, was the first to derive the theorem for the Dirac equations using a stretching method also referred to as a "scaling" [5] or "dilation" procedure [6].
The nonrelativistic wave equation for an N-electron atomic system, HΨ = EΨ, is often expressed in Cartesian coordinates along with approximate solutions, Ψ HF , obtained with the variational Hartree-Fock method [7]. In this approximation, the multielectron wave function is expressed as a determinant of one-electron wave functions or more generally as a linear combination of antisymmetrized products of orbitals, one for each electron, of the same form as wave functions of the one-electron case. As illustrated below, replacing 1/r ij with an expansion in terms of Legendre polynomials, namely, where r < = min(r i , r j ) and r > = max(r i , r j ), introduces Slater integrals into the energy expression [8,9]. The virial theorem is a special test of the wave function. When the spatial variables of a normalized wave function are inflated or scaled by the transformation r → λr for all electrons, then for a stationary solution, the energy is still stationary, namely, ∂( Ψ(λ)|H|Ψ(λ) )/∂λ = 0.
This stationary condition holds for the exact solution of the wave equations, but it also holds for a solution of the variational equation itself. This stationary condition that leads to the virial theorem (VT) ratio. Thus, it is a property of a stationary solution along with other theorems such as Brillouin's [10,11].
In atomic calculations, Ψ HF is expressed in spherical coordinates as a vector-coupled configuration state function (CSF) for groups of equivalent electrons with given total spinangular quantum numbers and parity in nonrelativistic theory or total J and parity in relativistic theory, where equivalent electrons have the same radial functions. In spherical coordinates, the wave function for a group of equivalent electrons is a product of radial and spin-angular factors. In nonrelativistic theory, the radial factor is simply Π i P a (r i )/r i , i = 1, w a , where a identifies the group of equivalent electrons through the nl quantum numbers, and w a is the occupation number of the subshell. In Dirac theory, the radial factor is a 2 × 2 matrix with diagonal elements that are similar products of large (P a (r i )) and small (Q a (r i )) components, along with a two-component vector of spin-angular components , where a refers to the nκ quantum numbers, with κ = −(l + 1) and +l for j = l + 1/2 and j = l − 1/2, respectively [12]. Multiplying the wave equation by Ψ † and integrating over all coordinates, we obtain an energy functional (expression) in terms of the unknown radial functions. In both approximations, the energy functional for Ψ (D)HF is a list of one-electron integrals (I(a, a)) and two-electron Slater integrals, so that where a refers to a subshell of w a equivalent electrons, and b also refers to a subshell. The direct F k (a, b) and exchange G k (a, b) integral are particular cases of the generalized Slater integral R k (ab, cd) defined in Section 6.1 Coefficients f ab,k and g ab,k of the direct and exchange radial integrals, respectively, result from the integration of Coulomb Operator (1) over spin-angular coordinates. This integration is highly selective [13] and limits the value of k that defines the Slater integrals for a given state.
In this section, we only consider the Dirac-Hartree-Fock (DHF) case. In multiconfiguration expansions where Ψ is a linear combination of CSFs, Slater integrals of different symmetries also appear. Scaling affects only the radial factor of the CSF wave function and not the spin-angular factors.
For multielectron atoms, the variational method for relativistic orbitals leads to a system of equations, one for each orbital a with quantum numbers nκ, whose radial functions are varied. In the numerical multiconfiguration Dirac-Hartree-Fock method (MCDHF), these equations have the following form: where c is the constant used for the speed of light, u a is a two-component vector of large (P a (r)) and small (Q a (r)) components of the radial functions, i.e., u a (r) = (P a (r), Q a (r)) t , V(r) is the potential, a is the orbital energy, and S is a 2 × 2 identity matrix in the r variable. In numerical methods, for cases with two or more electrons, potential V(r) includes only the direct interactions, whereas the exchange contributions are included in the two-component functions X(r) along with contributions of the type ab u b (r)δ κ a ,κ b from off-diagonal Lagrange multipliers that ensure the orthogonality of orbitals of the same κ symmetry. The numerical solution of these equations involves the matching of outward and inward integration procedures. For more details, see [6,12,14,15]. In B-spline methods in which the radial functions for large and small components of one-electron spinors are expanded in B-spline bases, all variational contributions to the energy expression are included in a matrix, orbitals rotated for stationary energy, and Lagrangian multipliers eliminated through the use projection operators [16,17]. Orbital energy a is then an eigenvalue of an interaction matrix.
It is not obvious from the form of Equation (5) that, as c → ∞, the solution approaches the nonrelativistic limit with Q a (r) → 0. This is clearer if factor c is included in the definition of the orbital, i.e., if Q a (r) was replaced by cQ a (r), so that the second column of the matrix was divided by c and the second equation was also divided by c. In the nonrelativistic limit, the equation becomes the radial Schödinger equation but with the difference that the Dirac equation includes a "mass energy" correction. Thus, the virial theorem (VT) ratio is V / T = −1, in the relativistic case, and −2 in the nonrelativistic case for exact solutions of the variational equations, where V is the potential energy, and T the kinetic energy.
The virial theorem (VT) is a special test of the computed wave function involving the scaling of the radial functions. When numerical methods are used, the tests simply regard how well the differential equations are solved, but they are also a test of how accurately the algorithms have determined Lagrangian multipliers, and with which assumptions or simplifications. There are cases where Lagrangian multipliers can be set to zero, the wave function remains unchanged, and energy is stationary. At the same time, there are cases where Lagrangian multipliers are sufficiently large so that radial functions have extra nodes making node counting an art. A good case for study are VT results from the B-spline Dirac-Hartree-Fock program DBSR_HF [16] and the nonrelativistic B-spline Hartree-Fock SPHF code [18] for the 1s2s 1 S case shown in Table 1. Orthogonal transformations can be expressed as orbital rotations. In the DBSR_HF code, the default is no rotations and off-diagonal Lagrangian multipliers set to zero; however, as seen in Table 1, this result has a rather poor ratio, indicating 2-3 digits of accuracy, whereas with rotation, the same numerical methods are almost as accurate as those for the nonrelativistic B-spline method. 5d 10 6s 2 6p 6 using the DBSR_HF program to high precision. To improve accuracy, the first nonzero grid point was set to be 0.00001/Z, the exponential grid-step parameter was reduced to he=0.125, and the convergence for the energy was set to scf_tol=1.d-16, the change in the largest value of the radial function to orb_tol=1.d-10, and the tolerance for the tail region to end_tol=1.d-10. The first two tolerances were for a relative change. The entire calculation was performed in double precision arithmetic with 15-16 digits of numerical accuracy. In other words, the requested accuracy was machine precision for the total energy. We see immediately that the VT sum grew in magnitude as the number of shells and total energy increased, except for the value for Kr with the 3d 10 subshell where the value was especially large. In atomic spectroscopy, the digits after the decimal point are important, as wavelengths depend on energy differences that are usually a fraction of the E h unit. The DBSR_HF also reports the VT ratio that is related to the accuracy of the total energy, which was small for Rn because the total energy was large.
Both Kim [19], and Lindgren and Rosén [20], used the r k < /r k+1 > operator and associated Slater integrals in their analysis of the relativistic Hartree-Fock equations, but relied on the 1/r ij operator for deriving the VT theorem. In this paper, we explore scaling to see how the VT results are achieved when Slater integrals are used. We first consider the scaling of the radial equations of associated subshells of equivalent electrons, including the scaling of Slater integrals, and then relate the sum of these equations to the scaling of the total energy. We refer to the stationary condition for the orbital equation as VT a , whereas VT refers to that of the total energy.

General Theory
In order to analyze the effect on the energy of a single electron in a potential from the scaling of the radial function, it is convenient to write Equation (5) in operator form, namely, where where T(r) is the kinetic energy operator, V(r) the potential energy operator, M is a constant matrix referred to as the mass operator, and S is the identity matrix.
Consider scaling perturbation r → λr of the radial function, so that Equation (6) becomes where a is now a function of λ. Before deriving an expression for the energy parameter, let us descale the entire radial equation by replacing r by r/λ, so that the equation becomes The scaled radial function went back to the original function, whereas the operators in the radial equation changed.
In Equation (9), function u a (r) is normalized; so, by multiplying on the left by u t a (r) and integrating, it follows that For a stationary solution u a (r) with respect to scaling parameter λ, the energy satisfies variational condition In order to continue the analysis, let us consider a continuous operator O(r) and define the expectation value of the operator to be (12) and O a ≡ O a (λ = 1). Then, according to Equations (8), (10), and the variational condition of Equation (11): Thus, we show that the scaling of the energy depends only on the scaling of T(r) and V(r).
The scaling of T(r) is simple. Clearly κ/r → λ κ/r as r → r/λ. Operator d/dr scales in the same manner, so that Hence, However, the scaling of V(r) is more complicated. In general, since the differentiation with respect to λ and the integration with respect to r commute, we know that Now let r = r/λ. Then, Hence, with expectation values, where we express the final scaling as V a and a correction, W a , the deviation from linear scaling. The above −rd/dr operator was first introduced by Fock [3,4], and we refer to it as the Fock rule. This rule is useful when scaling is not straightforward. The discussion has so far assumed a single electron in a potential. We extend this to a single electron in a Hartree potential for equivalent electrons in a given subshell, and lastly an orbital in a Hartree-Fock potential for multiple subshells of equivalent electrons, all of which are scaled. We show that, in general, the scaling condition of Equation (14) resulting from the stationary condition of orbital energy Equation (11) with respect to the variation of λ (at λ = 1)) has the form where R a may be zero for exact solutions in the case of a single orbital, but may not be zero in the case where the variation includes multiple subshells. Expressions for W a appearing in Equation (20) are determined as a correction arising from the nonhydrogenic potential.
We refer to this equation as the virial theorem of an orbital (VT a ) to distinguish it from the VT for the total energy for systems with multiple orbitals and show that ∑ a w a R a = 0 for a solution that satisfies the virial theorem. Here, w a is the occupation number of shell a.

Nonrelativistic Case
In nonrelativistic theory, small component Q a (r) = 0. Consequently, M = 0 and can be omitted. Though the values of Slater integrals are different, the scaling is the same. The nonrelativistic equation for an electron in a potential is Thus, the kinetic energy operator is As a result, the energy equation and the stationary condition for the dilation are the equivalent of Equations (13) and (20), respectively. For the total energy, the VT is usually written as

Point Nucleus: V (r) = −Z/r
In the special case of a single electron in a potential of an atom or ion with a point nucleus, V(r) = −Z/r, where Z is a constant representing the nuclear charge. Hence, like the kinetic energy in the relativistic case (see Equation (15)), V a (λ) = λ V a . Taking the partial derivative with respect to λ and then setting λ = 1 and R(r) = 0 for an exact solution for a single electron, Fock [3,4] obtained the following VT a , Essentially, the kinetic energy and the potential energy are equal in magnitude but opposite in sign, with the kinetic energy being positive, and with the sum being a test of accuracy of the computed solution. In this case, the potential function is homogeneous of degree −1 (i.e., scales as λ) [19] and W a = 0.
Consider some examples. The recently proposed integration method for the Dirac equation [15] based on Simpson's rules solves the hydrogenic equation for a point nucleus to almost machine precision for high Z, as shown in Table 3. There are two more significant digits for hydrogen-like fermium, Z = 100, than there are for the neutral hydrogen atom, which suggests that the sum represents an absolute rather than relative error. Table 3. Comparison of kinetic T 1s and potential V 1s energies for a 1s electron for Z = 1 and Z = 100 hydrogen-like atoms. For an exact solution, T 1s + V 1s = 0 but computationally T 1s + V 1s = R 1s .

General Case: V (r) = −Z(r)/r
In the case of multiple equivalent electrons, the effective nuclear charge varies as a function of r as in the Hartree-potential for a subshell or in the case of a finite nucleus where the potential near the nucleus is modified. For such a situation, where the potential can be represented as V(r) = −Z(r)/r, the Fock rule yields Thus, the scaling of this potential introduces a new type of term: This depends on a derivative of Z(r). Therefore, it is not affected by adding or subtracting a constant. If Z(r) = Z − s(r), then W a = − s (r) a , where s(r)/r represents the deviation of V(r) from the Coulomb potential −Z/r.
Given these quantities, the virial theorem of an electron in a potential remains that of Equation (20), but now W a is not zero. Adding and subtracting W a to orbital energy Equation (13), we obtain With Equation (20), the sum of the first three terms is R a , and we obtain the second orbital virial relationship, namely, Let us now consider two typical examples of potentials for a single subshell.

Finite Volume Nucleus
In the case of a finite nucleus, Matsuoka and Koga [21] took a slightly different point of view. Their interest was primarily in the effect of a finite nucleus on the virial theorem for the total energy of many-electron systems as a correction to energy. In relativistic calculations for heavy elements, variational calculations need to treat the finite volume as a correction that modifies the potential for all electrons [22,23].
Traditionally, for a finite nucleus, the potential is defined as follows: where R N is the nuclear radius. Hence, function Z (r) would be nonzero only in the region r < R N . Their formula for the correction to the virial theorem was equivalent to W a = − s (r) using the uniform model, integrating over (0, R N ).
In the case of a Fermi nucleus [24], the potential and low-order derivatives are continuous at r = R N . In particular, −s (r) → 0 as r → R N and differentiation and/or integration can be over the range (0, ∞). This form is useful when both one-electron and many-body effects are included in the potential.

The 1s 2 1 S 0 Case
As in the hydrogenic case, the variational equation for this two-equivalent-electron system is homogeneous, in that there is no exchange contribution in Equation (5). Unlike the finite-nucleus effect that modifies only part of the orbital range, the potential for the 1s orbital is modified over the entire range by a contribution from the Slater integral F 0 (1s, 1s) in the energy expression, but the virial equations have exactly the same form. The related contribution to W 1s can be computed directly by differentiating the potential for the 1s orbital. There is only one orbital, and in this case, R 1s = 0 to numerical accuracy. Table 4 shows the various contributions to the energy and the VT equations for these two cases -the 1s potential of hydrogen-like fermium with a Fermi nucleus, and the heliumlike fermium 1s 2 case with a point nucleus. The very common Fermi nucleus [22,24] used in GRASP [25] was used in the first case since unlike uniform distribution, the potential has continuous derivatives. The results here were obtained using the recently developed numerical procedures [15]. For completeness, we also report the total energy E(1s 2 ) = −2 1s − F 0 (1s, 1s).
These two cases are important test cases for any code of relativistic, computational, and atomic structures. There is only one orbital and the energy expression includes no exchange contributions. For a system with one orbital, R a = 0 for an exact solution. Table 4. Energy contributions (in E h ) to virial theorems for hydrogen-like (1-electron, Fermi nucleus) and heliumlike (2-electron, point nucleus) fermium (Z = 100, A = 257) ground states. Parameters for the Fermi nuclear charge distribution, ρ(r) = ρ 0 /(1 + e (r−b)/a ), were r rms = 5.8756 fm, a = 0.523387555 fm, b = 7.170561722 fm. All quantities are for a single 1s orbital.

DHF Equation for Equivalent Electrons
The presented theory has so far been discussed strictly as an electron (a) in a potential, but it can readily be extended to a subshell (e.g., α) of w a -equivalent electrons by multiplying the orbital properties by the occupation number. For example, for the 1s 2 subshell containing two 1s electrons, we have Multiplying Equation (20) by w 1s = 2, we obtain the VT α for the 1s 2 subshell, namely, In this form, the radial equations for the orbitals are derived from the variation of a specific orbital without dividing by the occupation number for the subshell, as is often performed for the derivation of the Hartree-Fock equations for multisubshell atoms.

Multielectron Systems with Two or More Subshells
In order to derive expressions for the various kinetic, potential, and other properties of subshells, let us now individually consider the scaling of each of the contributions from radial integrals and their coefficients that appear in the energy functional of a state, and their contribution to the radial equations, and then derive the scaling equation for the total energy by summing over all subshells. This approach includes both direct and exchange contributions to the potential. For an exact (self-consistent) solution for which ∑ a w a R a = 0, we now have ∑ α R α = 0 by Equation (32).
In our analysis, we consider Dirac-Hartree-Fock energy expressions that include only one-electron integrals I(a, a) = T a + V 1 a + M a and Slater integrals F k (a, b) or G k (a, b) that are 2-electron contributions to V 2 α where Given an orbital basis and an energy expression, the kinetic and potential energies for subshells can readily be computed directly without considering orbital potentials.
The one-electron integrals I(a, a) multiplied by their coefficients contribute directly to T α , V α , and M α , whereas Slater integrals contribute to two subshells, when β = α or twice to subshell α when β = α.

Slater Integrals
The relativistic Slater integral, usually denoted as R k (ab, cd), is a two-dimensional integral with two coordinates, e.g., (r, s). Let us denote radial factors in terms of vector products of large and small components ρ(ac; r) = u t a (r)u c (r), and ρ(bd; s) = u t b (s)u d (s).
In this definition, ρ is defined by the orbitals of the first coordinate and the contribution to the potential by the orbital of the second coordinate. Slater integrals have many symmetries arising from the symmetry of a product and the symmetry of the coordinate system. In the canonical form, orbitals are in a designated order, such that a <= c, b <= d, a <= b and, if a = b, then also c <= d. Consequently, a typical canonical integral is by the symmetry of the coordinate system. Because of this symmetry, we consider only the first case in our study. Each Y k function itself is the sum of two integrals. In particular, where one is integrated over s < r and one over s > r as indicated below: Because variable r only defines the range of integration of the X integrals, the sum of derivatives of X k < and X k > obeys a useful relation r −k X < (bd; r) + r k+1 X > (bd; r) = 0 . (36)

Scaling of the V k (bd; r) Operator
For the variation of a contribution to the potential from orbitals of the second coordinate, we need to consider the contribution to a potential , for example, Factors r −(k+1) and r k scale respectively as λ (k+1) and λ −k , and given Equation (36), direct differentiation is more straightforward than the Fock rule (18) and yields ∂V k (bd; r/λ) ∂λ λ=1 = V k (bd; r) + W k (bd; r) by definition Then, it follows that W k (bd; r) = kV k < (bd; r) − (k + 1)V k > (bd; r) (38)

Symmetry of Slater Integrals in DHF Energy Expressions
The above derivation was a general derivation for any combination of orbitals. In DHF calculations, only R k (ab, ab) or R k (ab, ba) symmetries occur. Then, if d = b and c = a, the contribution to scaling W α would be Here, we introduce subscript aa to denote the ρ(r) function required for the calculation of the expectation value. Scaling in the other coordinate yields a contribution to W β In the case of a G k (ab) Slater integral, there is only one possibility: Y k (ab; r) and ρ(ab; r), and integral would contribute to both W α and W β , multiplied by the coefficient of the Slater integral.
In this notation, a Slater integral can be expressed as an expectation value of a potential function, namely,

Total Energy Virial Equation
Computationally, the list of one-electron integrals and Slater integrals along with their coefficients defines the total energy of a state and also determines the subshell quantities T α , V α , M α , and W α . Then, are total values that contribute directly to the total energy and its dilation equation. However, the sum of potential energies of subshells is different. Because the Slater integrals with orbitals a, b from subshells α, β contribute to both V α and V β (or twice to V α when b = a), we have where V is the potential energy for the total energy and in both coordinates, and V 2 represents the contribution from the 2-body Slater integrals. Then, we have two equations: Consequently, T + V = R − U . Substituting into Energy Equation (45), we obtain mass dilation condition Thus, an exact solution requires that R = U = 0, in which case there are two virial conditions, namely, T + V = 0 and M = E.

Slater Integral Scaling Condition
Condition U = 0 is satisfied if for each Slater integral the dilation contributions from both coordinates, and the Slater integral is zero. There are three cases that need to be considered:

1.
F k (a, a): In Equation (38) the scaling is the same in both cooordinates and expressing F 0 (a, a) in terms of potentials; then, it follows that 2 W k (aa; r) aa + F k (a, a) = (2k by Equation (38).

2.
G k (a, b): In this case, the contribution is also the same for each coordinate, and the sum (involving both subshells) is 3. F 0 (a, b): In this case, the scaling depends on the coordinate, and the effect of scaling for orbital a differs from that for orbital b. However, the combined contribution to the dilation equation for the total energy should still be zero, namely, In this case, the scaling within a given coordinate does not need to be zero, although the sum for the two coordinates should be zero.
In summary, the scaling or dilation conditions for (Dirac)-Hartree-Fock calculations for the three types of Slater integrals are: Thus, the stationary conditions for dilation require a balance of the contributions to the potential from the inner (s < r) and outer (s > r) regions. Particularly intricate is the condition for the F k (a, b), b = a condition, involving two different charge densities.

Results and Analysis
Traditionally, the virial theorem for the total energy is computed (in sum form) simply as T + V , which would be R − U in the formalism of Equation (46). In order to gain insight into how the Slater integral method achieved linear scaling, we predicted the scaling contribution to the stationary condition of each integral and computed T + V + U , which showed the extensive balancing of contributions of an orbital potential from the charge distribution. However, the most extensive analysis of the orbital equations for an atomic system is the study of the scaling of orbital equations leading to the calculation of the total energy.  [25])). We report residuals for subshell equations, R α = T + V + W α , the dilation correction for the total energy U , and R , the residual of the dilation equation for the total energy. See text for details.  Table 5 shows results for the ground states of ions of Fm (Z = 100) as more closed shells are added with values for He-like, Be-like, and Ne-like studies using radial functions from DBSR_HF [16] for standard (not high-precision) calculations, and Ar-like with radial functions from GRASP [25]. In each case, radial functions were transformed to the numerical grid of the revised procedures [15], and the various components were computed and analyzed. Displayed is information about the scaling of each individual orbital subshells, the most important being R α . For He-like, the value is small but somewhat larger than the value reported in Table 4 because the present values are not high-precision results.
Here, ρ(1s, 1s; r) is both a weighting factor for the expectation value and the charge density function that determines a contribution to the potential as a function of r.
Another interesting case is Be-like 1s 2 2s 2 1 S. The F 0 (1s, 2s) Slater integral connects two different orbital equations and the dilation conditions (Equation (51) are distributed to both equal size and opposite sign. When summed for the total energy, these contributions are cancelled out. In larger systems with multiple subshells, similar cancellations are involved. Table 6. Comparison of virial theorem analysis for a fully variational calculation for Be 1s 2 2s 2 and a fixed core calculation where the 1s orbital is from the ion Be +2 1s 2 . Frequently, when spectral calculations are performed for multiple levels, a fixed core approximation is used where the orbitals for the core are not varied. Table 6 shows highprecision virial results for the case where the 1s orbital is fixed at the value for the Be +2 1s 2 , and only the 2s orbital is varied. In this fixed core approximation, the R α values no longer cancel. The fixed core approximation is often desirable in the case of heavy elements. Codes such as DBSR-HF should report the virial theorem results only for the orbitals that are varied, which means that the interaction of varied subshells with the fixed core subshells needs to be omitted in analysis.

(i) Fully Variational
Some spline methods for nonrelativistic variational equations, updating a few or all orbitals simultaneously, were investigated by Froese Fischer et al. [26]. This is a quadratically convergent Newton-Raphson method that has achieved better performance. In their study, the calculation for Be 1s 2 2s 2 converged in 4 iterations with a virial theorem deviating from −2 by 10 −13 , considerably fewer than the typical 12 SCF iterations. Newton-Raphson methods have not yet been applied to atomic relativistic equations 1 but are expected to perform similarly.

Concluding Remarks
Our study of the virial theorem shows how extensively the variational method relies on cancellation in computation. For the orbital equation, the positive kinetic energy and the negative potential energy representing the interaction between an electron and the nucleus are the largest contributors to the total energy, yet they nearly cancel, particularly as the nuclear charge increases. All calculations for this study were performed in the double precision of about 15 significant digits. What is important is the accuracy of the sum of one-electron integrals, T + V 1 . These two quantities should possibly be computed in quadruple precision. THe deviation of the dilation from linear scaling for the pairs of orbitals also needed to be balanced between contributions from an inner and outer region. Thus, the simultaneous updates of orbitals are better able to balance virial theorem conditions that accurate solutions of the variational problem need to accomplish. The performance of these methods as implemented in an SPHF [18] program was better for the Newton-Raphson method when orbital rotation is included as part of the simultaneous update process [28]. However, most tests focused on light atoms; what are needed are stable and accurate multiconfiguration methods for systems with several open shells, such as some cases found in the lanthanides and actinides. The virial theorem could be a useful tool for confirming a stationary solution for a case with off-diagonal Lagrangian multipliers in both nonrelativistic and relativistic frameworks.
Lastly, Slater integrals in Cartesian coordinates contain a singularity in the integrand when r ij = 0, whereas no such singularity occurs in spherical coordinates. Integrating over spin-angular coordinates to determine an energy functional by using Equation (1) limits the value of k that defines the Slater integrals for a state [13].
Author Contributions: Methodology, both authors; software, C. Froese Fischer; validation, all authors; investigation, both authors; writing, review and editing, both authors. All authors have read and agreed to the published version of the manuscript. Note 1 2nd-order optimization, using both gradient and Hessian information, is implemented in the Quantum Chemistry relativistic package DIRAC [27] and activated in case the regular SCF does not converge, but is so far available for closed-shell calculations.