Next Article in Journal
High Sensitivity Multi-Axes Rotation Sensing Using Large Momentum Transfer Point Source Atom Interferometry
Next Article in Special Issue
A Tribute to Oleg Zatsarinny (1953–2021): His Life in Science
Previous Article in Journal / Special Issue
Scientific Activities of Oleg Zatsarinny in the Ukraine
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Towards B-Spline Atomic Structure Calculations

by
Charlotte Froese Fischer
Department of Computer Science, University of British Columbia, Vancouver, BC V6T1Z4, Canada
Submission received: 13 July 2021 / Revised: 26 July 2021 / Accepted: 27 July 2021 / Published: 31 July 2021

Abstract

:
The paper reviews the history of B-spline methods for atomic structure calculations for bound states. It highlights various aspects of the variational method, particularly with regard to the orthogonality requirements, the iterative self-consistent method, the eigenvalue problem, and the related sphf, dbsr-hf, and spmchf programs. B-splines facilitate the mapping of solutions from one grid to another. The following paper describes a two-stage approach where the goal of the first stage is to determine parameters of the problem, such as the range and approximate values of the orbitals, after which the level of accuracy is raised. Once convergence has been achieved the Virial Theorem, which is evaluated as a check for accuracy. For exact solutions, the V/T ratio for a non-relativistic calculation is −2.

1. Introduction

In 1996 Oleg Zatsarinny published a program for computing matrix elements in a non-orthogonal radial basis in which configuration state functions (CSFs) were expanded in determinants [1]. A non-orthogonal basis can greatly reduce the size of multiconfiguration wave function expansions. The orthogonal versus non-orthogonal issue is a complex one, so I invited Oleg for a three-month visit, with the National Science Foundation support for members of the former Soviet Union. Before he left, I shared with him the B-spline library I had developed for calculations of high accuracy “good to the last bit”. Oleg improved and extended the library, which became the basis for the BSR code he published in 2006 [2] for both non-relativistic and the Breit–Pauli R-matrix theory. It was not until 2011 [3] that I published a non-relativistic B-spline Hartree-Fock program for bound states using the extended library. In the mean time, Oleg extended BSR to Dirac relativistic theory where the R-matrix approach relies on eigenvectors of a B-spline matrix at the boundary. However with the traditional Galerkin method applied to the pair of first-order differential equations, spurious solutions appeared for the R-matrix.
The presence of spurious solutions from the application of the Galerkin methods is well known—a Google search for “Galerkin spurious” has 165,000 responses from areas of theoretical physics, applied mathematics, and engineering. Inspired by the work of Igarashi [4,5], I had the idea that the problem was related to the fact that Dirac equations are a system of first-order equations. I showed that the problem also occurs for y = λ 2 y , with boundary conditions y ( 0 ) = 0 and y ( R ) = 0 , when the single second-order equation is replaced by a pair of first-order equations. Igarashi was focusing mostly on kinetically-balanced basis sets but tried also some other ideas like using B-splines of a different order. If we think of the Galerkin method approximating the small component by say, Q k s ( r ) = i c i B i k s ( r ) , where B i k s ( r ) is a B-spline of order k s (or piecewise polynomial of degree k s 1 ), then the large component is proportional to ( d / d r κ / r ) Q k ( r ) , namely piecewise polynomials of degree k s 2 with a B-spline basis of order k s 1 . Representing these functions by a higher-order polynomial presents difficulties that would not arise if a lower order of splines were used. In my test case, the eigenvalues from using different orders for pairs of equations were the same as those from a single second-order equation, indicating that the problem was not just a Dirac equation-related problem. Splines of different orders worked beautifully for the R-matrix approach as well as the Thomas–Reiche–Kuhn sum rule and all our other tests. I convinced Oleg that we should prepare a paper for Physical Review Letters, but the referee was not as excited by this result as I was. Oleg was not prepared to argue so the paper was published in Computer Physics Communication [6].
The present paper is dedicated to the memory of Oleg Zatsarinny who extended the bound state B-spline codes to the fully relativistic dbsr-hf version [7].

2. In the Beginning ...

In the early days of computing, atomic and molecular calculations relied either on finite-difference methods based on numerical procedures developed by Hartree [8], or analytic methods based on basis sets such as Slater-type orbitals [9]. Bachau et al. [10] in their excellent review on B-splines, refer to these methods as local versus global methods, respectively, in that finite difference methods relied on only a few values at adjacent points of a grid in specifying an approximation. Much changed when de Boor [11] published his book in 1985 about B-splines with a local but complete basis for piece-wise polynomial approximations along with Fortran programs for standard spline procedures. Among the first to use his codes were Johnson and Sapirstein [12], who used B-splines to generate their numerical basis set for perturbation theory that greatly advanced the theory. However B-splines also offer advantages for a wide range of applications as Bachau et al., show in their extensive review.
By the 1990s, the non-relativistic Hartree–Fock variational method had evolved into a multiconfiguration theory as implemented in the mchf77 [13] program. Later this theory was extended to include relativistic corrections through a Breit–Pauli approximation and, together with programs for atomic properties (hyperfine, isotope shifts, and transition probabilities), became an atomic structure package atsp2K [14]. All calculations were based on numerical procedures and supported a limited amount of non-orthogonality of the radial basis [15] of one or at most two overlap integrals.
With the development of faster computers along with considerable more memory, an interesting question arose—what is the best method for solving the mchf equations? The latter required node counting to control the convergence of an integro-differential equation with many solutions, dealing also with Lagrange multipliers for satisfying orthonormality.
This paper reviews the development of spline methods for variational approaches to solutions of the wave equation for atoms with the assumption that the reader is familiar with the basic properties of splines as presented by Bachau et al. [10]. Unlike the latter, the emphasis here is more on the computational methods and the programs that have emerged rather than their application.

3. The B-Spline Basics

A spline approximation of order k s is an approximation that is a piecewise polynomial of degree k s 1 (with k s coefficients) in intervals defined by “knots” that define a grid. In applications, the grid itself, in a sense, is arbitrary but the choice of grid can have a large influence on the accuracy of an approximation.

3.1. Spline Grid for Radial Functions

In atomic structure applications, it is convenient to define the grid in terms of the variable t = Z r (the hydrogenic case) and then transform to the current value of Z. Let h = 1 / 2 m . This is not strictly necessary, but it ensures that the arithmetic is exact in binary arithmetic. Then:
t i = 0 for i = 1 , , k s
t i = t i 1 + h for i = k s + 1 , , k s + m
t i = t i 1 ( 1 + h ) for i = k s + m + 1 , , n s + 1
t i = t n s + 1 for i = n s + 2 , , n s + k s
r i = t i / Z for i = 1 , n s + k s
Thus the grid is linear for m intervals near the origin after the t = 0 knots of multiplicity k s , then linear in log ( r ) . If continuum calculations are involved, the range may be extended to include equally-spaced values at the end of the exponential range, followed by final knots of t = t max of multiplicity k s . The above grid is most appropriate for a finite nucleus. The number of non-zero intervals between knots is n s k s + 1 where n s is the number of basis states. Note also that, the end of the exponential region, t = t n s + 1 defines the maximum range R of the bound-state orbitals, that are assumed to be zero for r > R .
In an early B-spline study, Froese Fischer and Parpia [16] studied the accuracy of the application of B-splines to the Dirac equation for He with a grid over the range ( 0 , 40 ) , comparing a grid for ρ = r 1 / 4 with a grid for ρ = l o g ( r ) . The former was considerably more accurate but is rarely used other than for plotting since r = 0 transforms to ρ = 0 in the first transformation, whereas in the latter it transforms to . However, an advantage of splines is that points need not be equally spaced so r = 0 can be chosen as a grid-point (as in Equation (1)) even when the remaining grid-points are on an exponential grid.

3.2. Integration Methods

The piecewise polynomial values are represented by k s values at the Gaussian points for Gauss–Legendre integration. Then,
0 1 f ( x ) d x = i = 1 k s g w ( x i ) f ( x i ) ,
where g w ( x i ) are the Gaussian weights at the Gaussian points x i . This formula is exact for integrands f ( x ) that are polynomials of degree 2 k s 1 .

3.3. Slater Integrals

A detailed analysis of algorithms for the calculation of Slater integrals was reported by Qiu and Froese Fischer [17], and will be summarized here.
Let a , b , c , d denote a set of n l quantum numbers associated with a radial function:
P ( a ; r ) = i a i B i ( r )
with similar expansions for b , c , and d. Then,
R k ( a , b ; c , d ) = i j i j a i b j c i d j R k ( i , j ; i , j ) ,
where,
R k ( i , j ; i , j ) = 0 R 0 R r < k r > k + 1 B i ( r 1 ) B j ( r 2 ) B i ( r 1 ) B j ( r 2 ) d r 1 d r 2
will be referred to as B-spline Slater integrals. Many symmetries exist. Furthermore,
R k ( i , j ; i , j ) = 0 if | i i | > k s , o r = 0 if | j j | > k s .
Symmetry and the “local” property of these integrals, greatly reduce the number of B-spline Slater integrals included in the summation.
The spline-approximation divides the integration range into sub-intervals and the two-dimensional integration into patches or “cells”. Thus the fundamental process is one of integration over a cell, r i v r 1 r i v + 1 and r j v r 2 r i v + 1 , namely,
r i v r i v + 1 r j v r j v + 1 r < k r > k + 1 B i ( r 1 ) B j ( r 2 ) B i ( r 1 ) B j ( r 2 ) d r 1 d r 2 .
However for an off-diagonal cell, the above two-dimensional integral is separable and equal to (for i v < j v ):
r i v r i v + 1 B i ( r 1 ) r 1 k B i ( r 1 ) d r 1 r j v r j v + 1 B j ( r 2 ) ( 1 / r 2 k + 1 ) B j ( r 2 ) d r 2 r k ( i , i ; i v ) × r ( k + 1 ) ( j , j ; j v ) ,
where r k ( i , i ; i v ) and r ( k + 1 ) ( j , j ; j v ) are referred to as moments.
The diagonal cells reduce to integrations over upper and lower triangles or, through an interchange of arguments, as: 
R k ( i , j ; i , j ; i v ) = R Δ k ( i , j ; i , j ; i v ) + R Δ k ( j , i ; j , i ; i v )
and
R Δ k ( i , j ; i , j ; i v ) = r i v r i v + 1 1 r 1 k + 1 B i ( r 1 ) B i ( r 1 ) d r 1 r i v r 1 r 2 k B j ( r 2 ) B j ( r 2 ) d r 2 .
Thus the computational complexity of an integration over a diagonal cell is O ( k s 2 ) and for the R k ( i , j ; i , j ) array, it is O ( n s k s 2 ) . Furthermore, when all B-splines are defined on an exponential grid where t i + 1 = ( 1 + h ) t i , (as in the middle region of the logarithmic grid) scaling laws can be applied, namely:
B i + 1 ( r ) | r k | B j + 1 ( r ) = ( 1 + h ) 1 + k B i ( r ) | r k | B j ( r ) B i + 1 ( r ) | 1 / r k + 1 | B j + 1 ( r ) = ( 1 + h ) k B i ( r ) | 1 / r k + 1 | B j ( r ) R k ( i + 1 , j + 1 ; i + 1 , j + 1 ) = ( 1 + h ) R k ( i , j ; i , j ) .
With such a grid, some Slater integrals would be computed using the basic integration procedures and others could be scaled.

4. Tensor Products of B-Splines as a Basis

Traditionally, the wave functions for multi-electron systems are expanded in terms of configuration state functions (CSFs) that are products of one-electron wave functions (or orbitals) where, for bound-state problems, the latter provides an orthonormal basis. An early spline study of the 1 s 2 , 1 s 2 s 1 S or 3 S , and  1 s 2 p 1 P or 3 P states of Helium [18], took a different approach and explored the direct use of a tensor product of B-splines as a non-orthogonal basis.
For a two-electron system, a  wave function for a state Ψ ( γ L S ) , where γ denotes the configuration, can be expressed in terms of pair functions, identified by their angular symmetry, namely:
Ψ ( γ L S ) = n n c n n n l n l L S .
Here,
n l n l L S = N 1 r 1 r 2 ( 1 P 12 ) P ( n l ; r 1 ) P ( n l ; r 2 ) l l L S ,
where N is a normalisation factor that may depend on symmetry, P 12 is a permutation operator, l l L S a is spin-angular factor that identifies the pair function, and P ( n l ; r ) and P ( n l ; r ) are radial functions. In other words, the pair function is a linear combination of configuration state functions of the same symmetry. Substituting into Equation (16) and noting that each term has the same spin-angular factor, that the double sum is an expansion in a basis of a function of two variables, we get a pair function p ( r 1 , r 2 ) , such that:
Ψ ( l l L S ) = N 1 r 1 r 2 ( 1 P 1 , 2 ) p l ( r 1 , r 2 ) l l L S .
In essence, p l ( r 1 , r 2 ) is a general two-dimensional function, usually expanded in terms of orbitals but one that could also be expanded in terms of tensor products of B-splines, namely B i ( r 1 ) B j ( r 2 ) . The advantage of this basis is the local support (non-zero region) compared with products of orbitals which extend over the entire two-dimensional region.
In the case of 1 P state, the possible orbital symmetries ( s p , p d , d f , f g ) can be represented as ( l , l ) where l = 0 , 1 , l m a x , with  l = l + 1 . Then the 1 s 2 p 1 P state can be expressed in terms of partial waves as:
Ψ ( 1 s 2 p 1 P ) = 1 2 r 1 r 2 l ( 1 P 1 , 2 ) p l ( r 1 , r 2 ) l l 1 P , l = l + 1 ,
where,
p l ( r 1 , r 2 ) = i j p i , j l B i ( r 1 ) B j ( r 2 ) .
The application of the Galerkin method leads to a system linear equations, where the equations for a specific basis, identified by the parameters ( i , j , l ) , 1 i , j N and l = 0 , , l m a x , are:
i j p i j l L i i l S j j + S i i L j j l E S i i S j j + l i j k c ( l l , l l ; k 1 ) p i j l R k ( i , j ; i i ) = 0 .
In the above, L l = ( L i i l ) and S = ( S i i ) are matrices with:
L i i l = B i ( r ) | 1 2 d 2 d r 2 Z r + l ( l + 1 ) 2 r 2 | B i ( r ) , S i i = B i ( r ) | B i ( r ) .
In the B-spline basis, the one-electron operators ( L l and S matrices) connect the coefficients within a pair function ( p l ( i , j ) for a given l) whereas the Slater integrals in the B-spline basis, connect the different pair functions. The coefficients of the Slater integrals c ( l l , l l ; k L ) arise from angular integrals related to the Slater integrals. Both can be treated as data for the calculation and depend only on the underlying grid. The equations themselves were solved iteratively with the energy E determined as a Rayleigh quotient. The solution is numerically intensive and has many options for parallelism.
Of interest in 1991 was the performance of the algorithm on parallel vector processors. An accuracy of a micro-Hartree was achieved for the energy. Figure 1 shows the matrix of coefficients p i , j l for B-spline tensor-product expansions for s p , p d , d f , f g (or l = 0 , 1 , 2 , 3 ) pair functions from left to right and top to bottom, respectively. Notice that each subplot has a different scale with the maximum value decreasing. For  s p symmetry, the matrix is approximately the cross-product of the expansion coefficients for the 1 s and 2 p Hartree–Fock orbitals in a B-spline basis. As l increases, the maximum coefficient (dark red in colour) moves closer to the diagonal r 1 = r 2 region. However more importantly, the significant components are concentrated in a smaller and smaller region. What this calculation clearly shows is that correlation is a “local” correction and that, as the orbital symmetry l increases, an oscillating orthonormal basis would not be an efficient basis.
However the example is also interesting from other perspectives. Like the Hylleraas method, there are no orbitals, nor are there orthonormality constraints and hence, no Lagrange multipliers. The wave function is a sum of pair functions rather than a linear combination of configuration state functions, and there is no matrix diagonalisation. In present grasp [20] or atsp calculations, the usual first step is to compute all angular data, which then needs to be read and stored in memory when needed. With pair functions the amount of angular data is greatly reduced and often could be computed as needed. It should also be remembered that orbitals are a theoretical concept.

5. Spline Galerkin and Inverse Iteration Methods

The development of B-splines (of an order greater than that of cubic splines) began before the LAPACK [21] routines were released. Available instead, were LINPACK [22] routines for solving systems of equations. Some early papers by Froese Fischer and Idrees [23,24], described a spline algorithm for solving the continuum functions that used the spline Galerkin method for deriving linear equations of the form:
( H E S ) c = A ( E ) c = 0
and inverse iteration for solving the equation for a given energy. The latter was similar to the power method for finding the eigenvector associated with the largest eigenvalue. For continuum solutions, the energy E is specified and what is needed is the eigenvector of the nearest (or smallest) eigenvalue to E. Like the power method, the method is iterative.
Let the matrix A be an N × N matrix and L and U lower and upper triangular matrices, respectively, of similar dimension, such that A = L U . and c i ( 0 ) = 1 , i = 1 , , N is a vector. Then, starting with m = 0 and incrementing by 1, until the vector c has converged (i.e., | c i m + 1 c i ( m ) | < 10 12 ), let:
L y = c i ( m ) U x = y c ( m + 1 ) = x / | | x | | ,
where x and y are vectors of length N. The method was tested for the hydrogen scattering problem and then photoionisation in He. It was also determined how orthogonality conditions could be included by extending the definition of A ( E ) and the vector c for the case,
Ψ ( 1 s k s 1 S ) = c 0 | 1 s 2 > + | 1 s k s > ,
as an example. Resonances, phase shifts, or photoionisation cross-sections were extracted from the results. A  more extensive investigation of resonance positions and widths for H and He was reported by Brage et al. [25,26], adding to the rich ‘flora’ of results by many different methods for these cases. Xi and Froese Fischer extended these results to three-electron He system in the investigation of cross-section and angular-distribution for the photodetachment of He 1 s 2 s 2 p 4 P o below the He ( n = 4 ) threshold [27] and also below the 1 s detachment level [28]. A multichannel theory was developed. Among the resonances found was a 2 s 2 p 2 4 P resonance state immediately below the 1 s threshold. A similar theory was applied to Be   1 s 2 2 s 2 p 2 4 P  [29].
The above methods had only one region. These simple approaches were extended by Oleg Zatsarinny to numerous continuum processes based on the R-matrix method [2,30] with its inner and outer regions and non-orthogonal orbitals. He referred to one region methods as “straight forward” methods.

6. Spline Methods for Bound State Problems

Several options are available for B-spline solutions that were not feasible for finite difference methods. Consider the simple equation for 1 s 2 that was studied by Froese Fischer and Guo [31]. The differential equation that needs to be solved is:
d 2 d r 2 2 r Z Y 0 ( 1 s , 1 s ; r ) ε P ( 1 s ; r ) = 0 ,
where P ( 1 s ; r ) = 0 when r = 0 or r and the radial function P ( 1 s ; r ) is expanded in a B-spline basis so that:
P ( 1 s ; r ) = c i B i ( r ) .
This problem can be linearised by computing Y 0 ( 1 s , 1 s ; r ) from current estimates, and then solved as a generalised eigenvalue problem for an improved estimate, with at best a linear rate of convergence. Or, we can think of the equations as non-linear equations, and solve for changes in the expansion coefficients and energy parameters, using the Newton–Raphson (NR) method with a quadratic rate of convergence (see Ref. [32] for details).
The many-electron variational methods are extensions of the Hartree–Fock methods [3,19] for which three categories of methods have been implemented and evaluated.

6.1. Generalised Eigenvalue Problem for a Single Orbital

In this approach, orbitals are improved one at a time according to a generalized eigenvalue problem,
( H a ε a a S ) a = 0 .
When two orbitals are constrained through orthogonality, as in the case of 1 s 2 s 1 S , then the iterative process will not converge without first rotating the two orbitals, say a and b for a stationary energy, a process implemented in the sphf program [3]. Projection operators may then be applied to eliminate the off-diagonal Lagrange multiplier ε a b . The matrix H a is then a full matrix when exchange contributions are present.

6.2. Multiple Orbitals and SVD

Another possibility is to solve for a set of orbitals at the same time using singular value decomposition (SVD) which is closely related to inverse iteration and is included in LAPACK [21].
Let A i be the expansion vector for orbital i. Then the system of equations can be written as:
F 11 F 12 F 1 m F 21 F 22 F 2 m F m 1 F m 2 F m m A 1 A 2 A m = 0 ,
where F i i contains the contributions from one-electron integrals, the direct Slater integrals for orbital i, as well as ε i i B , and  F i j ( i j ) contains the contribution from exchange integrals between orbitals i and j and possible orthogonality constraints. In this case, all matrices are banded. In this approach, the energy parameters are computed from current estimates of the orbitals.
The earlier SVD study [19] showed poor convergence in the case of 1 s 2 s 1 S but 1 s 2 2 s 2 converged linearly. For the latter, the off-diagonal energy parameter is zero and rotation of orbitals is not important. Applying a projection operator to a system of equations is equivalent to using current estimates of a radial function to determine off-diagonal energy parameters. It is possible that SVD equations should be extended to include the energy parameters (diagonal ε i i ) as well as off-diagonal ( ε i j ) unknowns so that an effective rotation of orbitals is part of the SVD solution.

6.3. Newton–Raphson with Quadratic Rate of Convergence

Consider the case of two orbitals a and b, with an orthogonality constraint between them. Then the unknowns are ( a , b , ε a a , ε b b , ε a b ) . The equations to be solved are the two orbital equations along with the three orthonormality conditions that are part of the energy functional. Let ( a , b ) be the current estimates that are used to evaluate the ε -matrix. If symmetry conditions are not satisfied exactly, we can use the average value,
ε a b = ε b a = b t H a a + a t H b b / 2 .
Then, by the Newton–Raphson method [3]. Δ a and Δ b are solutions of:
H a a ε a a B H a b ε a b B B a B b H b a ε a b B H b b ε b b B B b B a ( B a ) t ( B b ) t ( B b ) t ( B a ) t Δ a Δ b Δ ε a a Δ ε b b Δ ε a b = res a res b 0 0 0 .
In the above, H a ε a a S = res a , namely the amount by which the current estimates do not satisfy the equations.

7. The sphf and spmchf Programs

The B-spline methods clearly offer some advantages, even when they are more computationally intensive than finite difference methods. The eigenvalue method is the preferred method for singly occupied orbitals, particularly with large principal quantum numbers. However for a multiply occupied shell, the Newton–Raphson method is efficient in that the “self-energy” is readily accounted for which greatly improves the rate of convergence. When many subshells are present with different principal quantum numbers, convergence can be improved with the simultaneous improvement of orbitals. For this reason, both sphf [3] and spmchf (available at GitHub [33]) have “phases” and different “levels” of accuracy. All iterative methods need initial estimates. The first phase is stable even when estimates are of a poor accuracy and the grid that defines the B-spline basis is relatively coarse. Orbitals are updated sequentially using the generalised eigenvalue method on the first iteration but, in other iterations use NR for a single orbital when a subshell is multiply occupied. This first phase is referred to as the “SCF” phase. The codes also have the initial level of accuracy (with a coarse grid) and less accurate tests for convergence and a higher level of accuracy. When initial estimates are for the refined grid, the program assumes there is only one level of accuracy. In any event, the convergence tests may be reset by the user, over-riding program defaults. Unlike the sphf program, the spmchf program does not include orbital rotation in this phase. The next phase is thought of as a “clean-up” phase that raises the level of accuracy and improves the relationship of one subshell to the other. In the Be 1 s 2 2 s 2 case, if the 1 s orbital is too contracted the 2 s subshell will be too expanded. The NR phase deals with this issue that greatly affects the Virial Theorem (V/T) ratio, as well as orbital rotation when off-diagonal energy parameters are important.
Figure 2 shows the parameters for two levels in the calculation for Be 1 s 2 2 s 2 1 S . The default grids and convergence parameters are displayed for both accuracy levels. No initial estimates were provided in this case so the initial range was estimated to be fairly large. At the second level of accuracy, the range is known much more precisely. Notice the improved VT in going from the first to the second level of accuracy. The  1 s 2 s 1 S case is very similar except now off-diagonal Lagrange multipliers are needed for a stationary energy with respect to rotation and the requirement that ε a b = ε b a . Table 1 shows how during the SCF iterations, the values are opposite in sign and the calculations do not converge, but everything changes during the NR iterations and iterations converge rapidly. The default options of dbsr-hf do not include orbital rotations. For the 1 s 2 s 1 S case, it converges but to an incorrect value, made evident only through the Virial theorem (VT). Thus VT is an important check for this code.
The spmchf method differs not only in that expansion coefficients for configuration states (CSFs) need to be determined but also in that the the energy expression is no longer limited to direct ( F k ( a , b ) ) and exchange ( G k ( a , b ) ) Slater integrals. The process by which the energy expression relates to the matrix form a H a a requires that two occurrences of the orbital a be present in the matrix element [19]. In the MCHF approximation:
Ψ ( 3 s 3 d 1 D ) = c 1 Φ ( 3 s 3 d 1 D ) + c 2 Φ ( 3 p 2 1 D )
the interaction matrix element is 2 15 R 1 ( 3 s 3 d , 3 p 3 p ) . The spmchf treats the contribution to Equation (27) for 3 s or 3 d as a “residual” term (not included in the matrix), but Zatsarinny [34] pointed out that if the matrix element was treated as:
2 15 R 1 ( 3 s 3 d , 3 p 3 p ) S ( 3 s , 3 s ) S ( 3 d , 3 d )
then the matrix form could be retained. The factors, S ( 3 s , 3 s ) and S ( 3 d , 3 d ) are normalisation integrals. This modification has not yet been implemented in spmchf.
Figure 3 and Figure 4 show the spmchf solution for this example. Neither sphf nor spmchf require initial estimates so the first calculation was for a Hartree Fock calculation for Mg 3 s 3 d 1 D . The output was then used in the second multiconfiguration run as input to serve as initial estimates and the first level varied only the 3 p orbital, not present in the input. Then the second varied all orbitals with excellent results as reported.
At this time, the spmchf program for multiconfiguration wave functions has not yet been been fully tested. Through the use of the term=jj command-line option, the  fully relativistic dbsr-hf code can be used as an average energy level (EAL) approximation.

8. Concluding Remarks

The spline codes reviewed in this article were developed in the last 20 or so years. For bound state problems, the GRASP code, for example [20], based on finite difference methods and an orthonormal orbital basis, is extensively used for complete spectra including relatively highly excited levels [35] and certain heavy elements, primarily for highly ionised systems where valence correlation can be dealt with and the effect of core correlation on spectra and other properties is limited. spmchf with its greater flexibility still needs to be tested on lanthanides and actinides with two open shells ( n = 4 , 5 for lanthanides and n = 5 , 6 for actinides) where natural orbital transformation could prove useful in connection with orthonormal radial functions [36]. Multiple f-orbitals could also be a challenge for bsr in that the number of determinants can increase rapidly.
In 1984, it took 254 seconds to execute a numerical HF program for Ra 7 s 2 1 S on a VAX 11/780; in 1987, 2 seconds on a Cray X-MP; and in 1993, 1 s on a Dec Alpha, which was a popular computer at its time. Thus time (cost) is no longer an important factor but rather ease of use and reliability. The same is not true when accurate wave functions are needed for many-electron systems for different atomic properties. Important now is how well an algorithm can be made parallel, and how efficiently memory usage can be managed.
The present paper has been about non-relativistic calculations. Relativistic methods are similar, except that the decision to represent large and small components of radial functions by splines of different order, automatically implies that every Slater integral is the sum of four integrals. This decision was important in bsr [2] in that it eliminated spurious solutions in the R-matrix and was also used for dbsr-hf. However, is it the most efficient solution of the problem? For a grid of n v intervals, the number of independent basis functions is n v + k s 1 . Thus, by going from k s to k s 1 , the number of independent basis states decreases. The spurious solutions generally are higher in the spectrum. For bound state solutions, the asymptotic conditions are such that both the large and small components and their derivatives go to zero. In fact, in non-relativistic calculations, applying both conditions stabilises the solution at large r. Further study might be appropriate.

Funding

This research was funded by Canada NSERC Discovery Grant 2017-03851.

Data Availability Statement

Not applicable.

Acknowledgments

Much of the spline research described here was performed over many years at Vanderbilt University with funding support from the Division of Chemical Science, Geosciences, and Bioscience, Office of Basic Energy Sciences, Office of Science, and the U.S. Department of Energy. Some of the program development was performed while at the Atomic Physics Division, National Institute of Standards and Technology, Gaithersburg, and Maryland 20899-8422. The author very much appreciated Professor Michel R. Godefroid’s comments and suggestions for this paper.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Zatsarinny, O. A general program for computing matrix elements for atomic structure with nonorthogonal orbitals. Comp. Phys. Commun. 1996, 98, 235–254. [Google Scholar] [CrossRef]
  2. Zatsarinny, O. BSR: B-spline atomic R-matrix codes. Comput. Phys. Commun. 2006, 174, 273–356. [Google Scholar] [CrossRef]
  3. Froese Fischer, C. A B-spline Hartree-Fock program. Comput. Phys. Commun. 2011, 182, 1315–1326. [Google Scholar]
  4. Igarashi, A. B-Spline Expansions in Radial Dirac Equation. J. Phys. Soc. Jpn. 2006, 75, 114301. [Google Scholar] [CrossRef]
  5. Igarashi, A. Kinetically balanced B-spline expansions in radial Dirac equation. J. Phys. Soc. Jpn. 2007, 76, 05431. [Google Scholar] [CrossRef]
  6. Froese Fischer, C.; Zatsarinny, O. A B-spline Galerkin method for the Dirac equation. Comput. Phys. Commun. 2009, 180, 879–886. [Google Scholar] [CrossRef] [Green Version]
  7. Zatsarinny, O.; Froese Fischer, C. DBSR—A B-spline Dirac-Hartree-Fock program. Comput. Phys. Commun. 2016, 202, 287–303. [Google Scholar] [CrossRef] [Green Version]
  8. Hartree, D.R. The Calculations of Atomic Structures; Springer: New York, NY, USA, 1957. [Google Scholar]
  9. Roothaan, C.C.J. New Developments in Molecular Orbital Theory. Rev. Mod. Phys. 1951, 23, 69–89. [Google Scholar] [CrossRef]
  10. Bachau, H.; Cormier, E.; Decleva, P.; Hansen, J.E.; Martin, F. Applications of B-splines in atomic and molecular physics. Rep. Prog. Phys. 2001, 64, 1815–1942. [Google Scholar] [CrossRef]
  11. de Boor, C. A Practice Guide to Splines; Springer: New York, NY, USA, 1985. [Google Scholar]
  12. Johnson, W.R.; Sapirstein, J. Computation of Second-Order Many-Body Corrections in Relativistic Atomic Systems. Phys. Rev. Lett. 1986, 57, 1126. [Google Scholar] [CrossRef]
  13. Froese Fischer, C. A general multiconfiguration Hartree-Fock program. Comput. Phys. Commun. 1978, 14, 145–153. [Google Scholar] [CrossRef]
  14. Froese Fischer, C.; Tachiev, G.; Gaigalas, G.; Godefroid, M.R. An MCHF atomic-structure package for large-scale calculations. Comput. Phys. Commun. 2007, 176, 559–579. [Google Scholar] [CrossRef] [Green Version]
  15. Hibbert, A.; Froese Fischer, C.; Godefroid, M.R. Non-orthogonal orbitals in MCHF or configuration interaction wave functions. Comput. Phys. Commun. 1988, 51, 282–293. [Google Scholar] [CrossRef]
  16. Froese Fischer, C.; Parpia, F.A. Accurate spline solutions of the radial Dirac equation. Phys. Lett. A 1993, 179, 198–204. [Google Scholar] [CrossRef]
  17. Qiu, Y.; Froese Fischer, C. Integration by cell algorithm for Slater integrals in a spline basis. J. Comput. Phys. 1999, 156, 257–271. [Google Scholar] [CrossRef]
  18. Froese Fischer, C. Concurrent vector algorithms for spline solutions of the helium pair equation. Int. J. High Perform. Comput. Appl. 1991, 5, 5–20. [Google Scholar] [CrossRef]
  19. Froese Fischer, C. B-splines in variational atomic structur calculations. Adv. At. Mol. Phys. 2007, 55, 539–550. [Google Scholar]
  20. Froese Fischer, C.; Gaigalas, G.; Jönsson, P.; Bieroń, J. GRASP2018—A Fortran 95 version of the General Relativistic Atomic Structure Package. Comput. Phys. Commun. 2019, 237, 184–187. [Google Scholar] [CrossRef]
  21. LAPACK Library. Available online: http://www.netlib.org/lapack/ (accessed on 28 July 2021).
  22. LINPACK Library. Available online: http://www.netlib.org/linpack/ (accessed on 28 July 2021).
  23. Froese Fischer, C.; Idrees, M. Spline algorithms for continuum functions. Comput. Phys. 1989, 3, 53–58. [Google Scholar] [CrossRef]
  24. Froese Fischer, C.; Idrees, M. Spline methods for resonances in photoionisation cross sections. J. Phys. B At. Mol. Opt. Phys. 1990, 23, 679. [Google Scholar] [CrossRef]
  25. Brage, T.; Froese Fischer, C.; Miecznik, G. Non-variational, spline-Galerkin calculations of resonance positions and widths, and photodetachment and photoionization cross sections for H and He. J. Phys. B At. Mol. Opt. Phys. 1992, 25, 5289–5314. [Google Scholar] [CrossRef]
  26. Brage, T.; Froese Fischer, C. Spline-Galerkin methods for Rydberg series, including Breit-Pauli effects. J. Phys. B At. Mol. Opt. Phys. 1994, 27, 5467–5484. [Google Scholar] [CrossRef]
  27. Xi, J.; Froese Fischer, C. Cross section and angular distribution for the photodetachment of He1s2s2p4Po below the He n = 4 threshold. Phys. Rev. A 1996, 53, 3169–3177. [Google Scholar] [CrossRef]
  28. Xi, J.; Froese Fischer, C. Photodetachment cross-section of He- (1s2s2p4Po ) in the region of the 1s detachment threshold. Phys. Rev. A 1999, 59, 307–314. [Google Scholar] [CrossRef]
  29. Xi, J.; Froese Fischer, C. Cross section and angular distribution for photodetachment of Be− 1s22s2p24P. J. Phys. B At. Mol. Opt. Phys. 1999, 32, 387–396. [Google Scholar] [CrossRef]
  30. Zatsarinny, O.; Froese Fischer, C. The use of basis splines and non-orthogonal orbitals in R-matrix calculations: Application to Li photoionization. J. Phys. B At. Mol. Opt. Phys. 2000, 33, 313–341. [Google Scholar] [CrossRef]
  31. Froese Fischer, C.; Guo, W. Spline algorithms for the Hartree-Fock equation for the helium ground state. J. Comput. Phys. 1990, 90, 486–496. [Google Scholar] [CrossRef]
  32. Froese Fischer, C.; Guo, W.; Shen, Z. Spline methods for multiconfiguration Hartree-Fock calculations. Int. J. Quantum Chem. 1992, 42, 849–867. [Google Scholar] [CrossRef]
  33. spmchf. Available online: http://github.com/compas/spmchf (accessed on 28 July 2021).
  34. Zatsarinny, O.; Drake Umiversity. Private communication, 1 May 2016.
  35. Li, W.; Amarsi, A.M.; Papoulia, A.; Ekman, J.; Jönsson, P. Extended theoretical transition data in C i–iv. Mon. Not. R. Astron. Soc. 2021, 502, 3780–3799. [Google Scholar] [CrossRef]
  36. Schiffmann, S.; Godefroid, M.; Ekman, J.; Jönsson, P.; Froese Fischer, C. Natural orbitals in multiconfiguration calculations of hyperfine-structure parameters. Phys. Rev. A 2020, 101, 062510. [Google Scholar] [CrossRef]
Figure 1. A visualisation of the magnitude of the matrix of expansion coefficients p i , j l { i , j } = 1 , , 30 for the 1 s 2 p 1 P state of Helium. Shown, from left to right and top to bottom, are the expansion coefficients for the s p , p d , d f , and f g pair functions. The maximum values are 0.45, 0.010, 0.0035, 0.0020, respectively. Note both the changing region and decreasing maximum magnitude of each partial wave. (See also [19].)
Figure 1. A visualisation of the magnitude of the matrix of expansion coefficients p i , j l { i , j } = 1 , , 30 for the 1 s 2 p 1 P state of Helium. Shown, from left to right and top to bottom, are the expansion coefficients for the s p , p d , d f , and f g pair functions. The maximum values are 0.45, 0.010, 0.0035, 0.0020, respectively. Note both the changing region and decreasing maximum magnitude of each partial wave. (See also [19].)
Atoms 09 00050 g001
Figure 2. Computer output showing the convergence for the two levels of accuracy. The calculations are for Be 1 s 2 2 s 2 1 S .
Figure 2. Computer output showing the convergence for the two levels of accuracy. The calculations are for Be 1 s 2 2 s 2 1 S .
Atoms 09 00050 g002
Figure 3. Computer output of a Hartree–Fock calculation for the Mg [Ne] 3 s 3 d 1 D wave function.
Figure 3. Computer output of a Hartree–Fock calculation for the Mg [Ne] 3 s 3 d 1 D wave function.
Atoms 09 00050 g003
Figure 4. Computer output for multiconfiguration calculation for Ψ ( 3 s 3 d 1 D ) = c 1 Φ ( 3 s 3 d 1 D ) + c 2 Φ ( 3 p 2 1 D ) . Calculations are for Mg [Ne] 3 s 3 d 1 D using the radial functions from Figure 3 as initial estimates.
Figure 4. Computer output for multiconfiguration calculation for Ψ ( 3 s 3 d 1 D ) = c 1 Φ ( 3 s 3 d 1 D ) + c 2 Φ ( 3 p 2 1 D ) . Calculations are for Mg [Ne] 3 s 3 d 1 D using the radial functions from Figure 3 as initial estimates.
Atoms 09 00050 g004
Table 1. Table showing the convergence of off-diagonal energy parameters ε 1 s 2 s and ε 2 s 1 s as a function of the iteration (i) for the eigenvalue SCF iteration updating orbitals sequentially and the Newton–Raphson (NR) method updating both simultaneously. The calculations are for He 1 s 2 s 1 S . The values for i = 0 of the SCF iteration are computed from initial estimates. For all other iterations, values are determined after the iteration has completed. For exact solutions, ε 1 s 2 s = ε 2 s 1 s . Note that the SCF process did not converge.
Table 1. Table showing the convergence of off-diagonal energy parameters ε 1 s 2 s and ε 2 s 1 s as a function of the iteration (i) for the eigenvalue SCF iteration updating orbitals sequentially and the Newton–Raphson (NR) method updating both simultaneously. The calculations are for He 1 s 2 s 1 S . The values for i = 0 of the SCF iteration are computed from initial estimates. For all other iterations, values are determined after the iteration has completed. For exact solutions, ε 1 s 2 s = ε 2 s 1 s . Note that the SCF process did not converge.
SCF NR
i ε 1 s 2 s ε 2 s 1 s i ε 1 s 2 s ε 2 s 1 s
00.034327680.3574208340.150389170.14949971
1−0.006963790.1350364650.150933340.15089696
2−0.009694130.1410638460.150896960.15089696
3−0.009687230.1410418170.150896910.15089691
80.150897060.15089705
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Fischer, C.F. Towards B-Spline Atomic Structure Calculations. Atoms 2021, 9, 50. https://0-doi-org.brum.beds.ac.uk/10.3390/atoms9030050

AMA Style

Fischer CF. Towards B-Spline Atomic Structure Calculations. Atoms. 2021; 9(3):50. https://0-doi-org.brum.beds.ac.uk/10.3390/atoms9030050

Chicago/Turabian Style

Fischer, Charlotte Froese. 2021. "Towards B-Spline Atomic Structure Calculations" Atoms 9, no. 3: 50. https://0-doi-org.brum.beds.ac.uk/10.3390/atoms9030050

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop