Computer Physics Communications 174 (2006) 273–356

[Link]/locate/cpc

BSR: B-spline atomic R-matrix codes ✩

Oleg Zatsarinny
Department of Physics and Astronomy, Drake University, Des Moines, IA 50311, USA
Received 19 April 2005; received in revised form 26 October 2005; accepted 30 October 2005
Available online 20 December 2005

Abstract
BSR is a general program to calculate atomic continuum processes using the B-spline R-matrix method, including electron–atom and electron–
ion scattering, and radiative processes such as bound–bound transitions, photoionization and polarizabilities. The calculations can be performed
in LS-coupling or in an intermediate-coupling scheme by including terms of the Breit–Pauli Hamiltonian.

New version program summary

Title of program: BSR

Catalogue identifier: ADWY
Program summary URL: [Link]
Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland
Computers on which the program has been tested: Microway Beowulf cluster; Compaq Beowulf cluster; DEC Alpha workstation; DELL PC
Operating systems under which the new version has been tested: UNIX, Windows XP
Programming language used: FORTRAN 95
Memory required to execute with typical data: Typically 256–512 Mwords. Since all the principal dimensions are allocatable, the available
memory defines the maximum complexity of the problem
No. of bits in a word: 8
No. of processors used: 1
Has the code been vectorized or parallelized?: no
No. of lines in distributed program, including test data, etc.: 69 943
No. of bytes in distributed program, including test data, etc.: 746 450
Peripherals used: scratch disk store; permanent disk store
Distribution format: [Link]
Nature of physical problem: This program uses the R-matrix method to calculate electron–atom and electron–ion collision processes, with options
to calculate radiative data, photoionization, etc. The calculations can be performed in LS-coupling or in an intermediate-coupling scheme, with
options to include Breit–Pauli terms in the Hamiltonian.
Method of solution: The R-matrix method is used [P.G. Burke, K.A. Berrington, Atomic and Molecular Processes: An R-Matrix Approach, IOP
Publishing, Bristol, 1993; P.G. Burke, W.D. Robb, Adv. At. Mol. Phys. 11 (1975) 143; K.A. Berrington, W.B. Eissner, P.H. Norrington, Comput.
Phys. Comm. 92 (1995) 290].
 2005 Elsevier B.V. All rights reserved.

PACS: [Link]; [Link]; [Link]; [Link]

Keywords: Electron–atom scattering; Electron–ion scattering; Photoionization; Oscillator strengths; Polarizability; R-matrix; B-splines

✩
This paper and its associated computer program are available via the Computer Physics Communications homepage on ScienceDirect ([Link]
com/science/journal/00104655).
E-mail address: oleg_zoi@[Link] (O. Zatsarinny).

0010-4655/$ – see front matter  2005 Elsevier B.V. All rights reserved.
doi:10.1016/[Link].2005.10.006
274 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

Contents

1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 275
2. General theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 277
2.1. The close-coupling expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 277
2.2. The R-matrix method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 279
2.3. The external region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 281
2.4. Radiative processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282
2.5. Bound-state calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283
2.6. Relativistic corrections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284
3. General features of the BSR complex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 286
3.1. Preparation of target states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 286
3.2. Main programs and data files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 286
3.3. Others features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 286
4. Program BSR_PREP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 287
4.1. Outline of the BSR_PREP calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 287
4.2. Data files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 288
4.3. Input data in file target . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 289
5. Program BSR_CONF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 289
5.1. Outline of the BSR_CONF calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 289
5.2. Data files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292
5.3. Input parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293
5.4. Modification of file target . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293
6. Program BSR_BREIT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 294
6.1. Outline of the BSR_BREIT calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 294
6.2. Structure and data flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 295
6.3. Data files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 297
6.4. Input parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 297
6.5. Structure of the angular-coefficients databank . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 297
7. Program BSR_MAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 298
7.1. Outline of the BSR_MAT calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 298
7.2. Structure and data flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 301
7.3. Data files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 302
7.4. Input data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303
7.5. Output of Hamiltonian matrices and asymptotic coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304
8. Program BSR_HD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 305
8.1. Structure and data flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 305
8.2. Data files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 307
8.3. Input data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 308
8.4. Structure of the H file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 309
8.5. Output of the R-matrix solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 309
8.6. Output of bound solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 310
8.7. Output of weights . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 310
9. Program MULT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 310
9.1. Outline of the MULT calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 310
9.2. Structure and data flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313
9.3. Data files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 315
9.4. Input parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 315
9.5. Structure of the dipole databank . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 315
10. Program BSR_DMAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 316
10.1. Dipole transition matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 316
10.2. Oscillator strengths and decay probabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 318
10.3. Structure and data flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 318
10.4. Data files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 320
10.5. Input parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 321
10.6. Output of dipole matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 321
10.7. Output of radiative data in zf_res . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 321
11. The program BSR_PHOT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 322
11.1. Structure and data flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 322
11.2. Data files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323
11.3. Input parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324
12. Libraries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 275

12.1. Description of B-splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 325

12.2. General routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 327
12.3. Atomic orbital routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 331
12.4. Evaluation of two-electron integrals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332
12.5. The library ZCONF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343
12.6. The library ZCOM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 346
13. Additional programs-utilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 347
14. Compilation and test runs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 349
14.1. Compilation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 349
14.2. Test runs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 349
14.3. LS-coupling test run . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 350
14.4. jK-coupling test run . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 351
14.5. Photoionization calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352
14.6. Bound-state calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 355
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 355

1. Introduction

Collisions of electrons and photons with atoms and ions are among the most elementary processes that can be studied within
the framework of quantum mechanics. The latter processes include elastic scattering, inelastic excitation, and ultimately colli-
sional break-up in ionization. These collisions determine to a considerable extent the energy balance for different plasma and laser
devices, while the understanding of the underlying mechanisms allows one to study comprehensively the properties of atomic sys-
tems and their structure. Benchmark experiments and calculations have further advanced our basic knowledge of this field, while
at the same time providing urgently needed input data for modeling programs describing astrophysical and laboratory plasmas,
planetary atmospheres, and lasers. Due to the many difficulties associated with experimental investigations of these collisions,
such as the normalization of absolute cross sections and/or low count rates, and also due to the enormous amount of data needed
for the modeling programs, the vast majority of currently available input data for these programs consists of theoretical predic-
tions.
During the past decade, a number of general computer codes have been developed to generate accurate data for the various scat-
tering parameters. The R-matrix technique describing electron and photon scattering in complex atoms and ions was introduced
to atomic collision physics by Burke and collaborators [1,2] and now is the method most frequently applied for the determination
of various atomic properties. The compilation volume by Burke and Berrington [3] includes the major papers in R-matrix the-
ory until a decade ago, along with a comprehensive list of calculations using the method. Two of the major world-wide projects
successfully undertaken with extensive use of the general R-matrix package published by Berrington et al. [4] are the Opac-
ity Project ([Link] and the Iron Project ([Link]
[Link]). The above suite of codes, RMATRX-I, still uses the basic ansatz of the version first published in 1974 [5]. While
many subroutines were rewritten to increase efficiency and user-friendliness [6] and to include relativistic effects [7], these were
essentially updates to the original version. More radical changes were made in RMATRX-II, with particular emphasis placed
on improving the efficiency of the angular-momentum package. The latter code is now parallelized and resulted in the recent
publication of PRMAT [8]. Straightforward applications of this program, however, are more limited than those of RMATRX-I;
currently only electron collisions (no photon processes) can be handled within a non-relativistic framework. Consequently, work
on producing a robust parallel version of RMATRX-I is still in progress ([Link]
Relativistic effects can be included in RMATRX-I at the level of first-order perturbation theory by adding the one-body rela-
tivistic terms of the Breit–Pauli Hamiltonian (essentially the Pauli approximation) to the non-relativistic Hamiltonian that needs
to be diagonalized. Although fully relativistic versions of both atomic structure and collision packages exist (GRASP [9] and
DARC [10], respectively), applications of DARC have been comparatively rare and mostly restricted to highly ionized sys-
tems.
There are also other packages around, such as the published code RMATRX-ION [11] to handle general ionization as well
as ionization–excitation processes, another modified version of RMATRX-I ([Link] and the
intermediate-energy R-matrix code IERM [12] based on the theory outlined by Burke et al. [13]. In addition, many other re-
searchers developed their own codes, but these are generally unpublished and the developers were mostly interested in special
applications, particularly (quasi)-two electron problems such as electron scattering from alkali-like targets.
One of the principal ingredients of the above-mentioned programs is the usage of one set of orthogonal one-electron orbitals:
one subset is used to construct the target states in multi-configuration expansions and the other one is employed to represent
the scattered projectile. It simplifies the calculations and allows one to develop fast effective programs. On the other hand, it
276 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

leads to three major problems that often arise when the current suite of R-matrix codes is applied to truly complex targets. These
problems include (i) the difficulties in describing all target states of interest for a given calculation to sufficient accuracy, (ii) the
likely occurrence of unphysical structures, so-called “pseudo-resonances”, when an attempt is made to address the former problem,
and (iii) numerical difficulties due to an ill-conditioned orthogonalization procedure and the need to modify the so-called “Buttle
correction”.
Given the success of the R-matrix method for atomic collision processes, it is not surprising that many modifications of the
R-matrix codes and more detailed formulation of the algorithms have been suggested since the programs were first published. They
include but are not limited to the eigenchannel formulation [14], different choices of basis functions [15], particularly targeted to
eliminate or to reduce the sometimes problematic Buttle correction [16–18], or recipes to avoid the pseudoresonance problem that
occurs in many sophisticated applications [19,20].
The present program is another implementation of the R-matrix method with two significant innovations compared to the existing
codes:

• non-orthogonal orbitals can be used to represent both the bound and continuum one-electron orbitals;
• a set of B-splines defines the R-matrix basis functions.

The use of non-orthogonal bound orbitals generally allows us to achieve a much higher accuracy in the description of the target
states. Now they can be optimized in separate calculations, and high level of accuracy can be achieved with compact configuration
expansions. The use of non-orthogonal continuum orbitals can drastically reduce the pseudo-resonance problem. Indeed, in stan-
dard applications of the R-matrix method, certain (N + 1)-electron bound configurations must be included in the close-coupling
expansion in order to compensate for orthogonality constraints imposed on the continuum orbitals. However, it can be difficult to
keep the expansion fully consistent, and any inconsistency may lead to an undesirable pseudo-resonance structure. The use of a
non-orthogonal orbital set allows us to avoid the introduction of these (N + 1)-electron terms completely.
It follows from the above discussion that the non-orthogonal orbital technique may indeed overcome some of the major problems
associated with the standard application of the R-matrix method. The use of non-orthogonal orbitals was so far avoided because of
the fact that the most time-consuming part of atomic structure calculations is connected with the angular integration to construct
the Hamiltonian matrix elements. When a set of efficient general codes had been developed for this task for fully orthogonal
orbitals [21,22] and with some restricted non-orthogonality [23], it became possible to automate to a large extent this part of the
inner-region problem in conventional R-matrix calculations. This resulted in a wide application of this method. A key step for
the present program is the recent appearance of efficient codes for dealing with non-orthogonal orbitals in a general way [24,25].
Hence, it is now possible to use any non-orthogonal basis in the atomic structure calculations.
The choice of B-splines as basis functions has some advantages. They were introduced to atomic structure calculations in the
1980s and since became widely used due to their excellent numerical approximation properties (see, for example, the review by
Bachau et al. [26]). The B-splines are bell-shaped piecewise polynomial functions of order ks (degree ks − 1), defined by a given
set of points in some finite radial interval. Since the characteristic feature of the B-spline method is the solution of the Schrödinger
equation within a box, which is very similar to the inner region in the R-matrix method, it can be expected that B-splines will be
very effective in forming the R-matrix basis. Another important property of B-splines is that they form an effectively complete
basis on the interval spanned by the knot sequence. This completeness of the B-spline basis ensures that no Buttle correction to the
R-matrix elements is required. The use of B-splines as the R-matrix basis set was first outlined by van der Hart [27], who applied
the method to e–H scattering and obtained excellent agreement with existing benchmark results.
As before, the final step in the calculation requires the solution of the scattering problem in the external region and to match
the solutions on the boundary r = a. The present program deals only with the inner-region problem. Note, however, that once
the R-matrix has been constructed, no further changes are required. Any improvements achieved for solving the outer-region
problem, for example, through packages such as STGF (Seaton, unpublished) or FARM [28], can therefore be taken advantage of
immediately in the present approach as well. Also, reconstructing the wavefunction, as needed for calculations of photoionization
and photodetachment processes, proceeds along the same lines as in the standard implementation. The present programs can also
be used for bound-states calculations in the B-spline basis.
The atomic R-matrix programs described here have been successfully applied to several problems. First, Zatsarinny and Froese
Fischer [29] outlined a general B-spline R-matrix approach and demonstrated a successful implementation by treating photoioniza-
tion of lithium. Zatsarinny and Tayal [30,31] then applied the method to electron scattering from sulfur and oxygen, and Zatsarinny
et al. [32] achieved excellent agreement with experimental benchmark data for photodetachment of the He− (1s2s2p)4 Po metastable
state and the 1s-photodetachment of C− and B− [33]. These calculations were performed as a proof-of-principle to show a work-
ing prototype version. While the results are most impressive, they could be (and in some cases were) also generated with other
formulations of the R-matrix method, including the current Belfast implementation. During the past two years, the present code
was extensively tested for truly complex systems, including electron scattering from the noble gases [34,35], where considerable
improvement in comparison to results of standard R-matrix calculations was obtained.
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 277

This write-up is organized as follows. In the next section we will briefly summarize the most important aspects of the R-matrix
theory. The general features of the present BSR complex are given in Section 3, while each program in BSR is then described
in detail in Sections 4–11. Section 12 presents the libraries, common to all programs, including the B-spline library BSPLINE,
which is the key part of the present implementation. Additional program utilities, which can be used for preparing and processing
input/output data, are given in Section 13, and finally Section 14 describes the test runs.
Unless otherwise specified, all quantities are expressed in atomic units.

2. General theory

2.1. The close-coupling expansion

The problem of low-energy electron scattering from N -electronic atomic target is reduced to solving the time-independent
Schrödinger equation

(H − E)Ψα (Γ X, xN+1 ) = 0 (2.1)

with appropriate boundary conditions. The collision wave function Ψα (Γ X, xN+1 ) represents a fully antisymmetrized wave func-
tion of the system “target atom + projectile electron”, where X ≡ (x1 , x2 , . . . , xN ), xi ≡ (r i , σi ) with spatial (r i ) and spin (σi )
coordinates of the ith electron, Γ is a complete set of quantum numbers of the (N + 1)-electron system, and E is the total energy.
The subscript α characterizes the initial conditions and usually denotes the incoming scattering channel.
The Hamiltonian HN+1 which describes the scattering of an electron on an N -electron atomic target with the nuclear charge Z
has the form


N+1
1 Z


N+1
1
HN+1 = − ∇i2 − + , (2.2)
2 ri rij
i=1 i>j =1

where ri,j = |r i − r j |, r i and r j being the vector coordinates of electrons i and j . The origin of the coordinate frame is set at the
target nucleus, which is assumed to have an infinite mass. For the time being we have neglected all relativistic effects. We introduce
a set of target eigenstates, and possibly pseudostates Φi , and their corresponding eigenenergies Ei by the equation

Φi |HN |Φj  = Ei (Z, N )δij , (2.3)

where the integration is carried out over all the space and spin coordinates of the target electrons. Then the total energy is E =
Ei + ki2 /2, with Ei being the energy of the target in the state i while ki2 /2, represents the kinetic energy of the projectile electron.
The target states Φi are expanded in terms of single-configuration basis states ϕj by


Φi (x1 , . . . , xN ) = ϕj (x1 , . . . , xN )cij , (2.4)
j

where the coefficients cij are determined by diagonalizing the target Hamiltonian. The basis configurations ϕj are constructed
from a one-electron bound orbital basis, usually consisting of physical self-consistent field orbitals plus possibly additional pseudo-
orbitals. The latter are included to represent correlation effects. Note that we do not assume the one-electron basis to be orthogonal,
as usually imposed in scattering calculations. This allows us to optimize the bound orbitals in independent calculations for each
target state and to use term-dependent one-electron radial functions.
The solution of (2.1) has to satisfy the boundary conditions of an incoming wave in some scattering channel α and outgoing
waves in this and all other channels. In the close-coupling approximation, this solution is expanded in terms of a complete set of
N -electron target wave functions Φi . The corresponding expansion coefficients play effectively the role of wave functions for the
incident electron. In practice, one uses

n
  1

m
ΨαΓ (x1 , . . . , xN+1 ) = Â Φ̄iΓ x1 , . . . , xN ; r̂N+1 σN+1 Fiα (rN+1 ) + cj χjΓ (x1 , . . . , xN+1 ), (2.5)
rN+1
i=1 j =1

where  is the antisymmetrization operator with respect to the exchange of any pair of electrons while Fiα (r) is the radial compo-
nent of the scattered electron wave function when the target is in the ith state. Φ̄iΓ is called a channel function, which is obtained
by coupling the target state Φi with the spin-angle function of the scattered electron. Let us denote the quantum numbers of the
scattered electron in channel i as ki li mli msi . The Hamiltonian (2.2) is diagonal with respect to the total orbital momentum L, total
spin S, their projections onto the chosen axis ML and MS , and the parity π of the total (N + 1)-electron system. Therefore, in the
278 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

expansion (2.5) it is convenient to use the total momentum representation in which

Γ ≡ γ LSML MS π (2.6)
and the channel functions Φ̄iΓ are defined according to the following coupling:
 

   
Φ̄iΓ x1 , . . . , xN+1 ; r̂N+1 σN+1 = Li MLi li mli |LML Si MSi 12 mi SMS
MLi mli MSi mi
 
× Φi (x1 , . . . , xN )Yli mli r̂N+1 χ 1 mi (σN+1 ). (2.7)
2

Here Ylm is a spherical harmonic, χ(σ ) is a spin function, and we use the standard notation for the Clebsch–Gordan coefficients.
The function Fiα (xN+1 ) of the incident electron describes both open and closed channels. Very often the Fiα (xN+1 ) in (2.5) are
called channel orbitals. If (E − Ei ) is positive, the channel is said to be “open”; otherwise the channel is “closed”. Note, however,
that the function Fiα (xN+1 ) is not quadratically integrable for open channels. Closed-channel radial functions must satisfy the same
boundary conditions at r = 0 and the same orthogonality conditions as the open channel functions, but the closed-channel functions
are quadratically integrable.
The first term in expansion (2.5) should also include the integration over the continuous spectrum of the target, which corresponds
to the virtual (or real) target excitation into the ionization continuum. A direct inclusion of this term in (2.5) would tremendously
complicate the computational problem since the channel index becomes a continuous variable and the number of channels is not
denumerable. Very often, therefore, this term is omitted. However, the continuum part of the close-coupling expansion was found to
be very important at intermediate scattering energies. It can be simulated, to some extent, by the inclusion of bound pseudo-states.
The so-called correlation functions χi are quadratically integrable functions, usually constructed of the same one-electron or-
bitals as the target states Φi . The correlation functions ensure that the expansion (2.5) is complete in the bound orbital basis space,
even though the continuum radial functions are chosen to be orthogonal to the bound orbitals as discussed below. In the central field
approximation, the atomic orbitals are represented in the form
1
ϕj (x) = Ylj mj (r̂)χ(ms |σ ) Pnj lj (r). (2.8)
r
Then it is usually demanded that for lj = li the orthogonality condition
∞
Pnj lj (r)Fiα (r) dr = 0 (2.9)
0

is satisfied. This condition does not follow from general principles and is introduced only for simplification of numerical calcula-
tions. The introduction of the correlation functions χj (x1 , . . . , xN+1 ) means that in spite of implying conditions (2.9) the second
sum in the expansion (2.5) permits us to take into account the possibility of the virtual capture of electrons in the unfilled subshell.
The minimally necessary set of correlation functions χj in this case can be estimated as follows. The terms in the first sum in (2.5)
can be schematically represented as Ckj lj , where C denotes some target state configuration. Then, after imposing the orthogonality
conditions (Pnl |Fiα ) = 0, it is necessary to include the correlation function with the configuration Cnl, that could require recoupling
of angular momenta in the case of equivalent electrons.
In the present implementation, the conditions (2.9) are optional. The use of non-orthogonal continuum orbitals allows us to
avoid the introduction of the correlation functions χj (x1 , . . . , xN+1 ), or to use them directly for additional inclusion of short-range
correlation. In this case they also can be generated independently from the target states.
Coupled equations for the radial components of the functions Fi (r) representing scattered electrons and the coefficients cj can
be obtained by substituting expansion (2.5) into the Schrödinger equation and projecting onto the target functions Φi and the L2
functions χi . After separating out the spin and angular variables and eliminating the coefficients cj , we obtain the following set of
coupled integro-differential (“close-coupling”) equations
 2 

d li (li + 1) 2Z
2
− 2
+ + k i Fi (r) = 2
2
(Vij + Wij + Xij )Fj (r) (2.10)
dr r r
j

satisfied by the radial components Fi (r). In this equation, li is the orbital angular momentum of the scattered electron and Vij , Wij
and Xij are the partial-wave decompositions of the local direct, non-local exchange and non-local correlation potentials. These
potentials are too complicated to write down explicitly except for the simplest cases. Instead they are constructed by general
computer programs. Eqs. (2.10) can also contain terms which arise from the orthogonality constraints on the scattering radial
functions Fi (r).
A set of different computational methods for solving Eqs. (2.10) to yield the scattering matrices and amplitudes, which can
be compared with experiment, has been developed over the past decades. These methods now form the basis for a number of
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 279

computer program packages, some of which are widely used. As examples, we mention the linear algebraic equation method [36],
the non-iterative integral equation method [37], and the convergent close-coupling method [38]. All these methods can be referred as
straightforward methods where the solution of Eq. (2.10) is obtained in either configuration or momentum space for each collision
energy. A promising new approach for the direct solution of the close-coupling equations (2.10), based on the B-spline basis and
notable for its simplicity which is a key to a successful computational implementation, has been put forward by Froese Fischer and
Idrees [39]. The method determines the required solution within a finite boundary, with no assumed boundary conditions. This is
not a limitation, provided that the asymptotic region is reached, so that the solutions can be matched to a linear combination of
true asymptotic solutions. The core of the algorithm involves the evaluation of the Hamiltonian and overlap matrix elements in the
B-spline basis,

Hij = Ψi , H Ψj , Sij = Ψi , Ψj  (2.11)

and the extraction of the eigenvectors relative to the minimum modulus eigenvalues of the non-Hermitian, energy-dependent matrix

A(E) = H − ES (2.12)
at each prefixed energy E. Here the very effective inverse iteration procedure is used. The B-spline approach has been further
generalized to the single- and multichannel continuum case, mainly with application to photoionization of few-body system such
as H, He, H− and He− [40–42].

2.2. The R-matrix method

The R-matrix method is one of many methods for solving the close-coupling equations (2.10). The important difference from
the straightforward close-coupling formalism is a separate treatment of two regions: an inner region, in which all the electrons are
fully interacting with each other and possible external fields, and an outer region, in which the continuum electron only feels a local
potential. Here the coupled equations (with simple long-range potentials) are solved for each collision energy and matched, at the
boundary r = a, to the solution in the inner region. However, instead of solving a set of coupled integro-differential equations in the
internal region for each collision energy, the (N + 1)-electron wavefunction is expanded in terms of an energy-independent basis
set and treated similarly to the N electrons in atomic bound states. Consequently, general computer codes written for bound-state
atomic structure problems can be used, with only slight modifications, to generate the scattering algebra.
In the internal region, the (N + 1)-electron wavefunction at energy E is expanded in terms of an energy-independent basis set,
Ψk , as


ΨE = AEk Ψk . (2.13)
k

The basis states for given total angular momenta are constructed as

  1

ΨkΓ (x1 , . . . , xN+1 ) = Â Φ̄iΓ x1 , . . . , xN ; r̂N+1 σN+1 uj (rN+1 )aijΓ k + χiΓ (x1 , . . . , xN+1 )bik
Γ
. (2.14)
rN+1
ij i

The ui in Eq. (2.14) are radial continuum basis functions describing the motion of the scattering electron. They are non-zero on the
boundary of the internal region and thus provide the link between the solution in the internal and external regions. The quadratically
integrable functions χi have the same meaning as in Eq. (2.5) and are assumed to be fully confined to the internal region.
Consider now the solution of the Schrödinger equation in the internal region [0, a]. We note that the Hamiltonian HN+1 is not
Hermitian in this region due to the surface terms at r = a that arise from the kinetic energy operator. These surface terms can be
canceled by introducing the Bloch operator

N+1  
1 d b−1
LN+1 = δ(ri − a) − , (2.15)
2 dri ri
i=1

where b is an arbitrary constant (in the present implementation we use b = 0). Note that HN+1 + LN+1 is Hermitian for functions
satisfying arbitrary boundary conditions at r = a. We then rewrite the Schrödinger equation in the inner region as

(HN+1 + LN+1 − E)Ψ = LN+1 Ψ. (2.16)

This equation can be formally solved in terms of the R-matrix basis functions Ψk , which are obtained as eigenfunctions of the
equation
 Γ   
Ψi HN+1 + LN+1 ΨjΓ int = EiΓ ΨiΓ ΨjΓ int , (2.17)
280 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

where the integration over the radial variables is restricted to the internal region. The formal solution of Eq. (2.17) can then be
expanded as
 Γ
 1  Γ 
Ψ = Ψ Γ Ψk LN+1 Ψ Γ int . (2.18)
k
k
Ek − E
Γ

Projecting this equation onto the channel functions Φ̄iΓ and evaluating on the boundary of the internal region rN+1 = a yields

 
dFjΓ
FiΓ (a) = RijΓ (E) a − bFjΓ , (2.19)
dr rN+1 =a
j

where we have introduced the R-matrix with elements

1
wik wj k
Γ Γ
RijΓ (E) = , (2.20)
2a
k
EkΓ − E
the reduced radial wavefunctions
  
−1
rN+1 FiΓ (rN+1 ) = Φ̄iΓ Ψ Γ , (2.21)
and the surface amplitudes
  
a −1 wik
Γ
= Φ̄iΓ ΨkΓ r . (2.22)
N+1 =a

The primes on the brackets in Eqs. (2.21) and (2.22) mean that the integration is carried out over all the electronic space- and spin-
coordinates except for the radial coordinate rN+1 of the scattered electron. Eqs. (2.19) and (2.20) are the basic equations describing
the scattering of electrons by atoms or ions in the internal region. Together with the following relations for the coefficients AEk in
Eq. (2.13),

 
1 dFi 1
AEk = (Ek − E)−1 wik (a) a − bFi = (Ek − E)−1 w T R −1 F , (2.23)
2a dr r=a 2a
i
they allow us to establish the total wavefunction ΨE in the inner region for any value of the total energy E given the values of the
scattering orbitals on the boundary. The R-matrix given by Eq. (2.20) is obtained at all energies by diagonalizing HN+1 + LN+1
once for each set of conserved quantum numbers Γ to determine the basis functions Ψk and the corresponding eigenenergies Ek .
The logarithmic derivatives of the continuum radial wavefunctions Fi (r) on the boundary of the internal region are then given by
Eq. (2.19).
An important point in the R-matrix method is the choice of the radial continuum basis functions uj in Eq. (2.14). Although
members of any complete set of functions satisfying arbitrary boundary conditions can be used, a careful choice will speed up
the convergence of the expansion (2.13). In the standard R-matrix approach developed by the Belfast group [4] numerical basis
functions satisfying homogeneous boundary conditions at r = a were adopted. This approach gives accurate results, provided
that corrections proposed by Buttle [16] to allow for the omitted high-lying poles in the R-matrix expansion are included. The
shortcoming of this approach is that all basis functions have the same (usually, zero) logarithmic derivative at the boundary of
the internal region. This yields a discontinuity in the slope of the resulting continuum orbitals. In the standard approach, the
basis functions ui are forced to be orthogonal to the bound orbitals Pnl used for the construction of the target wavefunctions. To
compensate for the resulting restrictions on the total wavefunctions, the basis states Ψk must contain the correlation functions χi as
in Eq. (2.14). In the case of complex atoms, when extensive many-configuration expansions are used for accurate representations
of the target state wave functions, this may lead to a very large number of the correlation functions χi which must be included
in the close-coupling expansion to compensate for the orthogonality constraints. The present program has the option to imply a
fully non-orthogonal set of basis functions uj in Eq. (2.14). It allows us to use the close-coupling expansions (2.14) without any
correlation functions at all. Nevertheless, in practical calculations some restricted orthogonality constraints are usually imposed on
the basis functions uj , depending on the physical model under consideration.
The other key point of the present approach is the use of B-splines instead of the one-electron basis functions ui (r) in the
R-matrix representation (2.14) of the inner region. Full description of B-splines is given in Section 12. The B-splines possess prop-
erties as though they were especially created for the R-matrix theory. They form a complete basis on the finite interval [0, R], have
a universal nature, and are very convenient in numerical calculations because they allow us to avoid finite-difference formulae. Here
we have to distinguish between using B-splines as another basis for representation of the one-electron orbitals and using B-splines
for generating the complete pseudo-spectrum for some one-electron Hamiltonian, as is done in many atomic structure calcula-
tions [26]. The present program provides both options. In the first case, the coefficients aij k are found from the diagonalization
(2.17) of the full Hamiltonian in the B-spline representation. Such an approach is more suited for bound-state calculations. In the
second case, we first perform a preliminary diagonalization of the Hamiltonian blocks corresponding to one channel. This generates
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 281

a complete set of one-electron orbitals for each channel. We can transform the Hamiltonian matrix to the new representation based
on these one-electron orbitals, and now we are able to reduce the dimension of the full interaction matrix by dropping some of the
basis orbitals, depending on the problem under consideration.
The boundary conditions in the B-spline basis define only the first and the last basis functions, which are the only non-zero
terms, respectively, for r = 0 and r = a. The boundary conditions for the scattering function at the origin are satisfied in the form
F (0) = 0 by simply removing the first B-spline from the basis set. For the definition of the R-matrix at the boundary (2.20), the
amplitudes of the wavefunctions at r = a, wki , are required. These values are defined by the coefficients of the last spline, which is
equal to one at the boundary. The summation over the entire expansion (2.13) now gives the surface amplitudes in a straightforward
manner. Details for using B-splines in the construction and diagonalization of the interaction matrix (2.17) are given in Sections 7
and 8, respectively, along with a description of the corresponding programs.

2.3. The external region

The next step in the calculation is to solve the scattering problem in the external region and to match the solutions on the boundary
r = a in order to obtain the K-matrices, S-matrices, or phase shifts. Since the radius a is chosen such that electron exchange is
negligible in this region, we can expand the total wavefunction in the form

  −1 Γ
Ψ Γ (x1 , . . . , xN+1 ) = Φ̄iΓ x1 , . . . , xN ; r̂N+1 σN+1 rN+1 Fi (rN+1 ), rN+1 > a, (2.24)
j

where the Φ̄iΓ are the same set of channel functions as those retained in expansion (2.5) and the Fi (r) are the analytic continuations
for r > a of the reduced radial wavefunctions defined by Eqs. (2.21). The radial functions Fi (r) satisfy the set of coupled differential
equations
 

n

Λ aλ
d2 li (li + 1) 2(Z − N ) ij
− + + k 2
i F i
Γ
(r) = 2 F Γ (r), i = 1, n (r
a), (2.25)
dr 2 r2 r r λ+1 j
j =1 λ=1

where n is the number of channel functions retained in expansions (2.5) and (2.24) while li and ki2 are the channel angular momenta
and energies. The interaction between channels is defined by the long-range potential with coefficients

 Γ 
N
  
aijλ,Γ = Φ̄i x1 , . . . , xN ; r̂N+1 σN+1  rkλ Pλ (cos θkN+1 )Φ̄jΓ x1 , . . . , xN ; r̂N+1 σN+1 , (2.26)
k=1

where cos θkN+1 = r̂k . r̂N+1 , Pλ (x) is a Legendre polynomial. The integral in Eq. (2.26) is carried out over all electronic space and
spin coordinates except for the radial coordinate of the scattered electron. In practical calculations, the long-range potential coeffi-
cients are the by-product of generating the interacting matrix (2.17) and are defined by the coefficients of the relevant Slater integrals
R k , describing the direct interaction between channels in the internal region. Note that aij0 = N δij because of the orthogonality of
the target wavefunctions Φi . This relation is used to check the target representation in the present program.
Eqs. (2.25) can be integrated outward from r = a and fitted to an asymptotic expansion at large r as described in [28]. If all n
scattering channels are open, the asymptotic form of the radial wavefunctions Fi (r) may be written in the form

F (r) ∼ k −1/2 (F + GK), (2.27)

r→∞

where we have written the channel momenta, k, as a diagonal matrix. The diagonal matrices F and G correspond to regular and
irregular Coulomb (or Riccatti–Bessel) functions in each scattering channel. The asymptotic expression (2.27) defines the reaction
K-matrix, K, appropriate for standing-wave boundary conditions. The K-matrix is related to the scattering S-matrix by the equation

1 + iK
S= . (2.28)
1 − iK
It may be used directly to calculate cross sections and other scattering observables.
The long-range coefficients αij together with the target energies and the definition of the structure of the close-coupling equations
constitute the information needed to solve the scattering problem in the external region. This problem is well developed, and there is
a set of general computer codes for its solution. An efficient code is the outer region program FARM [28], which uses the R-matrix
propogator technique. In the present implementation, we use the program ASYPCK [43] for treating the external region. Note that,
in addition to the K-matrix, we also need the outer region solutions at the boundary r = a for the calculation of photoionization
cross sections.
282 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

2.4. Radiative processes

In R-matrix theory, the photoionization cross-sections can be defined through the dipole matrix between the initial state Φ0 and
the R-matrix basis states Ψk , provided that all radial orbitals of the initial state are well confined to the inner region. The total
photoionization cross section for a given photon energy ω is
  

4 2 2 ωC  − 
σ (ω) = π a0 α (Ψ D Φ0 )2 , (2.29)
3 2L0 + 1 j
j

where D is the electric dipole operator. It could be either in the length or the velocity form, with C = 1 in the length form, and
C = 4/ω2 in the velocity form, and the photon energy being in Ry. The index j runs over different possible solutions, and the
other quantities have their usual meaning. The solutions Ψj− in (2.29) correspond to asymptotic conditions with a plane wave in the
direction of the ejected electron momentum k and ingoing waves in all open channels. The corresponding radial functions F − are
related to the F (r) with the K-matrix asymptotic form (2.27) via

F − = −iF (1 − iK)−1 . (2.30)

Expanding Ψj− in terms of the R-matrix states as in Eq. (2.13) and using the expressions (2.23), we find that

1
(Ψk D Φ0 ) T −1 −
(Ψj− D Φ0 ) = w R F j (a) (2.31)
a Ek − E0 − ω k
k

where (Ψk D Φ0 ) are reduced matrix elements between the initial state and the R-matrix basis functions.
In addition to the total cross section, another important quantity is the anisotropy parameter β, which defines the angular distri-
bution of photoelectrons. For example, for linearly polarized incident radiation the angular distribution of photoelectrons is given
by
dσ σ
= 1 + βP2 (cos θ ) ,
dk̂ 4π
where θ is the angle of the ejected electron relative to the axis of polarization, while for unpolarized radiation it takes the form

dσ σ β
= 1 − P2 (cos θ ) ,
dk̂ 4π 2
where θ is the angle of the ejected electron relative to the incident radiation beam. β is determined from [2]:

15(−1)L0 +Lf

(2l + 1)(2l  + 1)(2L + 1)(2L + 1)
1/2
β= −
j |(Ψ j D Φ 0 )|2
jj 
   
l−l  −iσl +iσl  L l Lf L 1 L0 l l 2
×i e (Φ0 D Ψj− )(Ψj− D Φ0 ),
l  L 2 1 L 2 0 0 0
where Lo and Lf are the angular momenta of the initial atomic state and final target states, respectively, and the indexes j and j 
stand for the continuum channels, associated with the given target state. σl stands for Coulomb phase, see (2.34). Here and below
the standard notation for the nj symbols is used.
In order to use the expression (2.31), we need the values of the solutions Fi (a) at the R-matrix boundary. They can be obtained
by matching the general asymptotic solutions of Eq. (2.25) to the solutions in the internal region at r = a. This can be done with
relation (2.19), which in matrix form reads

F = aRF  − bRF (r  a). (2.32)

Note that in the outer region we have no independent physical solutions, where no is the number of open channels, which is defined
by all the target states accessible at a given excitation energy. To relate the n × n-dimensional R-matrix to the no × no K-matrix
defined in Eq. (2.27), we introduce n + no linearly independent solutions sij (r) and cij (r) of Eq. (2.25) satisfying the boundary
conditions
 
sij (r)   sin θi δij , i = 1, n, j = 1, no ,
cij (r) ∼ cos θ δ
i ij , i = 1, n, j = 1, no , (2.33)
 r→∞  exp(−φ )δ , , i = 1, n, j = no + 1, n,
cij (r) i i j −no
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 283

where θi and φi define the asymptotic phases in the open and closed channels, respectively:
 
1 z z
θi = ki r − li π − ln 2ki r + arg Γ li + 1 + i ,
2 ki ki (2.34)
z
φi = |ki |r − ln(2|ki |r).
|ki |
Now we can rewrite the asymptotic form of the scattering wavefunction in the general form
F = s + cK (r
a). (2.35)
Substituting this into Eq. (2.32) and solving for K, we get
K = B −1 A, (2.36)
where
   
b b
A = −s + aR s  − s , B = +c − aR c − c . (2.37)
a a
This completes the evaluation of the reactance matrix K and the value F (a), provided that the asymptotic solutions sij (r) and cij (r)
are known. As mentioned above, a number of computer packages are available for obtaining these solutions, such as the program
ASYPCK [43], which we use in the present package.
The initial state Φ0 in the present program can be obtained either in an independent MCHF calculation, or in the framework
of B-spline bound-states calculations discussed in the next section. In general, it requires the evaluation of dipole matrix elements
between states with non-orthogonal orbitals. Details of the calculation of dipole matrix elements in this case are presented in
Sections 9 and 10.

2.5. Bound-state calculations

Electron–atom collision theory is concerned with states of an (N + 1)-electron system for which N electrons are bound in an
atom or atomic ion and one electron can escape to infinity. Such states may be represented accurately using the close-coupling
expansion (2.5). This expansion is also suitable for representing states with one electron highly excited and the other N electrons
more tightly bound. When the close-coupling approximation is used to calculate bound states of atomic systems, it is referred to as
the frozen-core (FCS) approximation, and the Φi are usually labeled ‘core’ rather than ‘target’ functions.
The FCS method has several advantages. It can readily be extended to highly excited states. The multichannel form of Eq. (2.5)
allows us to include explicitly the interaction between different Rydberg series, as well as the interaction of the Rydberg series with
perturbers that can be represented in the second part of the expansion. The energies and wavefunctions can be computed efficiently
with an accuracy comparable to that obtainable using the best alternative methods, from which a large amount of radiative data
can be generated. The same expansions can be used for close-coupling collision calculations as for FCS bound-state calculations.
Often a comparison of the calculated bound-state energies with experimental energies provides a check on the accuracy of the
collision calculations. At the same time, one can take advantage of the extensive experience accumulated from close-coupling
collision calculations and the codes developed for such calculations, and one can apply them to the study of Rydberg series. The
first extensive FCS calculations for many-electron systems were made using the computer program IMPACT [36] that solved
the resulting integro-differential equations (see, for example, [44–46]). Further developments of the FCS method for the study of
Rydberg series were made in connection with the R-matrix method [47]. An example of such a calculation is given by Berrington
and Seaton [48].
The accuracy of both the above-mentioned methods is restricted through the use of finite-difference methods for computing
the radial functions. More accurate numerical results can be obtained by using a spline FCS method for Rydberg series, in which
the wavefunctions of the outer electrons are expanded directly in B-splines in some finite region r  a, with a sufficiently large
value of a. Such a method is implemented in the present package. The zero boundary conditions are imposed by deleting from the
expansion the first and the last splines, i.e. the only splines with a nonzero value at the boundaries. We also delete the next to last
spline, in order to guarantee a zero derivative at the border for all bound solutions.
The choice of B-splines as basis functions has some advantages. The completeness of the B-spline basis ensures that, in principle,
we can study the entire Rydberg series. The number of physical states we can obtain in one diagonalization is defined by the box
radius a, which can easily be varied in the B-spline representation. Exponential grid was found to be quite suitable for bound
states, and it allows us to use a rather large radius with a relatively small number of B-splines. For example, in order to obtain
Rydberg states up to n = 10 in neutral atoms, it is often sufficient to choose the box radius equal to 300 a.u. and the number
of splines as 45. If we aim to study Rydberg series up to n = 20, we should increase the border radius to 1200 a.u. With an
exponential grid, this increases the number of splines only to 51. Hence, the size of the interaction matrix, which is proportional
to the number of splines, does not increase considerably. Of course, these numbers somewhat depend on the nuclear charge Z. For
284 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

very high Z, it is advisable to add a few splines at very small radii to achieve an accurate representation of the orbitals near the
nucleus.
The wavefunctions in this method are obtained for all radii and for all Rydberg states under consideration. There is no need
to obtain an asymptotic solution and to match it to the inner-region solution as in the R-matrix method discussed above. This
considerably simplifies the calculations and the codes, but the use of a finite B-spline basis limits the upper n-value for the Rydberg
states. Consequently, this method can be efficient for the study of moderately excited states with values of n not too large, typically
in the range 10–30.
Our implementation of the spline method differs from a previous one [6] through the use of non-orthogonal orbitals, both
for the construction of target wavefunctions and for the representation of the outer electron. It provides us with a great deal of
flexibility in the choice of the ‘target’ wavefunctions, which can be optimized for each atomic state separately, and in the in-
troduction of different correlation corrections. In the present program, the core–core correlation may be taken into account by
using extensive multi-configuration target states. The core-valence correlation can be introduced, in principle, in two ways, either
by using a large set of excited target states in the close-coupling expansion or by introducing additional (N + 1)-electron states,
specially designed for this purpose. The convergence of the close-coupling expansion can be very slow, and hence the first ap-
proach is much more time consuming. Nevertheless, our experience shows that this method gives a more accurate description of
the core-polarization potential, and the (N + 1)-electron terms in Eq. (2.5) are better used only for the inclusion of the short-range
correlation.
Another difference compared to earlier implementations of B-splines concerns the relativistic corrections. Usually (see,
e.g., [49,50]) relativistic corrections are introduced via the Breit–Pauli CI approach after having calculated nonrelativistic wave-
functions using the B-spline expansions for the relevant states. We include the Breit–Pauli operators directly into the secular
equations (2.17). Hence, our B-spline expansions for Rydberg electrons include both the relativistic and correlation correc-
tions.

2.6. Relativistic corrections

With increasing nuclear charge Z, relativistic effects become important both in the target description and the scattering wavefunc-
tions, even for low-energy electron scattering. In the current code, relativistic corrections can be included through the Breit–Pauli
Hamiltonian.
The Breit–Pauli Hamiltonian [51] can be considered as a first-order correction to the non-relativistic atomic Hamiltonian. It
includes all relativistic terms up to order (αZ)2 , where α is the fine-structure constant. According to their different angular symmetry
and different influence on the atomic structure, we recognize three parts of the Breit–Pauli Hamiltonian, namely
HBP = HNR + HRS + HFS . (2.38)
Here HNR is the ordinary non-relativistic many-electron Hamiltonian

1
2
1
1
N N
HNR = − ∇i − Z + . (2.39)
2 ri rij
i=1 i=1 i<j

The relativistic shift operator HRS commutes with L and S and can be written as
HRS = HMC + HD1 + HD2 + HOO + HSSC , (2.40)
where HMC is the mass correction term

α 2
 2 + 2
N
HMC = − ∇i · ∇i , (2.41)
8
i=1
while HD1 and HD2 are the one-body and two-body Darwin terms, i.e. the relativistic correction to the potential energy,
 
α2 Z
2 1
N
HD1 = − ∇i , (2.42)
8 ri
i=1
 
α2
2 1
HD2 = ∇i . (2.43)
4 rij
i<j

Next, HSSC is the spin–spin contact term

8πα 2

N
HSSC = − (s i · s j )δ(r i · r j ) (2.44)
3
i<j
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 285

and, finally, HOO is the orbit–orbit term

 
α 2
1 (p i · pj ) r ij (r ij · pi ) · p j
N
HOO = − + . (2.45)
2 rij3
i<j
rij rij3

The fine-structure operator HFS describes interactions between the spin and the orbital angular momenta of the electrons. It does
not commute with L and S but only with the total electronic angular momentum J = L + S. The fine-structure operator consists
of three terms,
HFS = HSO + HSOO + HSS . (2.46)
Here HSO is the nuclear spin–orbit term

α2Z
1
N
HSO = (l i · s i ), (2.47)
2 ri3
i=1
HSOO is the spin-other-orbit term

α 2
r ij × p i
N
HSOO = − (s i + 2s j ), (2.48)
2
i =j
rij3

and HSS is the spin–spin term

N  
1 (s i · r ij )(s j · r ij )
HSS = α 2 3
(s i · s j ) − 3 2
. (2.49)
i<j
rij r ij

The present code has the option of directly including all the Breit–Pauli terms in the interaction matrix. For most applications,
however, it is sufficient to only include the one-electron terms and the two-electron spin-other-orbit interaction. In contrast to the
Belfast R-matrix code, the present code does not use recoupling to include the relativistic terms. Hence, storage and time problems
associated with a possibly very large number of non-relativistic Hamiltonian matrices to be calculated first can be avoided. An
exception is made for the mass correction term (2.41). Direct inclusion of this term in the interaction matrix often leads to non-
physical R-matrix poles with very low energies. For this reason, the program has the option to include this term as a first-order
perturbation, as described further in Section 8.
The non-fine-structure interactions commute with S 2 , S z , L2 and Lz and can therefore be considered in the |βLSML MS  rep-
resentation. The fine-structure interaction only commutes with S 2 , L2 , J 2 and J z , where J and J z are the total electronic angular
momentum and its projection on the quantization axis, respectively. Here we must use the |γ LSJ MJ  representation for the total
configuration state functions


Φ(γ LSJ MJ ) = (LML SMS |J MJ )Φ(γ LSML MS ). (2.50)
ML MS

When we include the fine-structure interactions, we must define the pair-coupling scheme in the close-coupling expansion (2.5).
The present program has three options for the definition of the pair-coupling scheme. The first option is LSJ-coupling. It is applied
only in the bound-state calculations as
Li + l = L, Si + 1
2 =S and L + S = J . (2.51)
Here Li and Si define the target term while l and 12 are the orbital angular momentum and the spin of the valence electron. This
scheme allows us to use simpler LS expansions for the target states and also reduces the number of channels (Rydberg series) in the
FCS expansion.
The most popular scheme for scattering calculations is the jK pair coupling
Ji + l = K and K + 1
2 = J, (2.52)
where Ji is the total electronic angular momentum of the target state. The target expansions should now be specified for each fine-
structure level. The corresponding Hamiltonian matrix (2.17) then becomes much larger. This results in more coupled channels in
the close-coupling expansion (2.5), therefore requiring a considerable increase in the computational effort.
The last coupling option is the jj pair coupling scheme
l+ 1
2 =j and Ji + j = J, (2.53)
which can be useful in case of a strong spin–orbit interaction.
286 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

3. General features of the BSR complex

3.1. Preparation of target states

In the present implementation, the preparation of the target states is the responsibility of the user. The input target description
includes both the target orbitals and configuration expansions. The present complex was designed, first of all, for a detailed and
accurate study of low-energy electron/photon scattering from neutral atoms, where we believe the accuracy of target states is
equally important as the scattering model. Especially it concerns the study of low-energy resonance structures. There is a set of
different atomic-structure programs which can be used for generating the atomic target states. The most popular ones are the
MCHF Atomic Structure Package by Froese Fischer et al. [52], the CIV3 program by Hibbert [53], and SUPERSTRUCTURE [54],
but others can be used as well. All these programs use different formats for the presentation of the one-electron orbitals and the
configuration expansions. As a basic representation for the different atomic states we choose the format of the MCHF package.
Here the configuration expansion for a specific atomic state is given in a so-called c-file, which is a list of configurations in the
spectroscopic notation along with the relevant expansion coefficients. The corresponding one-electron radial functions should be
provided in the B-spline representation (we will call such files bsw-files). These files can be generated from the w-files in the MCHF
package, or from the Slater-type representation used in CIV3, by using the interface programs described in Section 13. Hence, each
input target state in the present package should be represented by a pair: name.c for the configuration expansion and [Link] for
the relevant one-electron radial functions.
This preparation of the target states looks like a serious complication in comparison with the Belfast code, where one can simply
provide a single set of proper one-electron orbitals. On the other hand, the present method allows for a separate optimization of
different target states, or even to use non-orthogonal orbitals for one target state. It thereby enables the user to generate accurate
target states with a relatively small number of configurations and to directly include important effects such as term-dependence and
relaxation. The generation of the bsw-files for the one-electron radial functions frees us from the use of a specific atomic-structure
package. Furthermore, the user chooses the appropriate R-matrix radius at this stage and can perform additional checks of the target
states.

3.2. Main programs and data files

The program modules are run sequentially as outlined in Fig. 1. Besides a set of c- and bsw-files for the target description, the
user should prepare three input files ([Link], target, bsr_par), which contain the basic information about the given run and are
read by all programs. [Link] contains the B-spline grid and should be created at the time of the target preparation (a detailed
description of [Link] is given in Section 12.2). target contains the description of all target states and partial waves for the given
run (see program BSR_PREP and BSR_CONF for details). bsr_par defines the input parameters which control the execution of
the programs. All input parameters in bsr_par have the format name = value, and name should be placed at the beginning of the
line; otherwise the corresponding assignment will be ignored. This format is similar to the FORTRAN NAMELIST, but it does not
require specific NAMELISTs and is rather convenient in practical calculations. Furthermore, all input data form bsr_par can be
alternatively provided in the command line, in the same format name = value. In this case the data from the command line overwrite
the data from the input file. Data specific for a given partial wave are more convenient to be provided in the command line, while
data common for all partial waves should be provided in the input file.
Other disk files, outlined in Fig. 1 in bold, are produced by the programs at each stage and are normally deleted after use.
The final output files are the H and D files produced at the end of the BSR_HD and BSR_DMAT programs. They contain all
the information from the internal region that is required for the external region codes to run scattering or photoionization cal-
culations. The present package only deals with the internal region, and the user can then choose the program for the treatment
of the external region. In the case of photoionization calculations, two additional programs, MULT and BSR_DMAT, should be
run, in order to treat the dipole operator. The photoionization cross sections can then be produced by the program BSR_PHOT.
This program, however, is based on the differential solver of Crees [43]. While very general, it is not the fastest program avail-
able.
Another specific feature of the present file system is that the intermediate files are generated for each partial wave individually.
Such files have an extension which depends on the partial wave under consideration. For example, cfg.001, cfg.002, etc., are the
configuration files for the first, second and further partial waves. Thereafter, the files which are generated for the specific partial
wave will be denoted by nnn in the extension. Such a file system sometimes generates a large amount of intermediate files, but
allows us to run the calculations for different partial waves in parallel. It is especially convenient when the computer has several
processors.

3.3. Others features

The BSR package is written in FORTRAN 95 and uses double precision. It is designed to be transportable and has been imple-
mented on different computers and operating systems.
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 287

Fig. 1. Structure of the BSR package. Notation: capital names stand for the programs; the vertical dashed line indicates the sequential execution of the programs,
from top to bottom; bold names stand for the most important intermediate data files.

All unit numbers and file names are defined within the program. They are collected in the corresponding modules param, along
with all other parameters used in the specific program. Hence, if the user wants to change the default values for some parameter, it
only needs to be done in one place.
All principal arrays in the BSR package are allocatable, and they are initialized dynamically according to the input parameters.
Consequently, the user should not worry about the dimensions and their correspondence to the input parameters. On the other hand,
it limits the complexity of the physical model (e.g., number of scattering channels, number of orbitals, and so on) which we can
treat at the given computer by the size of the RAM memory. (Virtual memory will not help much here.)
The BSR package is organized as a set of main programs, shown in Fig. 1, utility programs, which help the user to prepare
and process the input/output data, and four libraries BSPLINE, ZCONF, ZCOM and LAPACK, which contain routines common
to all programs. The BSPLINE library, one of the most important parts of the package, contains routines for using B-splines
in the atomic structure calculations. The ZCONF library contains routines that deal with the configuration and close-coupling
expansions, and ZCOM contains routines for a variety of auxiliary computations that are common for different atomic-structure
calculations. The LAPACK library ([Link] was chosen for the linear-algebra calculations, such
as matrix diagonalization or solving a system of linear equations. For the compilation of the BSR program, a UNIX/LINUX user
can use the make-files provided along with the source files. Other users can use these files as a guide for compilation in a specific
environment.

4. Program BSR_PREP

4.1. Outline of the BSR_PREP calculations

This program provides initial preparations for BSR calculations. Since the input target states may be generated in independent
calculations, it is possible that the same spectroscopic notation was given to the different orbitals. For example, we can use dif-
ferent term-dependent radial functions for the same nl-orbital, or simply use the same spectroscopic notation for the correlated
functions in different states. In order to distinguish between such orbitals, we use an additional ‘set number’ or ‘set index’, e.g.,
nl1 , nl2 , . . . , nla , nlb and so on. It is difficult and very inconvenient for the user to keep track of all such cases by hand. The
288 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

BSR_PREP program therefore analyses the input target files and, if necessary, assigns an individual set number to each orbital.
The same concerns also the additional sets of (N + 1)-electron configurations included in the close-coupling expansion (2.14) for
a more complete description of short-range correlations. We will call these configurations ‘perturbers’, because they are often used
in bound-state calculations to describe the states that perturb the Rydberg series under consideration. These additional (N + 1)-
electron configurations should be provided in separate files, and they should be distinguished from the configurations which are
automatically generated for compensation of the orthogonality constraints imposed on the continuum orbitals (see the description
of the BSR_CONF program).
The list of input target states and perturbers are provided by the user in the file target. BSR_PREP sequentially reads all c- and
bsw-files and compares the one-electron radial functions. They are considered the same if their overlaps or first moments differ by
less than a given tolerance:
 
P1 |P2  − 1  ε,
 
P1 |r|P1  − P2 |r|P2   ε, (4.1)
 
P1 |1/r|P1  − P2 |1/r|P2   ε.

The radial functions are considered orthogonal if their overlap

 
P1 |P2   ε. (4.2)
For a given subset of orthogonal orbitals, BSR_PREP assigns the same set number. There are practically no restrictions on the
number of different non-orthogonal orbitals.
Then BSR_PREP rewrites the target c- and bsw-files with a new consistent spectroscopic notation. New target c-files receive
the names targ_001.c, targ_002.c, and so on, and they are ordered according to the target energies. All target radial functions are
placed in the single file [Link]. The same procedure is then applied to the perturber files, and each perturber is now rewritten
in the pair of files pert_nnn.c and pert_nnn.bsw, where nnn denotes the corresponding partial wave. All the above rewritings are
recorded in the file target, from which the user can obtain the relation between the new and the input names.
The only input parameter for BSR_PREP, which can be provided in the command line or in the input file bsr_par, is eps_ovr –
the value of the tolerance in Eqs. (4.1), (4.2). Its default value of 10−7 is usually sufficient for most applications.

4.2. Data files

target File type: formatted sequential input.

Written by user.
Read and modified by program BSR_PREP and BSR_CONF.
Description: list of target states and scattering channels.
Format: see Section 4.3.
bsr_par File type: formatted sequential input.
Written by user.
Description: input parameters for given run.
Format: see Section 3.2 and examples in Section 14.
[Link] File type: formatted sequential input.
Written by user.
Read and modified by routine define_grid from BSPLINE library.
Description: input parameters that define the B-spline grid.
Format: see Section 12.2.
name.c File type: formatted sequential input.
Provided by user.
Read by routine r_confc from ZCONF library.
Description: configuration expansion for one target state.
Format: see Section 12.5.
[Link] File type: unformatted sequential input.
Created by utilities w_bsw or slater_bsw (see Section 13).
Read by routine read_bsw.
Description: contains B-spline representation of one-electron radial functions.
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 289

targ_nnn.c File type: formatted sequential output.

Created by program BSR_PREP.
Read by program BSR_CONF.
Description: contains configuration expansion for one target state nnn, with consistent spectroscopic nota-
tion for one-electron orbitals.
[Link] File type: unformatted sequential input.
Created by program BSR_PREP.
Read by program BSR_MAT.
Description: contains all one-electron target orbitals in the B-spline basis.
pert_nnn.c File type: formatted sequential output.
Created by program BSR_PREP.
Read by program BSR_CONF.
Description: contains configuration expansion for perturber nnn, with consistent spectroscopic notation for
one-electron orbitals.
pert_nnn.bsw File type: unformatted sequential output.
Created by program BSR_PREP.
Read by programs BSR_MAT.
Description: contains one-electron orbitals in B-spline basis for perturber nnn.
bsr_prep.log File type: formatted sequential output.
Written by program BSR_PREP.
Read by user.
Description: running information.

4.3. Input data in file target

This file is created by the user and contains the description of the physical scattering model under consideration. An example of
an input target file is given in Fig. 2. In addition to the data, the file also contains dashed lines and comments placed behind the ‘!’
character. The comments are optional, whereas the dashed lines are used to delimit different types of data and are part of the format.
The first line is a title for the given run. Then, after a line-delimiter, the user provides some basic data such as the type of coupling
(LS or semi-relativistic JK, JJ), the number of electrons, and the atomic number. The format for all quantities is the same as for the
parameters in the file bsr_par: name = value, where name should be placed at the beginning of the line. The next block of data
is the list of target c-files, which should be placed after the number of target states, ntarg. Then the user provides the list of partial
waves, placed after the number of partial waves, nlsp. The description of each partial wave contains the total term (L, 2S + 1, π )
in LS coupling or (2J, 0, π) in semi-relativistic calculations. The perturber for a given partial wave, if any, is given by the name of
the corresponding c-file, otherwise the user should indicate no.
BSR_PREP adds some additional information to the target file. An example of target after a run of BSR_PREP is given in
Fig. 3. Now the target states are ordered by energy, and for each target state the program provides a new name targ_nnn.c, the term
(L, 2S + 1, π ), the energy, the number of configurations, and the number of new one-electron radial orbitals used in the description
of a given target state. All this information was extracted from the corresponding target c-files. Similar information is also provided
for the perturbers: a new standard name pert_nnn.c, the number of configurations and the number of new one-electron radial
orbitals used in the description of a given perturber. Also provided are the total number of target configurations, nct, the total
number of one-electron orbitals, nwt, and the list of orbitals with new spectroscopic notation, including the set numbers. In the
example shown, we used two different sets of orthogonal orbitals, for even-parity and odd-parity states, respectively. These sets
differ by their indexes 1 and 2. Note that the 1s orbital represents the closed shell 1s2 , and it is common for all configurations. These
core orbitals should be orthogonal to all orbitals, and hence they do not need any additional set index.

5. Program BSR_CONF

5.1. Outline of the BSR_CONF calculations

This program generates the close-coupling expansions (2.5) for the case under consideration, based on the information given
in the file target. First, the program reads the target configuration expansions from all targ_nnn.c files. For each partial wave,
the corresponding close-coupling expansions are then generated in LS-coupling, or in LSJ-, jK- or jj-coupling in the case of semi-
relativistic calculations, according to the input parameter coupling. The continuum configurations are simply generated from the
target configurations by adding the continuum orbital, in accordance with the angular momentum coupling rules. All possible
290 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

e + C
---------------------------------------------------------------
coupling = LS ! non-relativistic calculations
nz = 6 ! nuclear charge
nelc = 6 ! number of electron
---------------------------------------------------------------
ntarg = 36 ! number of target states
---------------------------------------------------------------
2p2_1S.c
2p2_1D.c
2p2_3P.c
2p3_3Do.c
2p3_5S.c
2p3s_3Po.c
2p3s_1Po.c
2p3p_1P.c
..........
3Po_ps4.c
3Do_ps3.c
---------------------------------------------------------------
nlsp = 64 ! number of partial waves
---------------------------------------------------------------
001 0 2 1 p_2Se.c
002 0 4 1 no
003 0 2 -1 no
004 0 4 -1 p_4So.c
005 1 2 1 p_2Pe.c
006 1 4 1 p_4Pe.c
007 1 2 -1 p_2Po.c
008 1 4 -1 no
009 2 2 1 p_2De.c
010 2 4 1 no
011 2 2 -1 p_2Do.c
012 2 4 -1 no
013 3 2 1 no
..........................
063 15 2 -1 no
064 15 4 -1 no
---------------------------------------------------------------
Fig. 2. Example of input file target prepared by the user.

continuum orbitals are considered, if only the additional restrictions are not indicated by input parameters max_ll, max_LT, or
max_ST. These parameters are mainly used in bound-state calculations, in order to get more compact expansions. Then, if there
is a pert_nnn.c file for a perturber, this expansion is also added without any change. Finally, if some orthogonal conditions are
imposed on the continuum orbitals, the program generates the additional (N + 1)-electron configurations, which are intended to
compensate for the orthogonality restraints imposed on the continuum orbitals. The resulting configuration expansion for each
partial wave is recorded in the file [Link], and the relevant channel information is added to the file target.
The continuum configurations in the close-coupling expansions in the [Link] files have the same expansion coefficients as the
corresponding target states, whereas the expansion coefficients for the (N + 1)-electron configurations are equal to unity. Further
modification of the expansion coefficients is introduced in the case of intermediate coupling. The configurations in the c-files are
written in consecutive LS-coupling of shells terms and have a definite total LS term. The intermediate coupling scheme is only
applied to the coupling of the target with the continuum electron. So, in order to guarantee that the given configuration represents
jK- or jj-coupling, we should apply the corresponding unitary transformation to the coupling between the target and the continuum
electron. The transformation coefficients are given by the following formulas:
 
TLS, jK = (l1 l2 )L, (s1 s2 )S J  (l1 s1 )j1 , l2 K, s2 , J
  
s2 +J −l2 −j1 1/2 L s1 K l2 l1 L
= (−1) [L, S, j1 , K] , (5.1)
s2 J S s1 K j1
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 291

e + C
------------------------------------------------------------------------
coupling = LS ! LS coupling
nz = 6 ! nuclear charge
nelc = 6 ! number of electrons
------------------------------------------------------------------------
ntarg = 36 ! number of target states
------------------------------------------------------------------------
2p2_3P.c targ_001.c 1 3 1 -37.77853957 258 10
2p2_1D.c targ_002.c 2 1 1 -37.72880614 263 0
2p2_1S.c targ_003.c 0 1 1 -37.67441594 132 0
2p3_5S.c targ_004.c 0 5 -1 -37.62901126 30 10
2p3s_3Po.c targ_005.c 1 3 -1 -37.50336422 318 0
2p3s_1Po.c targ_006.c 1 1 -1 -37.49458261 261 0
2p3_3Do.c targ_007.c 2 3 -1 -37.48150929 233 0
2p3p_1P.c targ_008.c 1 1 1 -37.46513100 169 0
........................................................................
3Po_ps4.c targ_035.c 1 3 -1 -36.72115047 318 0
3Do_ps3.c targ_036.c 2 3 -1 -36.69897493 233 0
------------------------------------------------------------------------
nct = 8091
nwt = 20
1s 2s1 2p1 3p1 3d1 3s1 9s1 9p1 9d1 9f1 2s2 2p2 3d2 3p2 3s2 4s2 9s2 9p2
9d2 9f2
------------------------------------------------------------------------
nlsp = 64 ! number of partial waves
------------------------------------------------------------------------
001 0 2 1 p_2Se.c pert_001.c 31 6
002 0 4 1 no no 0 0
003 0 2 -1 no no 0 0
004 0 4 -1 p_4So.c pert_004.c 28 9
005 1 2 1 p_2Pe.c pert_005.c 63 6
006 1 4 1 p_4Pe.c pert_006.c 40 6
007 1 2 -1 p_2Po.c pert_007.c 84 6
008 1 4 -1 no no 0 0
009 2 2 1 p_2De.c pert_009.c 55 6
010 2 4 1 no no 0 0
011 2 2 -1 p_2Do.c pert_011.c 141 10
012 2 4 -1 no no 0 0
013 3 2 1 no no 0 0
........................................................................
063 15 2 -1 no no 0 0
064 15 4 -1 no no 0 0
------------------------------------------------------------------------
Fig. 3. Example of the file target after a run of the BSR_PREP program.

 
l1 l2 L
TLS,jj = [L, S, j1 , j2 ]
1/2
s1 s2 S . (5.2)
j1 j2 J

Here the indices 1 and 2 denote the target and continuum orbital, respectively; LS is the total term of the given configuration;
J is the total angular momentum for the partial wave, and [L] stands for (2L + 1). The expansion coefficients for the continuum
configurations in the files [Link] are multiplied by the corresponding transformation coefficients, depending on the intermediate
coupling scheme under consideration. This procedure allows us to run all further calculations for generating and diagonalizing the
Hamiltonian matrix in the same way for all coupling schemes.
The most complicated point concerns the orthogonality constraints imposed on the continuum orbitals. The present package
allows the user to run calculations in different orthogonality modes, including the case when all orbitals are considered orthogonal
as in the standard R-matrix codes. The biggest advantage, however, can be gained by releasing some orthogonality constraints.
The orthogonality mode is defined by the input parameter jort. The default mode (jort = 1) is a partial orthogonality when the
continuum orbitals are orthogonal to a restricted number of bound orbitals and means the following:
292 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

(a) All continuum orbitals should be orthogonal to the closed-shell orbitals that are common to all states. Such shells will also be
referred to as core shells. Since these shells are closed, this does not lead to any additional (N + 1)-electron bound terms in the
close-coupling expansion in order to compensate these restrictions. The orthogonality of the scattering orbitals to closed-shell
core functions is physically justified. It avoids unphysical low-lying R-matrix poles corresponding to the capture of the incident
electron into these shells. Orthogonality to the closed core shell is automatically imposed in the present package, without any
additional input parameters.
(b) Some shells can be considered as ‘almost’ closed if their average occupation number in the given target expansion is close to
(4l + 2). For example, a configuration expansion for the ground state of neon, 1s2 2s2 2p6 , may contain only a few correlated
configurations with an excited 2s electron, and these configurations have relatively small expansion coefficients. In this case
we also can impose the orthogonality of the ks continuum electron to the 2s orbital without having to introduce additional
compensation configurations. Such cases are regulated by the input parameter c_orb. A shell is considered effectively closed,
when the average occupation number for the orbital is greater than c_orb × (4l + 2).

In the partial orthogonality mode, jort = 1, all additional orthogonality conditions should be indicated as kl|nl = 0. These
conditions can be placed in the input file bsr_par and they will be effective for all partial waves, or they can be placed at the end
of the [Link] file and then will only be effective for a given partial wave. If we introduce some orthogonal constraint, we should
generally also introduce the relevant (N + 1)-electron bound states in order to compensate for this constraint in the full functional
space. If we have a continuum configuration such as ‘target, kl’ and the continuum orbital kl is imposed to be orthogonal to the
bound orbital nl, then the compensation configurations have the structure ‘target, nl’, with proper recoupling for the equivalent
electrons. In doing this, however, we lost the information about the expansion coefficients for the individual target configurations.
As practical calculations show, not all target configurations here are of the same importance. The compensation configurations based
on the target configurations with small expansion coefficients may lead to an overestimate of the compensation effects. These, in
turn, may lead to overestimating the cross sections or to additional pseudo-structure. To avoid such situations, we introduce the
parameter c_orth, so that compensation configurations will only be generated for target states with expansion coefficients greater
than c_orth.
Another important situation is related to the case when different channels can generate the same configurations. Consider again,
for example, electron scattering from neon, with channels 2p6 kl and 2p5 nlkp. If the orbital kl has no orthogonal constrains, then this
channel can generate the 2p6 nl states. The same states will also be generated by second channel if the kp orbital is not orthogonal
to the 2p orbital. Such situation may lead to ‘overloaded’ overlap matrix in the generalized eigenvalue problem (2.17). As a result,
the overlap matrix will no longer be positive definite, and this prevents its diagonalization. To avoid such a situation, the user
should impose the orthogonality constraint kp|2p = 0 in the given example, without any compensation configurations, because
these configuration can still can be generated in the first channel. As a general advice, we suggest to impose all orthogonality
constraints which do not lead to compensation configurations. (The parameter icomp should be zero in this case.) It also simplifies
the generation of the Hamiltonian matrix and, in addition, avoids the numerical instability in the diagonalization of the generalized
eigenvalue problem (2.17).
Another important remark concerns the cases when orthogonality constraints are imposed to non-orthogonal but very similar
orbitals. For example, we can use term-dependent radial functions for the same nl orbital, or we can introduce relaxation effects in
the core-excitation. In the above example of e–Ne scattering, we can use two different radial functions for the 2p orbital in the 2p5
and 2p6 configurations to account for the relaxation effects. These functions are different but very similar. In this case we should
apply the kp|2p = 0 condition only to one of these orbitals.

5.2. Data files

bsr_par File type: formatted sequential input.

Written by user.
Description: input parameters for given run.
Format: see Section 3.2 and examples in Section 14.

target File type: formatted sequential input.

Written by user.
Read and modified by programs BSR_PREP and BSR_CONF.
Description: contains list of target states and scattering channels.
Format: see Section 4.3 and Section 5.4 below.
targ_nnn.c File type: formatted sequential input.
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 293

Created by program BSR_PREP.

Read by program BSR_CONF.
Description: contains configuration expansion for target states nnn.

pert_nnn.c File type: formatted sequential input.

Created by program BSR_PREP.
Read by program BSR_CONF.
Description: contains configuration expansion for perturber nnn.

cfg_nnn.c File type: formatted sequential output.

Created by program BSR_CONF.
Read by programs BSR_BREIT and BSR_MAT.
Description: contains configuration expansion for partial wave nnn.

bsr_conf.log File type: formatted sequential output.

Written by program BSR_CONF.
Read by user.
Description: running information.

5.3. Input parameters

Input parameters can be provided in the command line or in the input file bsr_par (data from the command line overwrite data
from the input file). Below we describe the data from bsr_par that are read by the program BSR_CONF. The default values for all
data are indicated in the brackets.

jort [1] orthogonality mode:

jort < 0 – full orthogonality;
jort = 0 – full non-orthogonality;
jort > 0 – partial orthogonality, default.
icomp [0] if = 1, the compensation configurations will be generated.
c_orth [0.75] tolerance for generation of compensation configurations, which will be generated only for target states with
expansion coefficients greater than c_orth.
c_orb [0.95] tolerance for closed shells. The continuum orbitals are supposed to be orthogonal to closed shells, and also
to almost closed shells, when the average occupation number for the orbital is greater than c_orb × (4l + 2).
max_ll [-1] restriction on the orbital angular momentum for the continuum orbitals
(the default value of −1 means no restrictions).
max_LT [-1] restriction on the total orbital angular momentum for the partial waves
(the default value of −1 means no restrictions).
max_ST [-1] restriction on the total spin for the partial waves
(the default value of −1 means no restrictions).
kl|nl = 0 additional orthogonality condition for the continuum orbital kl. If a set index indicated, this orthogonality
condition will be applied to specified orbitals; otherwise, it will be applied to all orbitals with the given l.
By default, all continuum orbitals are supposed to be non-orthogonal to the bound orbitals.

An additional input parameter coupling is read from file target. This parameter has three possible values LS, JK and JJ, and defines
the coupling for continuum states in the case of relativistic calculations. In this case, the input target states should be provided for
each J level, with the J -value being indicated in the first line of the corresponding c-file as 2J = value.

5.4. Modification of file target

BSR_CONF modifies the file target by adding information about the close-coupling expansions generated. An example of such
information is given in Fig. 4. First, for each partial wave, the program provides the number of channels, nch, and the number of
configurations, nc, in the corresponding [Link] file. The number of configurations is represented by two values, corresponding
to the continuum and the (N + 1)-electron parts of the given close-coupling expansion. Then, for each channel in a separate line,
we have the spectroscopic notation for the continuum orbital, its orbital momentum l, a pointer to the corresponding target state,
294 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

--------------------------------------------------
channels:
--------------------------------------------------
1. 001 nch = 22 nc = 5437 31
kd4 2 2 1 263 0
ks4 0 3 2 395 0
kp4 1 5 3 713 0
kp5 1 6 4 974 0
...................................................
kpf 1 35 22 5406 0
---------------------------------------------------
2. 002 nch = 11 nc = 3065 0
kp3 1 5 1 318 0
kd3 2 9 2 600 0
ks3 0 10 3 704 0
...................................................
kpa 1 35 11 3065 0
---------------------------------------------------
3. 003 nch = 13 nc = 2655 0
kp3 1 1 1 258 0
kd3 2 7 2 491 0
kp4 1 8 3 660 0
...................................................
kd9 2 36 13 2655 0
---------------------------------------------------
...................................................
---------------------------------------------------
64. 064 nch = 37 nc = 9689 0
ku1 15 1 1 258 0
kt1 14 5 2 576 0
kv1 16 5 3 894 0
kt2 14 7 4 1127 0
...................................................
kve 16 36 37 9689 0
---------------------------------------------------
max_ch = 66
max_nc = 15921
max_wf = 90
---------------------------------------------------
Fig. 4. Channel information added to target by the program BSR_CONF.

index of the channel, a pointer to the last configuration for the given channel in the file [Link], and the value of K, if we use the
jK intermediate coupling scheme. The spectroscopic notation for the continuum orbital is given only for the user’s information. As
additional information, the end of target provides the maximum number of channels, the maximum number of configurations, and
the maximum number of different one-electron orbitals used in the close-coupling expansion.

6. Program BSR_BREIT

6.1. Outline of the BSR_BREIT calculations

The BSR_BREIT program performs the angular integrations necessary to express the matrix elements of the Breit–Pauli Hamil-
tonian as a linear combination of radial integrals. Any amount of non-orthogonality between the orbitals may be present, leading
to overlap factors in the matrix elements. The program can effectively reuse data obtained previously by creating a databank for
angular coefficients. BSR_BREIT is a new optimized version of the general program BREIT_NO [24]. It can be used in many other
aspects of atomic-structure calculations. The calculations of angular integrals follow the method based upon the representation
of configuration wave functions through Slater determinants. The details of the present approach can be found in the long write-
up of the program BREIT_NO. Below we only provide a short description of the present version, outlining the most significant
modifications made in comparison to the original BREIT_NO.
The matrix element in the most general case of fully non-orthogonal orbitals has the form
     
 
Φ {nl}γ LSJ OΦ  {n l  }γ  L S  J  = ak R k (i, j ; i  , j  ) × D {nl}; {n l  } , (6.1)
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 295

where ak are the numerical coefficients, which depend only on the angular symmetry of the configuration state functions (CSFs)
involved, R k is the corresponding radial integral dependent on the operator O under consideration, and D({nl}, {n l  }) is the overlap
factor, which depends only on the sets of radial orbitals used in the construction of the wavefunctions Φ and Φ  , respectively.
Hereafter, the notation angular symmetry (AS) is used to denote the full set of angular momenta ({l}γ LS(L S  )), which includes
the momenta of spin–orbitals {l}, the terms of subshells αi Li Si , and the intermediate terms Li Si as well as the total momenta LS.
The symbol γ defines the coupling of momenta of all occupied subshells according to the given scheme
       
(core, L1 )L1 , L2 L2 , L3 L3 . . . L, (core, S1 )S1 , S2 S2 , S3 S3 . . . S, L + S = J . (6.2)

It is convenient to introduce the notation generalized shell term for (αL S) L S  , which also includes the corresponding intermediate
momenta L S  . It should be emphasized that the coefficients ak do not depend directly on the set of the principal quantum numbers
{n} and {n }. The overlap factor D({nl}, {n l  }) is generally a product of determinants corresponding to the matrices of one-electron
overlap integrals nl|n l; it depends only on the orthogonality conditions imposed on the one-electron radial functions. For example,
in the special case of orthogonal radial orbitals, the overlap factors are trivially reduced to 1 or 0, and we obtain the weighted sum of
radial integrals R k , which is ordinarily used in atomic structure calculations. For partial orthogonality conditions, some D-factors
disappear while others are simplified. Thus, once we obtain the matrix element in the case of fully non-orthogonal orbitals, we
can also obtain the matrix elements for all other configurations with the same AS and for all other orthogonal conditions, only by
analyzing the overlap factors D({nl}, {n l  }) and replacing the set {n, n } by the new one. Such a scheme is realized in the present
code. In large-scale calculations, it can lead to a considerable reduction of the execution time, because in practical calculations the
set of configurations contains many configurations with the same angular symmetry.
On the other hand, rather than record the angular coefficients ak for a given set of CSFs, we can save only the expansions for the
matrix between ASs with full non-orthogonality conditions. This data may be considered as a databank for angular coefficients ak .
They allow us to obtain the expansions of matrix elements for any set of CSFs constructed with the ASs included in the databank.
In case of new symmetries, we have to perform only the additional calculations and thus the databank can be easily extended.
The next idea concerns the structure of the angular coefficients ak themselves. Because the calculations of angular integrals
follow the method based upon the representation of configuration wavefunctions through Slater determinants, the coefficients ak can
be separated into a factor that depends only on the Slater determinants involved and a factor dependent on the coupling scheme (6.2).
Since many atomic states contain the same Slater determinants, it seems to be much more efficient to first calculate the matrix
elements between separate Slater determinants before multiplying them with the coupling-scheme factor depending on the AS
under consideration. The Slater determinants involved are defined only by the sets {l q } where q is the shell occupation number.
Hereafter, we will refer to the set {l q } as a configuration symmetry (CS). Hence, in order to avoid the recalculation of matrix
elements between the same Slater determinants, it is most effective to first define the list of all possible Slater determinants and then
to obtain the angular coefficients for the matrix elements between these Slater determinants.
The above ideas have been implemented in the present code. It requires additional sorting of the input configuration list according
to their ASs as well as to CSs. The principal features of the present version are given below.

6.2. Structure and data flow

The block diagram of the program BSR_BREIT, along with the data flow, is shown in Fig. 5. The only input data file is a c-
file with a list of the configurations under consideration. The default name for this file is [Link], or [Link], if a set of c-files for
different partial waves nnn should be handled. The range of required partial waves is given by the input parameters klsp1 and klsp2
(see Section 6.4). BRS_BREIT requires only a limited number of input parameters, most of which can be defined in the command
line in the format ‘name = value’. There is only one optional position parameter case.c which defines the name of specific input
c-file. Another input parameter is oper, which defines a set of operators from the Breit–Pauli Hamiltonian to be considered. All
input parameters are read from the command line by the routine read_arg.
Then the program executes fully independent calculations for each partial wave. First, the R_conf routine reads the input config-
urations in spectroscopic notation, decodes them to the internal ‘integer’ representation and tests the CSFs as to the correctness of
angular coupling, the number of electrons, parity and AFTER conditions. R_conf also allocates (or reallocates) all the main arrays
with dimensions specific for a given partial wave. Then the CSFs are sorted according to their angular and configuration symme-
tries. At this stage the program also checks if it is a continued calculation, i.e. if there exists a relevant databank from previous
calculations. The program analyses information contained in the databank and defines the subset of configurations and operators
that need additional consideration, or it decides that information in the databank is sufficient for the given input configuration list.
In the next block, pre_detexp, the determinant expansions for all configurations are generated and recorded in a scratch file. This
block is the main modification made in comparison to the previous version, BREIT_NO. It allows one to avoid a large amount of
recalculations for repeating determinant coefficients. The determinant expansion for a given angular symmetry is generated in the
routine det_expn, on the basis of angular coupling coefficients, calculated in detc_sh, and coefficients for determinant expansions
of individual shells, stored in the blockdata blk_det.
296 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

Fig. 5. Block diagram for the program BSR_BREIT and data flow (see text).

The calculations of a full set of angular coefficients and overlap factors for input angular symmetries are running in the block
conf_loop. First, routine det_breit computes the angular coefficients for matrix elements between all Slater determinants involved.
This information is accumulated in the module zoef_list. The cases of overlaps, one-electron and two-electron operators are treated
separately by the routines zno_0ee, zno_1ee, and zno_2ee, respectively. Routine zno_breit contains expressions for all Breit–Pauli
two-electron matrix elements in the nlms-representation. All these expressions are presented in the BREIT_NO write-up [24].
The corresponding angular coefficients for the integrals involved are accumulated in the module boef_list, in order to avoid their
recalculation. It is controlled by the routine check_boef.
Then, in routine term_loop, the angular coefficients between two given angular symmetries are calculated based upon the matrix
elements between the determinants involved, stored in module zoef_list, and the determinant expansion coefficients, which are read
from the scratch file nud (see block pre_detesp). These data are accumulated in the module coef_list. When all determinants for
the given angular symmetries are exhausted, the final angular coefficients are recorded by the routine add_res in the scratch file
nui, and module coef_list is initialized for accumulation of data from other matrix elements. The above calculations are repeated
in the block conf_loop for all combinations of input angular symmetries.
Finally, the new data and data from previous calculations are stored in file int_res.nnn. This file is then renamed to the initial
name int_bnk.nnn. Such a procedure prevents loosing the data in the initial databank if an error occurs when running BSR_BREIT.
The format of the output databank is given in Section 6.5. The output databank has the default name int_bnk for default input c-file
[Link]. For the specific input c-file case.c, the output databank has the name [Link].
All calculations in BSR_BREIT are carried out in full non-orthogonal mode, without specifications of principal quantum num-
bers for the orbitals involved. Instead, the positions of the corresponding shells in the configurations are used to specify the
one-electron orbitals involved. Consequently, the data obtained can be used further for all atomic states with the given angular
symmetries. In the present version we also abandoned the creation of an additional list of angular coefficients for the given set of
one-electron orbitals (the file [Link] in BREIT_NO). Instead, the matrix elements for specific orbitals are calculated dynamically
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 297

at the stage when the Hamiltonian matrix is generated in the program BSR_MAT. This approach has two benefits. First, we do
not loose the time for sorting the data in the [Link] file. Second, we save disk space, because the int_bnk files are usually much
smaller.

6.3. Data files

[Link] File type: formatted sequential input.

Created by program BSR_CONF.
Read by programs BSR_BREIT and BSR_MAT.
Description: contains configuration expansion for partial wave nnn.

int_bnk.nnn File type: unformatted sequential output.

Written by program BSR_CONF.
Read by BRS_MAT.
Description: databank for angular coefficients.

bsr_breit.log File type: formatted sequential output.

Written by program BSR_BREIT.
Read by user.
Description: running information.

6.4. Input parameters

Input parameters can be provided only in the command line in the format ‘name = value’. All data have the default values
indicated in the brackets.

klsp1 [0] first partial wave under consideration; default value klsp1 = 0 means that input configuration file is [Link]
and output data will be placed in file int_bnk.
klsp2 [klsp1] last partial wave under consideration.
oper [1110000] character*7, where each position can be 0 or 1, and indicate the operator under consideration:
oper(1) – overlaps
oper(2) – kinetic energy
oper(3) – two-electron electrostatic
oper(4) – spin–orbit
oper(5) – spin-other-orbit
oper(6) – spin–spin
oper(7) – orbit–orbit
Default value oper = 1110000 stands for non-relativistic calculations.
mk [9] maximum multipole index.
mrecl [100000] maximum record length (in bytes) in unformatted files. This parameter was introduced to insure the porta-
bility of the unformatted files on different computers and platforms.

6.5. Structure of the angular-coefficients databank

Summary of the output records in the int_bnk file.

1. nsymt,nsymc
nsymt – number of different angular symmetries
nsymc – number of different configuration symmetries
2. ASYM(1:nsymc)
List of configuration symmetries (character*26). Each configuration symmetry ASYM(.) contains (l(i),q(i),i=1,no),
LS in format (8(a1,i2),2a1), where no – number of shells; l(i) – electron orbital momentum in ith shell; q(i) –
number of electrons in ith shell; LS – total term in spectroscopic notation.
3. BSYM(1:nsymt)
298 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

List of angular symmetries (character*40). Each angular symmetry BSYM(.) contains a list of the generalized terms of shells
(A,L,S, LI,SI, 1:no) in format (8(i1,a1,2i1,a1)). Each generalized term includes: A – seniority number; L –
orbital momentum of the shell; S – spin; LI – intermediate orbital momentum; SI – intermediate spin.
4. IT_CONF(1:nsymt)
Indicates the corresponding configuration symmetry for a given angular symmetry, i.e. it defines the relation between the ASYM
and BSYM arrays.
5. nmatr
6. IT_OPER(1:7,1:nmatr)
nmatr = nsymt × (nsymt+1)/2 – dimension of the IT_OPER array.
IT_OPER(.,ij) indicates the operators that have already been considered by the program for angular symmetries i,j such
that ij = i × (i+1)/2+j, i > j.
IT_OPER(1:7,.) indicates the operator under consideration, in the same way as the input parameter oper (see Section 6.4).
7. ndet,kdet
8. KPD(1:ndet)
9. IPD(1:ndet)
10. NPD(1:kdet)
ndet – total number of overlap determinants.
kdet – size of the NPD array, containing the description of all overlap determinants.
KPD(i) – dimension of the ith overlap determinant, kd.
IPD(i) – pointer on the ith overlap determinant in array NPD, ip.
NPD(ip+1:ip+kd) – contains the description of the ith determinant as NPD(.) = i1 *bd + i2 , where bd = 215 –
packing basis for overlap determinants; i1 ,i2 – pointers to shells in the involved configurations.
11. ndef,kdef
12. KPF(1:ndef)
13. IPF(1:ndef)
14. NPF(1:kdef)
ndef – number of different overlap factors.
kdef – size of the NPF array, containing the description of all overlap factors.
KPF(i) – number of determinants in the ith overlap factor, kd.
IPF(i) – pointer on the ith overlap factor in array NPF, ip.
NPF(ip+1:ip+kd) – contains the description of the ith overlap factor as NPF(.) = ipd*bf + nd, where bf = 16 –
packing basis for overlap factors; ipd – pointer to the overlap determinant; nd – its power.
15. C,ij,int,idf repeat up to the end of file
C – angular coefficient for matrix element between angular symmetries i and j, recorded as ij = i × bc +j, where bc =215
– packing basis for configurations.
int – pointer on the relevant radial integral in the packing form int = m × b8 +k × b4 + i1 × b3 +i2 ×b2 +i3 ×b+i4 ,
where m – type of integral; k – multipole index; i1 ,i2 ,i3 ,i4 – pointer on the involved orbitals; b = 10 – packing basis.
idf – pointer on the corresponding overlap factor.

NOTE: all long arrays are divided into records with maximum length mrecl, in order to avoid a stack overflow on different platforms.

7. Program BSR_MAT

7.1. Outline of the BSR_MAT calculation

This program generates the interaction matrices (2.17) in the B-spline representation. It can be done either in LS-coupling or in
some intermediate coupling scheme with direct inclusion of selected Breit–Pauli terms. The eigenproblem (2.17) in the B-spline
basis leads to a generalized eigenvalue problem of the form
H c = ESc. (7.1)
The solution vector, c, can be written as
c = [a 1 , . . . , a Nc , b]T , (7.2)
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 299

where each a i is a column vector of B-spline coefficients for the given channel functions Φ̄iΓ
a i = [a1i , a2i , . . . , ans i ]T (7.3)
and b is the column vector of correlation functions coefficients
bj = [b1 , b2 , . . . , bNp ]T . (7.4)
Here Nc is the number of channels while Np is the number of correlation functions. Schematically, the interaction matrix H has
the form
 
H (11) H (12) . . . H (1Nc ) h(1p)
 H (21) H (22) . . . H (2Nc ) h(2p) 
 
 
 .. .. .. .. .. , (7.5)
 . . . . . 
 
 H (Nc 1) H (Nc 2) . . . H (Nc Nc ) h(Nc p) 
h(1p)T h(2p)T . . . h(Nc p)T h(pp)
where h(pp) is an Np × Np matrix that comes from the bound–bound interaction, the h(ip) are ns × Np matrices representing
the interaction between the ith channel and the bound states, and H (ij ) are ns × ns matrices for channel–channel interaction.
Furthermore, S in Eq. (7.1) is the overlap matrix. In the usual case of orthogonal conditions imposed on the scattering orbitals, it
reduces to a matrix consisting of diagonal blocks of the banded overlap matrix B between individual B-splines Bi
Bij = Bi |Bj , (7.6)
one for each channel. In the more general case of non-orthogonal orbitals, however, it has a more complicated structure (see Table 1).
Each individual matrix element in (7.5) is expressed by the program BSR_BREIT in the form

 
 
ci L(a, b) × D {nl}; {n l  } + ci R λ (a, b; c, d) × D {nl}; {n l  } , (7.7)
where the ci are numeric coefficients which depend only on the angular symmetry of the CSFs involved, L stands for one-electron
integrals, R λ are Slater integrals, and D({nl}, {n l  }) is the overlap factor, which depends only on the orthogonality of the radial
orbitals used in the construction of the CSFs. In general, the overlap factor is the multiplier of determinants of matrices consisting
of one-electron overlaps
 
 n1 l|n1 l n1 l|n2 l . . . n1 l|nq l 
 
 n l|n l n l|n l . . . n l|n l 
 2 1 2 2 2 q 
 . (7.8)
 ... ... ... . . . 

 n l|n l n l|n l . . . n l|n l 
q 1 q 2 q q
The scattering orbitals can appear in the radial integrals L and R λ , and also in the determinant factors. To derive the final expression
for the matrix elements in the B-spline basis, let us first simplify the overlap factors (7.8) by expanding them over the row (or
column) that contains one-electron overlaps with the scattering orbitals. At most, there can be one such row or/and column. The
residual overlap determinants depend only on the known bound orbitals and they can be calculated in a standard manner. It is
convenient to redefine the angular coefficients ci , multiplying them by these residual overlap factors. As a result, the initial overlap
factor D({nl}, {n l  }) in expression (7.7) is reduced to one- or two-electron overlaps of the form
kl|nl, kl|k  l or kl|nln l  |k  l  , (7.9)
where |kl stands for the radial function of the continuum orbital in the inner region


ukl (r) = ai Bi (r). (7.10)
i
Each term in (7.7) gives rise to a matrix in the spline basis. To elucidate the structure of these matrices, let us introduce, in addition
to the overlap B-matrix defined in Eq. (7.6), ns × ns matrices with elements
d2 l(l + 1) 2Z
L(..)ij = Bi | 2
− + |Bj ,
dr r2 r


(7.11)
R λ (.ab.)ij  = ai  bj  R λ (ij ; i  j  ), R λ (.a.b)ij = ai  bj  R λ (ij ; i  j  ),
ij jj 

where
  λ
R λ (ij, i  j  ) = Bi (r1 )Bj (r2 )U12 Bi  (r1 )Bj  (r2 ) . (7.12)
300 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

Here U12λ is the radial part of the λth term in the multipole expansion of 1/r , and a , b are the expansion coefficients in the spline
12 i j
basis for the bound orbitals Pa (r), Pb (r) as in (7.10). The four-dimensional array, R λ (ij ; i  j  ), effectively represents the Slater
integral in the B-spline basis. The integrals are zero if either |i − i  |
ks or |j − j  |
ks . Methods for computing these quantities,
which depend only on the basis, are given in Section 12. In order to describe the channel-bound interaction, let us introduce the
following vectors with elements


ns
B(.a)i = aj Bij , (7.13)
j =1


ns
L(.a)i = aj Lij , (7.14)
j =1


ns
R λ (.abc)i = aj bi  cj  R λ (ij ; i  j  ). (7.15)
j i  j  =1

With this notation, Table 1 lists the contributions of the different terms from the energy expression (7.7), which is usually produced
by the angular integration codes, into the interaction matrix (7.5), along with an indication of its structure. Table 1 also indicates
the contribution to the total overlap matrix arising from the matrix element on the right-hand side of Eq. (7.1). Note that the target

Table 1
Contribution of different terms (without indication of angular coefficients) in the interaction matrix. The symbols a, b, c, d stand for the bound orbitals, while ki
indicates the continuum orbital in channel i
Term Contribution Remarks
Channel–channel interaction
R λ (a, ki ; b, kj ) R λ (a.b.) Banded matrix,
R λ (ki , a; kj , b) R λ (.a.b) direct interaction
R λ (ki , a; b, kj ) R λ (.ab.) Full matrix,
R λ (a, ki , ; kj , b) R λ (a..b) exchange interaction
R λ (ki , a; b, c) d|kj  R λ (.abc).B(.d) Full matrix,
R λ (a, ki ; b, c) d|kj  R λ ([Link]).B(.d) contribution due to non-orthogonality
R λ (a, b; ki , c) d|kj  R λ (ab.c).B(.d)
R λ (a, b; c, ki ) d|kj  R λ (abc.).B(.d)
ki |ab|kj R λ (or L) B(.a).B(.b) ∗ R λ (or L) Full matrix,
contribution due to non-orthogonality
ki |kj R λ (or L) B(..) ∗ Ei Banded matrix, only for i = j
L(ki , kj ) L(..) Banded matrix, only for i = j
L(ki , a)b|kj  L(.a).B(.b) Full matrix,
contribution due to non-orthogonality
ki |kj  B(..) Banded matrix, only for i = j ,
contribution to overlap matrix
ki |ab|kj  B(.a).B(.b) Full matrix,
contribution to overlap matrix
Channel–bound interaction due to non-orthogonality
R λ (ki , a; b, c) R λ (.abc) Vector
R λ (a, ki ; b, c) R λ ([Link])
R λ (a, b; ki , c) R λ (ab.c)
R λ (a, b; c, ki ) R λ (abc.)
L(ki , a) L(.a) Vector
ki |aR λ (or L) B(.a) ∗ R λ (or L) Vector,
contribution due to non-orthogonality
ki |a B(.a) Vector,
contribution to overlap matrix
Bound–bound interaction
R λ (a, b; c, d) Scalar
L(a, b) Scalar
c Scalar, contribution to overlap matrix
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 301

wave functions are assumed to form an orthogonal set and diagonalize the target Hamiltonian. Hence, of all the matrix elements
containing kl|k  l, only the diagonal ones, with k = k  , contribute to the interaction matrix. The same concerns the matrix elements
with L(kl, k  l).

7.2. Structure and data flow

The block diagram of the program BSR_MAT, along with the data flow, is shown in Fig. 6. First, the program reads the data com-
mon to all partial waves. It includes the parameters for the calculation (routine read_arg), the B-spline parameters (define_spline),
and the description of the target states (read_target). Then the program executes fully independent calculations for each partial
wave. First, the program reads the configuration list from the relevant c-file [Link] and allocates all main arrays with dimensions
specific for a given partial wave. The most critical dimension here is the size of the interaction matrix, which usually occupies most
of the memory. The size of the interaction matrix is mainly defined by the product of the number of scattering channels and the
number of B-splines, because usually the (N + 1)-electron bound channels in the close-coupling expansion are either absent or
their number is relatively small in the present approach. In order to speed up the calculation, we place the entire interaction matrix
in the RAM memory. On the other side, this limits the number of scattering channels which we can include in the physical model
by the size of the computer memory. (Note that virtual memory will not help much here.)
The following calculations of the matrix elements are performed for each type of integrals sequentially. This part is controlled
by the routine state_res, whose flow diagram is placed at the bottom of Fig. 6. Each time this routine reads the relevant int_bnk
file from the beginning, chooses the integrals under consideration and provides their contribution to the total interaction matrix.

Fig. 6. Block diagram for the program BSR_MAT and data flow (see text).
302 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

Such a procedure reduces the number of integrals which should be considered at one time, and this increases the effectiveness
of their processing. One-electron and two-electron integrals are considered in different ways. For one-electron integrals, such
as kinetic-energy L-integrals (7.11) or spin–orbit Z-integrals (2.47), the program first calculates all possible integrals and their
B-spline representation for all one-electron orbitals in the configuration list, employing the routines gen_Lval and gen_Zval.
The relevant data are stored in the modules L_core and Z_core. Their storing requires only two-dimensional arrays and can be
handled even for a large number of one-electron orbitals. For the two-electron integrals, their storing requires four-dimensional
arrays. This cannot be realized in the case of a large number of different one-electron orbitals, which often occur in the present
approach. Hence, the two-electron integrals are calculated dynamically when they are required. In contrast to the standard R-matrix
approach, the present implementation uses a larger number of non-orthogonal orbitals for a more accurate representation of the
target states. Typical calculations may involve up to several hundred of different one-electron radial functions. The present code
has no practical limit on the maximum number of different one-electron orbitals that may be used to describe the target and the
scattering states.
The list of specific integrals needed for the given configurations is generated dynamically from the information in the databank
int_bnk. The integrals are then sorted according to their structure (see Table 1) and stored in the module cmdata in different ordered
blocks, using routines check_det, idef_type, add_coef, and add_cdata. The number and size of the individual blocks is defined by
the input parameters and could be adjusted to the speed and memory of a given computer. When all blocks are full, they are merged,
if necessary, with the routine merge_cdata and then are processed by the routine gen_matrix to get the B-spline representation
for specific integrals (using routines O_data, L_data, Z_data and I_data; the latter deals with two-electron integrals). All these
routines extensively use subroutines from the BSPLINE library. The corresponding contributions are then added to the interaction
matrix, located in module bsr_matrix.
The next step is the correction of the interaction matrix according to the orthogonal conditions imposed on the scattering orbitals,
using the routine BS_orth. This is done according to the procedure suggested by Bentley [55]. Suppose a continuum orbital, ukl (r),
in some ith channel is orthogonal to the bound orbital Pnl (r). This can be accounted for by modifying the Hamiltonian according
to

H → (1 − BccT )H (1 − ccT B), (7.16)

where c is the full solution vector (7.2) with all zero elements, except that the channel i is replaced by the B-spline expansion {pi }
of the orbital Pnl (r). This leads to transformation of only those blocks H (i, j ) in (7.5) which have one index equal to i. Recall that
the orthogonal conditions to the common core orbitals are implied automatically in the present implementation. All other conditions
should be specified through input data in the bsr_par or [Link] files.
At the last step, the interaction matrix, along with the long-range potential coefficients, is recorded in the relevant bsr_mat.nnn
file. The above calculations are repeated for each partial wave under consideration.

7.3. Data files

bsr_par File type: formatted sequential input.

Written by user.
Read by routine read_arg.
Description: input parameters for given run.
target File type: formatted sequential input.
Written by user and modified by BSR_PREP and BSR_CONF programs.
Read by routine read_target.
Description: contains description of target states and scattering channels.
[Link] File type: formatted sequential input.
Written by user.
Read by routine define_grid from BSPLINE library.
Description: input parameters that define the B-spline grid.
[Link] File type: formatted sequential input.
Created by program BSR_CONF.
Read by routine read_conf.
Description: list of configurations for given partial wave nnn.
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 303

bnk_int.nnn File type: unformatted sequential input.

Created by program BSR_BREIT.
Read by routine state_res.
Description: databank for angular coefficients.
[Link] File type: unformatted sequential input.
Created by program BSR_PREP.
Read by routine read_bsw.
Description: target’s one-electron orbitals in B-spline basis.
pert_nnn.bsw File type: unformatted sequential input.
Created by program BSR_PREP.
Read by routine read_bsw.
Description: perturber’s one-electron orbitals in B-spline basis, optional.
mat_log.nnn File type: formatted sequential output.
Written by program BSR_MAT.
Read by user.
Description: running information.
bsr_mat.nnn File type: unformatted sequential output.
Created by program BSR_MAT.
Read by program BSR_HD.
Description: Hamiltonian matrix and asymptotic coefficients.
int_mat.nnn File type: formatted sequential output.
Created by program BSR_MAT.
Read by user.
Description: list of one- and two-electron integrals. This output is only for debug purposes and controlled
by the parameter nud in the bsr_par input file.
7.4. Input data

Input data can be provided in the command line or in the input file bsr_par (the data from the command line overwrite the
data from input file). Below we describe those data from bsr_par which are read by program BSR_MAT. All data have the default
values indicated in the brackets. All default values, along with unit numbers and default file names, are placed in the module
param_mat.

klsp1 [1] first partial wave under consideration.

klsp2 [klsp1] last partial wave under consideration.
mk [7] maximum multipole index.
mb [5000] the size of one block in module cmdata used for the accumulation of a given type of integrals.
nb [500] number of blocks in module cmdata.
mrel [0] indicates the terms of the Breit–Pauli Hamiltonian which should be included
mrel=1 – include only one-electron scalar relativistic integrals.
mrel=2 – include additionally one-electron spin–orbit interaction.
mrel=3 – plus two-electron spin-other-orbit interaction.
mrel=4 – plus spin–spin interaction.
mrel=5 – plus orbit–orbit interaction.
mso [0] controls the inclusion of one-electron spin–orbit interaction. If set to 0, then a decision will be made ac-
cording to the mrel parameter; if set to −1(1), then the spin–orbit interaction will be excluded (included)
anyway.
msoo [0] controls the inclusion of the two-electron spin-other-orbit interaction in the same way as the mso parameter.
mss [0] controls the inclusion of the spin–spin interaction in the same way as the mso parameter.
moo [0] controls the inclusion of the orbit–orbit interaction in the same way as the mso parameter.
304 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

imvc [0] indicates the mode for processing the mass-velocity term:
imvc=+1 include mass-velocity term directly in the Hamiltonian (7.5).
imvc=-1 include the mass-velocity term only for the bound orbitals.
For scattering orbitals this interaction can then be included in the
BSR_HD program as a first-order correct
imvc=-2 the same as imvc=-1 but additionally exclude the mass-velocity term
for bound orbitals with principal quantum number nmvc.
This mode was introduced for consistency in cases when some
bound orbitals were generated in previous bound-state BSR calculations
without inclusion of relativistic corrections.
nmvc [0] the principal quantum number of orbitals for which the mass-velocity term is not included (see the imvc=
-2 mode above).
mrecl [100000] maximum record length (in bytes) in unformatted files. This parameter was introduced to insure the porta-
bility of the unformatted files on different computers and platforms.
nud [16] unit number for the debug int_mat.nnn file. If nud = 0 no debug file will be created.
eps_c [1.d-10] tolerance for angular coefficients combined with numerical overlap factor. Coefficients less than eps_c are
ignored.
eps_det [1.d-10] tolerance for the value of the numerical overlap factor. Overlap factors less than eps_det are taken as zero.
eps_soo [1.d-2] tolerance for the angular coefficients for the two-electron Breit–Pauli terms. The treatment of these terms
may take most of the time. On the other hand, they give only small corrections, and thus this parameter
allows us to consider contributions only from the principal configurations.
nphys [0] number of physical target states. If nphys = 0, all states are considered physical.
iphys [0] mode for processing the non-physical (pseudo) target states. If iphys > 0, all interactions between
pseudostates are ignored. If iphys = 1, the pseudostates interact only with the ground state.
maxnc [100000] size (in bytes) of buffer for coefficients read from int_bnk. The buffer was introduced in order to reduce the
number of read-write operations. It is especially important if several bsr_mat programs for different partial
waves are run with a single disk.

7.5. Output of Hamiltonian matrices and asymptotic coefficients

The main result of the bsr_mat program is the generation of the interaction matrix (2.17) in the B-spline representation (7.5).
The latter is stored in the bsr_mat.nnn files for each partial wave separately. Below is a summary of the output records in these
files.

1. ns,nch,kcp
ns – number of B-splines.
nch – number of channels.
ncp – number of bound (N + 1)-electron configurations.
The total dimension of interaction matrix is nhm=ns × nch+ncp.
2. om(i, 1 : i), repeat for each row i = 1,nhm – total overlap matrix.
3. hm(i, 1 : i), repeat for each row i = 1,nhm – total interaction matrix.
4. mk – maximum multipole index.
5. CF(1:nch,1:nch,0:mk) – matrix of asymptotic coefficients (2.26).
6. t(ns+1)
Array t(1:ns+ks) is the set of knot points used for the generation of B-splines.
Value t(ns+1) defines the maximum radius for the B-splines, i.e. the R-matrix radius.
7. th(1:ntarg,1:ntarg)
ntarg – number of target states.
th – interaction matrix between target states.
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 305

th is a by-product of the calculations of the total interaction matrix. This matrix should be diagonal with the target energies on
the diagonal. Comparison with the input target energies allows us to control the accuracy of the calculations.
8. to(1:ntarg,1:ntarg)
to – overlap matrix between target states. It should be equal to the unit matrix. A deviation from the unit matrix indicates that
the input target file is not correct.
9. EC – energy (in a.u.) of the common closed shells.

8. Program BSR_HD

8.1. Structure and data flow

This stage of the BSR complex deals with the final diagonalization of the Hamiltonian matrix. The eigenproblem (2.17) in the
B-spline basis leads to a generalized eigenvalue problem (7.1). For the solution of this problem, we use the standard routines from
the LAPACK library ([Link] If desired, the user can replace these routines with others.
The block diagram of the program BSR_HD, along with the data flow, is given in Fig. 7. First, the program reads the data common
to all partial waves. This includes the parameters of the calculation (routine read_arg), the B-spline parameters (define_spline),
and the description of the target states (read_target). Then the program executes fully independent calculations for each partial
wave under consideration.
At this stage, the program first reads the interaction matrix from the corresponding file bsr_mat.nnn, checking the option to use
experimental target energies instead of theoretical values. This option is controlled by the input parameter iexp. If this parameter
is greater than 0, the program reads experimental threshold energies from the file thresholds and calls the routine exp_thresh
for the corresponding modification of the interaction matrix. This modification is not as straightforward as in case of orthogonal
orbitals; in general, it needs an additional diagonalization procedure. For example, in order to shift all eigenvalues by the same
value Ec , we should add to the Hamiltonian matrix H the correction Ec S, where S is the total overlap matrix. Depending on the
orthogonality conditions imposed on the continuum orbitals, the total overlap matrix has a different structure. In case of orthogonal
orbitals, or for high partial waves, when the orbital angular momenta of the continuum orbitals are greater than the momenta of
the bound orbitals, the overlap matrix consists simply of diagonal blocks of the banded B-spline overlap matrix B (7.6). In this
case, the experimental-thresholds correction consists simply of adding (Eexp − Etheor )B to each diagonal block H (ii) in matrix
(7.5). This option is realized when iexp = 1. In case iexp = 2, the correction to the diagonal block H (ii) is (Eexp − Etheor )S(i, i),
where S(i, i) is the corresponding diagonal block in the total overlap matrix. This correction can be applied in the case of restricted
non-orthogonality, when the non-diagonal blocks in the total overlap matrix are small in comparison to the diagonal blocks. The

Fig. 7. Block diagram for the program BSR_HD and data flow (see text).
306 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

general procedure is applied if iexp = 3 and consists of the following steps. First, we introduce the orthogonal matrix O which
diagonalizes S according to
O T · S · O = s, (8.1)
where s is diagonal. We now define the square root matrix S 1/2 as
S 1/2 = O · s 1/2 · O T . (8.2)
Then the eigenproblem (7.1) can be written as
S −1/2 H S −1/2 (S 1/2 C) = (S 1/2 C)E, (8.3)
where E is the diagonal matrix of eigenvalues and C is the matrix of solutions. Now we can apply the usual procedure for adjusting
the target energies, i.e. add to each diagonal element of the modified Hamiltonian matrix the corresponding value of (Eexp − Etheor ).
If we denote these corrections by the diagonal matrix D, we have
(S −1/2 H S −1/2 + D)(S 1/2 C) = (S 1/2 C)E. (8.4)
Returning now to the original notation, we obtain
(H + S 1/2 DS 1/2 )C = SCE, (8.5)
i.e. the final correction is defined by the second term in the brackets. The above general procedure requires an additional diagonal-
ization of the overlap matrix and may be time-consuming in large-scale cases. Hence, it is recommended to apply this option only
for low partial waves where we expect a large contribution to the interaction matrix due to one-electron overlaps (see Table 1).
The next step is the diagonalization of the interaction matrix. It can be done in the full or simplified mode that is defined by the
input parameter idiag. In full diagonalization mode, idiag = 1, the only further modification of the interaction matrix is due to the
boundary conditions. Since the first B-spline is the only spline with nonzero value at r = 0, the zero boundary condition at r = 0
is imposed by removing the first B-spline from all expansions. In the code it means that we remove in the Hamiltonian matrix all
rows and columns related to the first spline. A further modification is related to the known dependence proportional to r l+1 for the
one-electron radial functions at small r. Recall that the first ks splines Biks behave as r i−1 at small r. Hence we can guarantee the
r l+1 behavior of the one-electron radial functions simply by removing the first (l + 1) splines from expansions. This can be done
only for orbitals with angular momentum l  ks − 2. For the remaining orbitals the r l+1 behavior is only approximately reproduced
by a suitable combination of the first splines. In bound-state calculations with the input parameter itype = −1, we additionally
remove the last two splines from the B-spline expansion, in order to guarantee that the resulting solutions have zero values and zero
first derivates at the border r = a.
When we diagonalize the Hamiltonian matrix in the B-spline basis, we get a very wide spectrum of eigenvalues, which spread
from threshold energies up to positive energies on the absolute energy scale. This is related to the effective completeness of the
B-spline basis, and hence no Buttle correction is required in this case. On the other hand, the B-spline expansions are comparatively
large, and in practical calculations it is not reasonable to treat all eigenvalues on equal footing. For the high-lying solutions, which
provide only a small contribution to the final R-matrix, we may use an approximate treatment. One of such treatments is realized
in the routine diag2_mat, which is called when the input parameter idiag = 2. In this case, we first diagonalize only the diagonal
blocks in the matrix (6.5) which are connected to the given scattering channels. Then the interaction matrix is transformed to this
new basis, and in the final expansion we retain only those solutions with energies below some value Emax. This parameter is an
input parameter and defined by the user on the basis of the physical model under consideration. Then the new simplified matrix is
diagonalized, and the resulting solutions are very close to the solution of the full matrix for energies below Emax. The remaining
solutions are approximated by those obtained at the diagonalization of the diagonal blocks. They can be considered as some kind of
Buttle correction. Such simplified treatment can reduce the size of the interaction matrix by up to factor of 2, and hence considerably
reduce the computational time.
The two-step diagonalization of the interaction matrix also allows us to introduce additional corrections at the level of first-order
perturbation theory. When we have zero-order solutions, ci , obtained from the diagonalization of the diagonal blocks, we can add
to the Hamiltonian matrix any perturbation as ci  |V |ci , where V is the perturbation potential. Such a possibility is realized in
the present implementation in the case of the mass-velocity correction for scattering orbitals, which is controlled by the parameter
jmvc. Note that the direct inclusion of the mass-velocity operator (2.41) in the Hamiltonian matrix leads to unphysical solutions
with very low energies.
After diagonalizing the Hamiltonian matrix, the results are treated in different ways, depending on the type of calculation
specified by the input parameter itype, namely, scattering calculations, itype = 0, photoionization calculations, itype = 1, or bound-
state calculations, itype = −1. For scattering calculations, including photoionization, the program first generates the corresponding
H file in the routine h_out. The output H file has a standard structure, first adopted in the Belfast R-matrix code. For a detailed
description of this file, see Section 8.4. Note that the surface amplitudes, which are needed to generate the final R-matrix, are simply
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 307

defined in the present implementation by the expansion coefficients of the last B-spline, which is the only spline with nonzero value
at the boundary. Another important difference from the earlier codes is that we do not provide any Buttle correction in the H file.
The user should remember that the resulting H file can only be used with the parameter ibuttle = 0 in further calculations, for
example, when running the FARM program.
In photoionization calculations, the program also records the full set of inner-region solutions in file [Link]. These data are
then used by program BSR_DMAT to generate the corresponding dipole matrix.
The last option in the BSR_HD program is a bound-state calculation. The resulting bound-state solutions are generated by the
routine b_out. In general, we can get a huge number of solutions, not all of which are physical or even helpful in further calculations.
A set of input parameters is therefore introduced in order to avoid the generation of a large number of output [Link] files with
excess information. The parameter msol limits the total number of output solutions, the parameters Emin and Emax restrict the
solutions according to their energies, the parameter it_max provides the highest target threshold for output solutions, and, finally,
the parameters n_min and n_max restrict the output solutions according to their principal quantum numbers. The principal quantum
numbers for bound-state solutions are defined according to their binding energies,

z2
Eb = . (8.6)
(n − µ)2
For a more accurate determination of the n values, the user may provide specific quantum defects µl through the parameters qd_0,
qd_1, qd_2 and qd_3 for orbital angular momenta l = 0, 1, 2, 3, respectively.
The routine b_out also tries to provide the spectroscopic assignment for the output solutions. To do this, the weights of different
channels and configurations are needed. Since the solutions vectors (7.2) have the metric

cT Sc = 1, (8.7)

the problem of weight determination is not trivial in the present case. For the determination of the weights we use the vector S 1/2 c,
where S 1/2 is defined in (8.2). It again requires the diagonalization of the total overlap matrix. The weights for different channels
and (N + 1)-electron configurations are determined in the subroutine w_out and saved in the file [Link] for further possible usage.

8.2. Data files

The following is a summary of input and output data files of BSR_HD. The unit numbers and file names are defined within the
program. Some files have an extension that depends upon the partial wave under consideration.

bsr_par File type: formatted sequential input.

Written by user.
Read by routine read_arg.
Description: input parameters for given run.
target File type: formatted sequential input.
Written by user and modified by BSR_PREP and BSR_CONF programs.
Read by routine read_target.
Description: contains the description of the target states and scattering channels.

[Link] File type: formatted sequential input.

Written by user.
Read and modified by routine define_grid from the BSPLINE library.
Description: input parameters that define the B-spline grid.
[Link] File type: formatted sequential input.
Created by program BSR_CONF.
Read by routine read_conf.
Description: list of configurations in spectroscopic notation for the partial wave nnn.

[Link] File type: unformatted sequential input.

Created by program BSR_PREP.
Read by routine read_bsw.
Description: target one-electron orbitals in the B-spline basis.
308 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

bsr_mat.nnn File type: unformatted sequential input.

Created by program BSR_MAT.
Read by routine read_bsrmat.
Description: Hamiltonian matrix and asymptotic coefficients.
thresholds File type: formatted sequential input.
Written by user.
Read by routine read_arg.
Description: experiment target energies, optional.
bsr_hd.nnn File type: formatted sequential output.
Written by program BSR_HD.
Read by user.
Description: running information.
[Link] File type: unformatted sequential output.
Created by routine h_out.
Read by program FARM, STG4, BSR_PHOT or other programs that deal with the outer region.
Description: diagonalized Hamiltonian matrix data for given partial wave.
[Link] File type: formatted sequential output.
Created by routine b_out.
Read by program BSR_DMAT and BSR_MCHF.
Description: bound-state solutions for given partial wave in the B-spline close-coupling representation (7.2).
[Link] File type: unformatted sequential output.
Created by routine r_out.
Read by program BSR_DMAT.
Description: full set of R-matrix solutions in the inner region.
[Link] File type: unformatted sequential output.
Created by routine w_out.
Read by routine b_out and program BSR_PHOT.
Description: weights of different channels in the inner-region solutions. Used for classification and sorting
of bound-state solutions.

8.3. Input data

Input data can be provided in the command line or in the input file bsr_par (the data from the command line overwrite the
data from the input file). Below we describe those data from bsr_par which are read by the program bsr_hd. All data have the
default values indicated in the brackets. All default values, along with unit numbers and default file names, are placed in the module
bsr_hd.

klsp1 [1] first partial wave under consideration.

klsp2 [klsp1] last partial wave under consideration.
itype [0] mode of calculations:
itype= 0 – scattering calculations; the [Link] file is generated.
itype= 1 – photoionization calculations; [Link] file is generated additionally.
itype=-1 – bound-state calculations; [Link] file is generated.
iexp [0] controls the treatment of experimental thresholds.
idiag [1] diagonalization mode:
idiag=1 – full diagonalization of interaction matrix; routine diag1_mat.
idiag=2 – simplified diagonalization; routine diag2_mat.
jmvc [0] pointer to the inclusion of the mass-velocity term as a first-order perturbation.
Parameters for bound-state calculations:
msol [0] restricts the number of bound-state solutions for output in [Link] file.
msol=0 means all solutions.
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 309

Emin [0] only solutions with E > Emin are considered (Emin = 0 means no restrictions)
Emax [0] only solutions with E < Emax are considered (Emax = 0 means no restrictions).
it_max [1] the highest target threshold for bound-state calculations.
n_min [1] the minimum principal quantum number of bound-state solutions.
n_max [99] the maximum principal quantum number of bound-state solutions.
qd_0 [3.0] quantum defect for s-orbitals.
qd_1 [2.5] quantum defect for p-orbitals.
qd_2 [2.0] quantum defect for d-orbitals.
qd_3 [0.5] quantum defect for f-orbitals.
The above values for quantum defects are used to classify the bound solutions.

8.4. Structure of the H file

We keep the same structure for the H files as in the Belfast R-matrix code [4]. The [Link] files from the BSR_HD program
contain information only for one partial wave. In order to generate the full [Link] file, the user should run the utility program
SUM_HH (see Section 13). Below is a summary of the output records in the [Link] files.

10. nelc,nz,lrange,km,ntarg,RA,BSTO
nelc – number of electrons in the target.
nz – atomic number.
lrange – maximum value of (l + 1) for the continuum orbitals.
km – maximum multipole index for the asymptotic coefficients.
RA – R-matrix radius.
BSTO {0.0} – the value of the logarithmic derivative to be imposed on the continuum orbitals.
The variables lrange and BSTO were included only for consistency with earlier versions.
11. etarg (1:ntarg)
12. ltarg (1:ntarg)
13. starg (1:ntarg),ptarg(1:ntarg)
etarg – target energies, in a.u.
ltarg – total L for target states (2J in case of intermediate coupling).
starg – total S for target states (0 in case of intermediate coupling).
ptarg – parity of the target states.
14. (0.d0,0.d0,0.d0,i=1,lrange)
This record in the standard [Link] file should contain the Buttle corrections, which are not used in the present implementation.
15. lpar,spar,ipar,nch,nhm,more
lpar – total L for given partial wave (2J in case of intermediate coupling).
spar – total S for given partial wave (0 in case of intermediate coupling).
ipar – total parity for given partial wave.
nch – number of scattering channels.
nhm – number of inner-region solutions.
more {0} – indicates no other information, except for the given partial wave.
16. NCONAT(1:ntarg) – indicates the number of channels associated with each target state.
17. lch(1:nch) – angular momenta of channel orbitals.
18. CF(1:nch,1:nch,1:mk) – matrix of asymptotic coefficients (2.26).
19. eval(1:nhm) – R-matrix poles in descending order (values Ek in Eq. (2.20)).
20. wmat(1:nch,1:nhm) – corresponding surface amplitudes (values wik in Eq. (2.20)).

8.5. Output of the R-matrix solutions

Summary of the output records in the [Link] files.

310 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

1. nhm,khm,kch,kcp,ns
2. eval(1:khm)
3. a(1:nhm,:)
repeat 3 khm times
nhm – dimension of the Hamiltonian matrix.
khm – number of solutions.
kch – number of channels.
kcp – number of (N + 1) configurations.
ns – number of B-splines for each channel.
eval(.) – R-matrix eigenvalues.
a(1:nhm,.) – corresponding eigenvectors.

8.6. Output of bound solutions

Summary of the output records in the [Link] files

1. ns,nch,ncp,nhm,nbound,L,S,P
2. i,LABEL
3. eval(i),Ebind(i),eff_n(i)
4. a(1:nhm,i)
repeat records 2–4 for each bound state.
ns – number of B-splines.
nch – number of channels.
ncp – number of (N + 1)-electron configurations (perturber).
nhm – dimension of solution vector in the B-spline representation.
LABEL – spectroscopic notation for solution i.
eval – energy in a.u.
Ebind – binding energy.
eff_n – effective principal number.
a – solution vector in the B-spline representation, see Eq. (7.2).

8.7. Output of weights

Summary of the output records in the [Link] files.

1. nch, ncp, khm

2. W(1:nch+ncp)
repeat record 2 for each solution
nch – number of channels.
ncp – number of (N + 1)-electron configurations (perturber).
khm – number of solutions.
W – array of weights for each channel and (N + 1)-electron configuration.

9. Program MULT

9.1. Outline of the MULT calculations

The MULT program performs the angular integrations necessary to express the matrix elements of the multipole tran-
sition operators as a linear combination of radial integrals. Any amount of non-orthogonality between the orbitals may be
present, leading to overlap factors in the matrix elements. The program can effectively reuse data obtained previously by cre-
ating a databank for angular coefficients. MULT is a new, completely revised and optimized version of the general program
ZAP_NO [24] and can be used in many other aspects of atomic-structure calculations. MULT can also be considered as a
straightforward extension of the program MCHF_MLTPOL [56] to the case of arbitrary orthogonal conditions between or-
bitals in the initial and final states. (In the following write-up, we try to keep the same notation as in MCHF_MLTPOL.)
The calculations of the angular integrals follow the method based upon the representation of configuration wavefunctions
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 311

through Slater determinants. The details of the present approach can be found in the long write-up of the program ZAP_NO.
Below we only provide a brief description of the present version, outlining the main modifications made in comparison
to ZAP_NO.

9.1.1. Definitions
The program deals with the transition operators Oµ[λ] for electric and magnetic transitions, which can be represented in tensor
forms, respectively, as


Eµ[λ] = e r λ (i)Cµ[λ] (i) (9.1)
i

and
!
 1 [1] 1
[λ]
Mµ[λ] = 2β λ(2λ − 1) r λ−1 (i)C [λ−1] (i) × l (i) + gs s [1] (i) . (9.2)
λ+1 2 µ
i

Here β = eh̄/2mc is the Bohr magneton, gs is the electron g-factor, and the summation is over the atomic electrons. In the Breit–
Pauli approximation, the initial-state and final-state wavefunctions are superpositions of LS-configuration state functions, i.e.


Ψ (J M) = ci Φ(γi Li Si J M). (9.3)
i

The dependence of the transition matrix elements on the magnetic quantum numbers can be found by applying the Wigner–Eckart
theorem,
 
[λ]   J −M J λ J
J M|Oµ |J M  = (−1) J O [λ] J  , (9.4)
−M µ M 
and the reduced matrix element (RME) J O [λ] J   itself can be written as


J O [λ] J   = ci ci  γi Li Si J O [λ] γi  Li  Si  J  . (9.5)
ii 

The operator Eµ[λ] in (9.1) associated with the electric transition Eλ is completely spin-independent, and we can use the “uncou-
pling” formula for its reduced matrix elements
 
 L S J
γ LSJ E [λ] γ  L S  J   = δSS  (−1)L+S+J +λ [J, J  ]1/2 γ LS E [λ] γ  L S  . (9.6)
J  λ L
The treatment of magnetic multipole transitions of any order is more complex. Usually, the magnetic operator (9.2) is split into
two parts with different spin and spatial dependence:
! 1 1
Mµ[λ] = 2β λ(2λ − 1) MA[λ] [λ]
µ + gs MBµ , (9.7)
λ+1 2
with

" #[λ]
MA[λ]
µ = r λ−1 (i) C [λ−1] (i) × l [1] (i) µ (9.8)
i

and

" #[λ]
MB[λ]
µ = r λ−1 (i) C [λ−1] (i) × s [1] (i) µ . (9.9)
i

The uncoupling formulae are different for these parts:

 
[λ]     L+S+J  +λ  1/2 L S J
γ LSJ MA γ LS J =δ SS  (−1) [J, J ] γ LS MA[λ] γ  L S   (9.10)
J λ L
and
 
L S J
γ LSJ MB[λ] γ  L S  J   = [J, J  , λ]1/2 L S J γ LS MB[λ] γ  L S  . (9.11)
λ−1 1 λ
312 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

The remaining problem is to evaluate the RME of the electric and magnetic multipole operator of any order between arbitrary
LS-coupled configurations. This is done in the present programs based upon the representation of configuration wavefunctions
through Slater determinants. This allows us to introduce a completely unrestricted version of a non-orthogonal scheme.

9.1.2. Calculation of RME between LS configuration wavefunctions

The present procedure is based upon the expansion of the N -electron configuration wave function in terms of Slater determinants

 
Φ(γ LSML MS ) = {u}|γ LS |u1 . . . uN , (9.12)
{u}

where ({u}|γ LS) denotes the expansion coefficient, and the individual determinant |u1 . . . uN  has the form
 
 t1 |u1  t1 |u2  . . . t1 |uN  
 
1  t2 |u1  t2 |u2  . . . t2 |uN  
|u1 . . . uN  = . (9.13)
N !  . . . ... ... . . . 
 t |u  t |u  . . . t |u  
N 1 N 2 N N
Here t = (r, s) represents the spin and position of the individual electrons and u stands for the one-electron quantum numbers,
|u = |nlmµ. Thus, the determination of matrix elements between CSFs is reduced to that of matrix elements between separate
Slater determinants, which, in turn, are reduced to one-electron integrals between spin–orbitals in the nlmµ-representation. In
comparison to the existing codes based upon the Racah technique, the resulting formulae explicitly include the magnetic quantum
numbers and contain additional summation over Slater functions and over all electrons but not shells. This increases considerably
the number of operations required, and the speed of the calculations is much reduced. On the other hand, the final formulae are
much simpler and more convenient for programming. The principal advantage, however, is that the nlmµ-representation procedure
facilitates the straightforward extension to the case of non-orthogonal orbitals to be made in the most general way. The main
problem here is to obtain the determinantal representation of arbitrary atomic wavefunctions, including configurations with several
open shells. An effective way to solve this problem has been proposed in the program ZAP_NO. It is based upon the preliminary
creation of tables of expansion coefficients for the individual shells. The determination of the coefficients for the CSFs is then
reduced to the vector coupling of shell terms. This procedure is used in the present code. For details we refer the reader to the
description of the program ZAP_NO.
When we get expansions (9.12) for the initial and final LS states, the total matrix elements for transition operators are re-
duced to matrix elements between Slater determinants. Consider a matrix element between two Slater determinants, |u1 . . . uN  and
|v1 . . . vN , for a transition operator O k = N λ λ
i=1 o (i), where o (i) is the one-electron operator with multipolarity λ, acting on elec-
tron i. Let us define the determinant Duv = det{ui |vj } of the matrix formed from all overlap integrals between the one-electron
orbitals ui and vj (i.j = 1, . . . , N) and let the Duv (ρ, ρ  ) denote the first cofactors of the determinant Duv , corresponding to the
matrix of overlap integrals with orbitals uρ and vρ deleted from the initial and final sets, respectively. Then, it can be shown [57]
that


u1 u2 . . . uN |Oµ[λ] |v1 v2 . . . vN  = (−1)ρ+ρ uρ |o[λ] |vρ  Duv (ρ, ρ  ). (9.14)
ρ,ρ 

This formula, in connection with irreducible tensor algebra, allows us to reduce the configuration matrix elements to one-electron
matrix elements between spin–orbitals nlmµ, multiplied by overlap factors in the case of non-orthogonal radial orbitals. Applying
the Wigner–Eckart theorem (9.4), the one-electron matrix elements in (9.14) can be reduced to one-electron RMEs. Here, however,
we should consider two cases, one related to the electric and MA operators, and one to the MB operator. In the first case, we have
no spin-dependence, and hence
 
[λ]     l λ l
nlmms |oµ |n l m ms  = (−1) l−m
nl O [λ] nl  δ(ms , ms ). (9.15)
−m µ m
For the MB operator (9.9), on the other hand, we should apply the formula for the matrix element of an irreducible tensor prod-
uct [58] of two operators, which only depend on the variables of the first and second subsystem, respectively:
   1 
[λ]     l λ − 1 l 1 1
nlmms |mbµ |n l m ms  = (−1) l−m
(−1) 1/2−ms 2 2
−m µl m −ms µs ms
× C(λ − 1, µl , 1, µs ; λ, µ)nl mb[λ] nl  , (9.16)
where C denotes a Clebsch–Gordon coefficient. Furthermore, µ, µl and µs are magnetic numbers defined by the total magnetic
numbers of the configurations under consideration: µl = ML − ML , µl = MS − MS , µ = µl + µs . In our case of electric and
magnetic transitions (9.1) and (9.2), the one-electron RME can be of three different types:
i e[λ] j  = i r λ C [λ] j  = (li C [λ] lj )Rijλ , (9.17)
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 313

i ma [λ] j  = i r λ−1 {C [λ−1] × l [1] }[λ] j 

$  
λ−1 1 λ
= (−1)li +λ+lj lj (lj + 1)(2lj + 1)(2λ + 1) (li C [λ] lj )Rijλ−1 , (9.18)
lj li lj
and
%
[λ] [λ−1] [1] 3
i mb j  = i r λ−1
C s j = (li C [λ−1] lj )Rijλ−1 , (9.19)
2
where Rijλ is the radial transition integral
∞
Rijλ = d r r λ Pni li (r)Pnj lj (r). (9.20)
0
The present algorithm for the calculation of matrix elements consists of two independent steps: the determination of the deter-
minantal representation for two input configuration wavefunctions with certain magnetic quantum numbers ML and MS , and the
calculation of matrix elements between all Slater states obtained by the above formulae. Our final task is to obtain the reduced
matrix elements for configuration wavefunctions, γ LS O [λ] γ  L S  . They are obtained on the basis of the calculated matrix ele-
ments γ LSML MS Eµ[λ] γ  L S  J  ML MS  by a reciprocal application of the Wigner–Eckart theorem. The corresponding formulae
can be obtained from Eqs. (9.15) and (9.16) by a straightforward replacing of the one-electron quantum numbers by the total ones.
Note that the amount of calculations in the present approach crucially depends on the number of Slater states involved, which, in
turn, depends on the chosen quantum numbers ML and MS . In the present implementation, we use ML = L and MS = S, which
guarantees a minimum number of relevant Slater states.

9.1.3. Output results

The MULT program provides, as a final result, the expressions for the reduced matrix elements of the electric and magnetic
operators (9.1) and (9.2) in the form
  & &  
 
Φ {nl}γ LS &O &Φ  {n l  }γ  L S  = ak d k (i, i  ) × D {nl}, {n l  } , (9.21)
where ak are the numeric coefficients which depend only on the angular symmetry of the configuration state functions involved,
d k is the corresponding radial integral dependent on the operator O under consideration, and D({nl}, {n l  }) is the overlap factor,
which depends only on the sets of radial orbitals used in constructing the wavefunctions Φ and Φ  , respectively. It should be
emphasized that the coefficients ak do not depend directly on the set of principal quantum numbers {n} and {n }. They only depend
on the orthogonality conditions imposed on the one-electron radial functions. For example, in the limiting case of orthogonal radial
orbitals, the overlap factors are reduced trivially to factors 1 or 0, and we obtain the weighted sum of radial integrals d k , which
is ordinarily used in atomic structure calculations. In the case of intermediate orthogonality conditions, some D-factors disappear
while others are simplified. Thus, once we have obtained the matrix element in the case of fully non-orthogonal orbitals, we can
also obtain the matrix elements for all other configurations with the same angular symmetry and for all other orthogonal conditions,
only by analyzing the overlap factors D({nl}, {n l  }) and replacing the set {n, n } by the specific one. Such a scheme is realized
in the present code. In large-scale calculations, this can lead to a considerable reduction of execution time, because in practical
calculations the set of configurations contains many configurations with the same angular symmetry.
On the other hand, rather than record the angular coefficients for the given set of CSFs, we can save only the expansions for the
matrix elements with different angular symmetries, supposing that all radial functions in the initial and final sets are non-orthogonal.
These data may be considered as a databank for angular coefficients. They allow us to obtain the matrix element expansions for
any set of CSFs constructed with the angular symmetries included in the bank. For new symmetries, we only have to perform the
additional calculations and the databank can be extended.

9.2. Structure and data flow

The block diagram of the program MULT, along with the data flow, is given in Fig. 8. The only input data files are c-files
with a list of configuration state functions for the initial and final states involved. The MULT program considers only one type of
transition operator in a single run. The type of the operator is given as E1, E2, . . . , M1, M2, . . . , where the first letter indicates if it
is an electric or magnetic transition while the second letter defines the multipole index. The parameter E0 indicates the calculations
of overlap matrix elements. Since MULT only requires a small number of input parameters, all of them must be defined in the
command line as arguments.
The r_conf routine reads the input configuration lists in spectroscopic notation, decodes them to the internal ‘integer’ represen-
tation, and tests the CSFs as to the correctness of angular coupling, number of electrons, parity and AFTER conditions. r_conf also
allocates (or reallocates) all main arrays with dimensions specific for the given partial wave. Then the CSFs are sorted according
314 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

Fig. 8. Block diagram for the program MULT and data flow (see text).

to their angular and configuration symmetries. At this stage the program also checks if this is a continuing calculation, i.e. whether
there exists a relevant databank from some previous calculations. The program analyses the information contained in the databank
and defines the subset of configurations that need additional consideration, or it decides that the information in the databank is
sufficient for the given input configuration lists.
In the next block, pre_detexp, the determinant expansions for all configurations are generated and recorded in the scratch
file. This block is the main modification made in comparison to the previous version, ZAP_NO [24]. It avoids the large amount
of recalculations of repeating determinant coefficients. The determinant expansion for a given angular symmetry is generated in
the routine det_expn, on the basis of angular coupling coefficients, calculated in detc_sh, and the coefficients for determinant
expansions of individual shells, stored in the blockdata blk_det.
The calculations of the full set of integral coefficients and overlap factors for the input angular symmetries are run in the
block conf_loop. First, the routine det_mult computes the radial expansions for matrix elements between all Slater determinants
involved. This information is accumulated in the module zoef_list. The cases of overlaps, zero-orbital or one-orbital difference in
the Slater determinants under consideration are treated separately by the routines zno_overlap, zno_000, and zno_001, respectively.
The routine rad_matr contains expressions for one-electron matrix elements in the nlmµ-representation (see Section 9.1.2). The
function idet_fact defines the corresponding overlaps factors Duv (ρ, ρ  ) in the expression (9.13).
Then, in the routine term_loop, the angular coefficients between two given angular symmetries are calculated on the basis of
angular coefficients between the individual determinants involved (stored in module zoef_list) and determinant expansion coeffi-
cients, which are read from the scratch file nud (see block pre_detesp). These data are accumulated in the module coef_list. When
all determinants for the given angular symmetries are exhausted, the final angular coefficients are recorded by the routine add_res
in the scratch file nui, and module coef_list is initialized to accumulate data from other matrix elements. The above calculations
are repeated in the block conf_loop for all combinations of input angular symmetries.
Finally, the new data and data from previous calculations are stored in the file mult_res. This file then is renamed to the initial
name mult_bnk. This procedure prevents loosing the data in the initial databank if an error occurs when running the MULT
program. The format of the output databank is given in Section 9.5.
All calculations in MULT are carried out in full non-orthogonal mode, without specifying the principal quantum numbers for the
orbitals involved. Instead, the positions of the corresponding shells in the configurations are used to specify the one-electron orbitals
involved. Hence, the data obtained can be used further for all atomic states with the given angular symmetries. In the present version
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 315

we also abandoned the creation of an additional list of angular coefficients for a specific set of one-electron orbitals. Instead, the
matrix elements for specific orbitals are calculated dynamically when the required transition data are generated.

9.3. Data files

name1.c File type: formatted sequential input.

Created by user.
Description: contains configuration expansion for the initial state.
name2.c File type: formatted sequential input.
Created by user.
Description: contains configuration expansion for the final state.
mult_bnk File type: unformatted sequential output.
Written by program MULT.
Read by BRS_DMAT.
Description: databank for angular coefficients of the dipole operator.
[Link] File type: formatted sequential output.
Written by program MULT.
Read by user.
Description: running information.

9.4. Input parameters

Input parameters must be provided in the command line.

name1.c name for c-file of initial state.

name2.c name for c-file of final state.
AA character*2, type of operator and multipole index: E1, E2, . . . , M1, M2, . . . .

9.5. Structure of the dipole databank

Summary of the output records in the mult_bnk file.

1. ktype,kpol
ktype – ‘E’ or ‘M’, for electric and magnetic transitions, respectively.
kpol – multipole index.
16. nsymt,nsymc
nsymt – number of different angular symmetries.
nsymc – number of different configuration symmetries.
17. ASYM(1:nsymc)
List of configuration symmetries (character*26). Each configuration symmetry ASYM(.) contains (l(i),q(i),i=1,no),
LS in format (8(a1,i2),2a1), where no – number of shells; l(i) – electron orbital momentum in ith shell; q(i) –
number of electrons in ith shell; LS – total term in spectroscopic notation.
18. BSYM(1:nsymt)
List of angular symmetries (character*40). Each angular symmetry BSYM(.) contains a list of the generalized terms of shells
(A,L,S,LI,SI,1:no) in format (8(i1,a1,2i1,a1)). Each generalized term includes: A – seniority number; L –
orbital momentum of the shell; S – spin; LI – intermediate orbital momentum; SI – intermediate spin.
19. IT_CONF(1:nsymt)
Indicates the corresponding configuration symmetry for a given angular symmetry, i.e. it defines the relation between the ASYM
and BSYM arrays.
316 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

20. IT_OPER(1:nsymt × (nsymt+1)/2)

IT_OPER(ij) indicates the operators, which have been already considered by the program for angular symmetries i,j,
such that ij = i × (i+1)/2+j, i > j.
21. ndet,kdet
22. KPD(1:ndet)
23. IPD(1:ndet)
24. NPD(1:kdet)
ndet – total number of overlap determinants.
kdet – size of the NPD array, containing the description of all overlap determinants.
KPD(i) – dimension of the ith overlap determinant, kd.
IPD(i) – pointer on the ith overlap determinant in array NPD, ip.
NPD(ip+1:ip+kd) – contains description of the ith determinant as NPD(.) = i1 *bd + i2 , where bd = 215 – pack-
ing basis for overlap determinants; i1 ,i2 – pointers to shells in the involved configurations.
25. ndef,kdef
26. KPF(1:ndef)
27. IPF(1:ndef)
28. NPF(1:kdef)
ndef – number of different overlap factors.
kdef – size of the NPF array, containing the description of all overlap factors.
KPF(i) – number of determinants in the ith overlap factor, kd.
IPF(i) – pointer on the ith overlap factor in array NPF, ip.
NPF(ip+1:ip+kd) – contains description of the ith overlap factor as NPF(.) = ipd*bf + nd, where bf = 16 –
packing basis for overlap factors; ipd – pointer to the overlap determinant; nd – its power.
29. C,ij,int,idf repeat up to the end of file
C angular coefficient for matrix element between angular symmetries i and j, recorded as ij=i × bc +j,
where bc = 215 – packing basis for configurations.
int pointer on the relevant radial integral in the packing form
int = m × b8 +k × b4 +i1 × b3 +i2 × b2 +i3 × b+i4 , where m – type of integral; k – multipole index;
i1 ,i2 ,i3 ,i4 – pointer on the involved orbitals; b = 10 – packing basis for orbital indexes.
idf pointer on the corresponding overlap factor.

Remarks.

• Only one type of transition can be in the given mult_bnk.

• The γ  L S  O[λ] γ LS reduced matrix elements can be obtained from the γ LS O[λ] γ  L S   matrix elements stored in
 
mult_bnk by using the factor (−1)L+S−L− S .
• The present package assumes the L + S coupling order. If one uses atomic wave functions with the S + L coupling order (as
in the MCFH_ASP package by Froese Fischer et al. [52], for example), one needs to include the correction factors (−1)L+S−J
in both sides.
• The reduced matrix elements in the mult_bnk are presented as a linear combination of radial integrals, but not as a summation
over weighted one-electron reduced matrix elements, as is done in the MCHF_MLTPOL program [56]. The user can switch
mult_bnk to the MCHF representation by changing the value of the internal parameter mltpol = 0 to mltpol = 1.

10. Program BSR_DMAT

The BSR_DMAT program prepares the dipole matrix, [Link], for subsequent use in the BSR_PHOT program for photoioniza-
tion calculations, or performs the calculations of oscillator strengths between B-spline bound-state solutions. For completeness,
BSR_DMAT also has the option to calculate oscillator strengths between MCHF solutions, and between B-spline bound states and
MCHF solutions.

10.1. Dipole transition matrix

For photoionization calculations we need the D-files, which contain the dipole matrix elements between the initial state and the
R-matrix solutions for a given partial wave. It is convenient to represent the dipole matrix elements by using the dipole transition
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 317

matrix, which consists from the dipole matrix elements between the basis states. Suppose we have two sets of solutions of the
Hamiltonian matrix in different bases. In our case, it can be either the configuration basis from the MCHF calculations, or a close-
coupling B-spline expansion (7.2) from a BSR calculation. Then the matrix element for the given solutions a and b is


a|E [λ] |b = ai dij bj = a T db, (10.1)
i,j

where {ai } and {bi } are the corresponding expansion coefficients in the respective bases. Once we generate the dipole transition
matrix d, the matrix elements for any other solutions is simply defined by convoluting d with the corresponding expansion vectors.
The logical structure of BSR_DMAT is such that it first generates the relevant dipole transition matrix, and then, if indicated,
calculates the corresponding set of oscillator strengths. This enables us to use the same code for photoionization and bound-state
calculations.
Consider the structure of the transition matrix in the B-spline representation. This structure is very similar to the structure of the
Hamiltonian matrix considered in the write-up of the BSR_MAT program in Section 7. In the present close-coupling approach, the
B-spline bound-state solution vector, c, can be written as
c = [a 1 , . . . , a Nc , b]T , (10.2)
where each a i is a column vector of B-spline coefficients for the given channel functions Φ̄iΓ
a i = [a1i , a2i , . . . , ans i ]T (10.3)
and b is the column vector of correlation functions coefficients
bj = [b1 , b2 , . . . , bNp ]T , (10.4)
where Nc is the number of channels while Np is the number of correlation functions. Then the transition matrix d has the form
 
D(11) D(12) . . . D(1Nc ) d(1p)
 D(21) D(22) . . . D(2Nc ) d(2p) 
 
 
 .. .. .. .. .. , (10.5)
 . . . . . 
 
 D(Nc 1) D(Nc 2) . . . D(Nc Nc ) d(Nc p) 
d(1p)T d(2p)T . . . d(Nc p)T d(pp)
where d(pp) is an Np × Np matrix that comes from the bound–bound interaction, the d(ip) are the ns × Np matrices, representing
the interaction between the ith channel and the bound states, and the D(ij ) are the ns × ns matrices for the channel–channel
interaction. When we work with MCHF solutions, the transition matrix is reduced to the single last row, if the initial state is given
in the MCHF representation, or simply to the d(pp) matrix, if both solutions are MCHF configuration expansions.
Each individual matrix element in (10.5) is expressed by the program MULT in the form

 
ci R λ (a, b) × D {nl}; {n l  } , (10.6)

where the ci are numeric coefficients that only depend on the angular symmetry of the CSFs involved, R λ stands for the correspond-
ing one-electron dipole integrals (9.20), and D({nl}, {n l  }) is the overlap factor which depends only on the orthogonality of the
radial orbitals used in the construction of the basis CSFs. In general, the overlap factor is the multiplier of determinants of matrices
consisting of one-electron overlaps between radial functions in the initial and final states (see Eq. (7.8)).
The B-spline orbitals can appear in the radial integrals R λ , as well as in the determinant factors. To derive the final expression
for the matrix elements in the B-spline basis, let us first simplify the overlap factors (7.8) by expanding them over those rows which
contain one-electron overlaps with the B-spline orbitals. At most, there can be two such rows. The residual overlap determinants
depend only on the known bound orbitals and they can be calculated in a standard manner. It is convenient to redefine the angu-
lar coefficients ci by multiplying them by these residual overlap factors. As a result, the initial overlap factor D({nl}, {n l  }) in
expression (10.6) is reduced to one-electron or two-electron overlaps of the form
kl|nl, kl|k  l or kl|nln l  |k  l  , (10.7)
where |kl stands for the radial function of the B-spline orbital


ukl (r) = ai Bi (r). (10.8)
i
Each term in (10.6) gives rise to some matrix in the spline basis. To elucidate the structure of these matrices, we introduce the
following ns × ns matrices with elements
B(..)ij = Bi |Bj , R(..)ij = Bi |r|Bj . (10.9)
318 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

Table 2
Contribution of different terms (without indication of angular coefficients) to the transition matrix. The symbols a, b stand for the bound orbitals, while the symbol ki
indicates the B-spline orbital in channel i
Term Contribution Remarks
Channel–channel interaction
R(ki , kj ) R(··) Banded matrix
ki |ab|kj R λ B(·a) · B(·b) ∗ R λ (or L) Full matrix,
contribution due to non-orthogonality
ki |kj R λ B(··) ∗ R λ Banded matrix
R λ (ki , a)b|kj  L(·a) · B(·b) Full matrix,
contribution due to non-orthogonality
Channel–bound interaction
R λ (ki , a) R λ (·a) Vector
ki |aR λ B(·a) ∗ R λ Vector,
contribution due to non-orthogonality
Bound–bound interaction
R λ (a, b) Scalar

To describe the channel–bound interaction, we introduce the following vectors with elements


ns

ns
B(.a)i = aj Bij , R(.a)i = aj Rij . (10.10)
j =1 j =1

Finally, for bound–bound interaction we have

ns
R(ab) = ai Rij bj . (10.11)
i,j =1

With this notation, Table 2 represents the contributions of the different terms from the matrix element expression (10.6), which
usually produces the angular integration codes, to the transition matrix (10.5), along with an indication of its structure.

10.2. Oscillator strengths and decay probabilities

For completeness of the description, we provide below the main formulae used in the present code for the calculation of radiative
properties. The line strength for a given transition is defined as

 
Sij = Sj i = Ji Mi |O [λ] |Jj Mj 2 . (10.12)
µ
m s

It is given in atomic units, which are e2 a02λ and β 2 a02(λ−1) for electric and magnetic transitions, respectively. Here β = eh̄/2mc is
the Bohr magneton. The transition probability can be written as:
 
2(2λ + 1)(λ + 1) 1 1 Ej − Ei 2λ+1
Aj i = Sij . (10.13)
λ[(2λ + 1)!!]2 gi h̄ h̄c
The oscillator strength (unitless) is related to the transition probability by
mc mc3 gj Aj i
gi fij = g j Aij λ 2
= , (10.14)
8π 2 e2 2e2 ω2
independent of the nature of the radiative transition.

10.3. Structure and data flow

The block diagram of the program BSR_DMAT, along with the data flow, is given in Fig. 9. Since BSR_DMAT requires only
a small number of input parameters, all of them must be defined in the command line as arguments. The input parameters for the
BSR_DMAT program are the names of the c-files with a list of configuration state functions for the initial and final states involved,
and the mode of calculation ctype (see Section 10.5). The input c-files may be the ordinary c-files from MCHF calculations, or
[Link] files generated in BSR calculations. The mode of calculation defines whether this is the generation of the transition matrix
for subsequent photoionization calculations or the calculation of oscillator strengths between bound states. The ctype parameter
also defines the representation format for the initial and final states, which can be in the LS- or LSJ-coupling schemes. BSR_DMAT
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 319

Fig. 9. Block diagram for the program BSR_DMAT and data flow (see text).

reads a set of other input data files. All of them have default names, or names which are directly connected to the names of the input
c-files.
First, the program defines the set of LS terms and one-electron orbitals under consideration (routines read_term and read_orbit).
The orbitals in the initial and final sets are supposed to be non-orthogonal. Then the program reads the one-electron radial functions
in the B-spline representation (routine read_bsw). They can be either the user-prepared [Link] files, or the [Link] files
generated by the BSR_HD program. In the latter case, the program also reads the target file with a description of the relevant
close-coupling expansions. Finally, the program generates all the B-spline arrays needed in the routine define-spline. Note that all
arrays involved are automatically allocated in accordance with the input data.
In the next step, routine check_bnk, the program analyses the information contained in the databank mult_bnk, and determines
whether the databank contains all angular coefficients needed for the given case. If not, the program stops, suggesting to rerun the
MULT program. Note that the multipole index and the type of transition under consideration are defined in mult_bnk.
The routine D_matr controls the calculation of the transition matrix. First, the routine reads all overlap factors D({nl}, {n l  })
from the databank in routine read_dets before calculating all necessary R λ integrals in routine quadr. All relevant B-spline matri-
ces presented in Table 2 are then generated by the routine bxv from the B-spline library. The list of specific integrals (10.6) for the
configurations under consideration is generated dynamically from the information in mult_bnk. The overlap factor D({nl}, {n l  })
in expression (10.6) is reduced to one-electron or two-electron overlaps of the form (10.7) by the routine check_det. The integrals
are then sorted according to their structure (see Table 2) and stored in the module cmdata in different ordered blocks (routines
add_coef and add_cdata). When all coefficients from mult_bnk are exhausted, or when all blocks are full, they are merged if nec-
essary with the routine merge_cdata before being processed by gen_matrix to get the B-spline representation for specific integrals
(routine d_data). The corresponding contributions are then added to the interaction matrix, located in module dmatrix.
The next step of the calculation depends on the input mode. In photoionization calculations, the dipole matrix (10.1) is calculated
on the base of the R-matrix solutions placed in the file [Link], and then is recorded in the file [Link]. For calculations of oscillator
strengths, however, the program calls the routine LS_case, or J_case for processing transitions between LSJ states. The resulting
line strengths, oscillator strengths, transition energies, and transition probabilities are recorded in the file zf_res (see Section 10.7
below).
320 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

10.4. Data files

name1.c or [Link] File type: formatted sequential input.

Created by MCHF or BSR_CONF programs.
Read by routine check_bnk.
Description: contains the configuration expansion for the initial state.
name2.c or [Link] File type: formatted sequential input.
Created by MCHF or BSR_CONF programs.
Read by routine check_bnk.
Description: contains the configuration expansion for the final state.
mult_bnk File type: unformatted sequential input.
Written by program MULT.
Read by routine check_bnk.
Description: databank for the angular coefficients of the dipole operator.
[Link] File type: formatted sequential input.
Created by program BSR_HD.
Read by BSR_DMAT program.
Description: bound-state solutions for a given partial wave in the B-spline representation (10.2).
target File type: formatted sequential input.
Written by user and modified by BSR_PREP and BSR_CONF programs.
Read by routine read_target.
Description: contains description of the target states and scattering channels.
[Link] File type: formatted sequential input.
Written by user.
Read by routine define_grid from BSPLINE library.
Description: input parameters that define the B-spline grid.
[Link] File type: unformatted sequential input.
Created by program bsr_prep.
Read by routine read_bsw.
Description: target one-electron orbitals in the B-spline basis.
pert_nnn.bsw File type: unformatted sequential input.
Created by program bsr_prep.
Read by routine read_bsw.
Description: perturber one-electron orbitals in the B-spline basis, optional.
[Link] File type: unformatted sequential input.
Created by user.
Read by routine read_bsw.
Description: one-electron orbitals in the B-spline basis for the initial state, optional.
[Link] File type: unformatted sequential input.
Created by user.
Read by routine read_bsw.
Description: one-electron orbitals in the B-spline basis for the final state, optional.
[Link] File type: formatted sequential output.
Written by program BSR_DMAT.
Read by BSR_PHOT.
Description: dipole matrix.
zf_res File type: formatted sequential output.
Written by program BSR_DMAT.
Read by user.
Description: oscillator strengths for the dipole operator given in mult_bnk.
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 321

bsr_dmat.log File type: formatted sequential output.

Written by program BSR_DMAT.
Read by user.
Description: running information.
10.5. Input parameters

Input parameters must be provided in the command line.

name1.c or [Link] name of c-file for the initial state.

name2.c or [Link] name of c-file for the final state.
ctype1 character*1, denotes the representation mode for the initial state:
c – initial state is defined by c-file from MCHF calculations.
l – initial state is defined by l-file from MCHF calculations.
j – initial state is defined by j -file from MCHF calculations.
b – initial state is defined by [Link] file from BSR calculations.
ctype2 character*1, the same as ctype1, but for the final state. An additional mode is p, which indicates
that only the dipole transition matrix should be generated in the output file [Link] for subsequent
photoionization calculations. All other modes assume the calculation of oscillator strengths, with
the output put in the file zf_res.

10.6. Output of dipole matrix

Summary of the output records in the [Link] file.

2. LT1,ST1,IP1,E1 – term and energy of the initial state.

3. LT2,ST2,IP2,nstate2 – term and number of final states.
4. (cl(i),cv(i),i=1,nstate2) – dipole matrix in length and velocity form.

10.7. Output of radiative data in zf_res

The output format in zf_res is the same as in the MCHF Atomic Structure Package [52]. An example of output is given in
Fig. 10. Each transition is represented by a block, which consists of two empty lines as a delimiter, two lines with a description
-1 -37.77853956 1s(2).2s(2).2p(2)3P2_3P
-1 -37.48150929 1s(2).2s_2S.2p(3)2D3_3D
65187.63 CM-1 1534.03 ANGS(VAC) 1534.03 ANGS(AIR)
E1 S = 4.85375D+00 GF = 9.61141D-01 AKI = 1.81637D+08
6.44181D+00 1.27561D+00 2.41066D+08

-1 -37.48150929 1s(2).2s_2S.2p(3)2D3_3D
-1 -37.45433756 1s(2).2s(2).2p_2P.3p_3P
5963.23 CM-1 16769.43 ANGS(VAC) 16764.85 ANGS(AIR)
E1 S = 6.26544D+00 GF = 1.13495D-01 AKI = 2.99142D+05
1.17451D+00 2.12756D-02 5.60767D+04

2 -940.13842975 2p(2)3P2
2 -932.54958058 2p(4)3P2
1665544.77 CM-1 60.04 ANGS(VAC) 60.04 ANGS(AIR)
E2 S = 3.55924D-05 GF = 2.76118D-08 AKI = 1.70307D+04

0 -940.46450303 2p(2)3P2
2 -940.13842975 2p(2)3P2
71564.16 CM-1 1397.35 ANGS(VAC) 1397.35 ANGS(AIR)
M1 S = 1.80864D+00 GF = 5.23415D-06 AKI = 5.96024D+03
Fig. 10. Example of the output zf_res file.
322 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

of the initial and final atomic states, one line with the transition energy in cm−1 and Angstroms, and finally a line with the line
strength, the oscillator strength, and the transition probability in the length form. For E1 transitions, there is an additional line with
the results in the velocity form.
Each atomic state is represented by its spectroscopic label, its absolute energy in a.u., and the value of 2J , in case of a transition
between LSJ states. The f -values are always for transitions from the lower to the higher state, whereas the transition probabilities
are always reported from the higher to the lower state.
The file zf_res is opened in the BSR_DMAT program with access mode append. Hence, the results from different calculations
can be accumulated.

11. The program BSR_PHOT

The BSR_PHOT program runs the photoionization calculations for one partial wave, based on the [Link] and [Link] files prepared
by the programs BSR_DMAT and BSR_HD, respectively. The present implementation of the R-matrix method for photoionization
calculations differs from the previous one in that the initial state may be calculated independently from the scattering solutions,
with any amount of non-orthogonality between the initial and final one-electron radial functions. We assume that the initial state is
well confined to the inner region.

11.1. Structure and data flow

The block diagram of the program BSR_PHOT, along with the data flow, is shown in Fig. 11. All input/output files have
prescribed internal names (see Section 11.2). First, the program reads the input parameters from the file bsr_phot.inp. The input
electron energies are given in Ry and defined in blocks with specific initial energy, energy step, and final energy. The user can define
up to ten blocks for a given run that allows to calculate the cross sections with different ‘intensity’, i.e. some energy regions can be
investigated in more detailed, with small energy increments. All output files are opened in the position mode APPEND, and hence
all sequential results are added to the output without deleting results from previous calculations. This allows the user to accumulate
data from different runs.
The radiative data for a given case are defined in the file [Link], together with the terms for the initial and final states, as well as
the energy of the initial state. Based on the term of the final state, the program looks for the corresponding scattering and R-matrix
data in the file [Link], which may contain data only for one partial wave. Finally, the program reads the internal solution weights

Fig. 11. Block diagram for the program BSR_PHOT and data flow (see text).
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 323

from the file [Link]. This optional input is guided by the input parameter nwt (see Section 11.3). These data make it possible
to define the composition of the final total wavefunction for a given energy, which may be used for detecting and classifying of
possible resonance structure.
At the next stage the program executes fully independent calculations of the photoionization cross section for each input energy.
The algorithm is given in Section 2.4. In order to determine the photoionization cross section, we need the set of linearly independent
solutions (2.33) in the external region. The determination of these solutions is beyond the scope of the present package, which deals
with the internal region solutions in the general case. For the treatment of the external region, we currently choose the program
ASYPCK [43]. This is a general program, which can determine the solution for both neutral and ionized atoms. It works well over
a wide range of energies, including near-threshold regions. ASYPCK is called as a subroutine, which is convenient for its use by
another program.
The input data for the external region consist only of the description of the scattering channels and asymptotic coefficients (2.26).
Only four main parameters, which are given in Section 11.3, guide the execution of ASYPCK. For a more detailed description of
the ASYPCK package, the user is referred to its write-up and that of the program IMPACT [36]. The ASYPCK program is called in
the present calculation by using the interface subroutine zaface (see Fig. 11). We did not introduce any modifications to ASYPCK
itself, and hence the user should check if the dimension parameters in ASYPCK are big enough for a given run. Over the past
years, a set of more effective algorithms were developed for the treatment of the external region, such as the flexible asymptotic
R-matrix package FARM [28] or a new version of the program STGF [59]. In principle, ASYPCK can easily be replaced by any
other algorithm, provided the latter is organized as a subroutine call.
The photoionization calculations for a given energy consist of the following steps. The routine zopen defines the number of open
channels and the corresponding channels energies. Then routine zaface calls the program ASYPCK to determine the asymptotic
radial functions (2.27) at the R-matrix radius. The routines zrmat, zkmat and get_eph calculate the R- and K-matrices, together
with the eigenphase sum at the energy of interest. Then the total dipole matrix elements (2.31), and the partial and total cross
sections are calculated by the routine [Link]. Optionally, the routine ak_coef is called if one needs the channel composition for
the total wave function.
The results are saved in a set of output files. The file [Link] contains the total cross section and the accuracy of the calcu-
lations, defined via the average asymmetry of the K-matrix. The results in this file are sorted according to energy by the routine
sort_photo. The file [Link] is in simple table form and can be used for a quick check of the calculations. The full results
including partial cross sections, eigenphase sums, channel weights and full dipole matrix elements (2.31) are recorded in the file
bsr_phot.nnn. These data may then be used for the calculation of differential cross sections or the asymmetry parameter β (see the
description of utility-program photo_tab in Section 13).

11.2. Data files

bsr_phot.inp File type: formatted sequential input.

Created by the user.
Read by the routine read_arg.
Description: input parameter for given run.
[Link] File type: unformatted sequential input.
Created by BSR_DMAT program.
Read by the routine read_Ddat.
Description: radiative dipole transition matrix for given initial state and final scattering partial wave.
[Link] File type: unformatted sequential input.
Created by the BSR_HD program.
Read by routine read_Hdat.
Description: diagonalized Hamiltonian matrix data for the inner region.
[Link] File type: unformatted sequential input.
Written by the program BSR_HD.
Read by the routine read_Wdat.
Description: weights of different channels in the inner-region solutions.
[Link] File type: formatted sequential output.
Created by the program BSR_PHOT.
Read by the user.
Description: electron and photon energies, total cross-section in the length and velocity forms, eigenphase
and accuracy of the calculations, defined as the average asymmetry of the K-matrix.
324 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

bsr_phot.nnn File type: formatted sequential output.

Created by the program BSR_PHOT.
Description: partial cross sections, eigenphases, dipole matrix elements (2.31) and channel weights for each
input energy.
bsr_phot.log File type: formatted sequential output.
Created by the program BSR_PHOT.
Read by the user.
Description: running information.

11.3. Input parameters

The input parameters are provided in the file bsr_phot.inp. All input data have the format ‘name = value’, and should be placed
at the beginning of the line; otherwise they will be ignored. The default values of the data are indicated in the brackets. All default
values, along with the unit numbers and the default file names are placed in the module mod_phot. Some of the input parameters
deal with the asymptotic package ASYPCK.

klsp index of partial wave (nnn).

elow Define the range of input electron energies (in Ry) as the initial energy (elow), the energy increment (estep),
estep and the final energy (ehigh). This can be given for up to ten different input energy intervals.
ehigh
awt [100] atomic number, used to convert Ry to eV.
ibug [0] defines the level of debug printing in the file bsr_phot.out.
e_exp [0.0] if = 0, define the ionization (detachment) energy.
nwt [0] if > 0, then the weights of each channel in the total solutions are also calculated, based on the data in the
file [Link].
Parameters for ASYPCK program:
AC [0.001] accuracy for the asymptotic solutions.
DR [0.1] in ASYPCK, the solutions are first obtained at R1=R-DR and R2=R+DR, and then by interpolation – in R.
mfg [300] an asymptotic expansion of the solution will be sought at R + DR*MFG.
iauto [2] if > 0, then ASYPCK may automatically increase parameter MFG.

12. Libraries

The BSR package includes the four libraries BSPLINE, ZCONF, ZCOM and LAPACK, which contain routines common
to all programs. The ZCONF library contains specific routines, which deal with the configuration and close-coupling expan-
sions, while ZCOM contains common routines used for a variety of auxiliary computations. The LAPACK library ([Link]
[Link]/lapack/[Link]) was chosen to provide the necessary linear algebra calculations, such as diagonalizing matrices or
solving a system of linear equations.
The BSPLINE library is the central part of the present complex. It contains routines for using B-spline bases in atomic structure
calculations. Extensive survey of recent applications of B-splines in atomic and molecular physics is presented in the review by
Bachau et al. [26]. The present library represents the accumulated developments made by different people over the last twenty years.
BSPLINE includes original FORTRAN IV routines from the monograph by de Boor [60], who first presented algorithms for the
construction and manipulation of B-splines. These were modified by Froese Fischer and her collaborators who started the design
and implementation of the present library, and numerous investigations of different numerical algorithms for using B-splines in
atomic structure calculations were published [39,49,50,61–64]. This spline library has been modified and extended by the present
author. Features were added to accommodate the more complex requirements of the R-matrix method, and the programming was
upgraded to F95 standard.
Most of the routines in BSPLINE can be used in many other atomic structure calculations. For this reason we decided to present
a rather detailed description of the BSPLINE library, including calling sequences for the main routines and the internal structure of
the modules. The FORTRAN module is a program unit, which contains specifications and definitions that can be made accessible to
other program units. It was found convenient to cluster the principal parameters, variables and arrays in different modules according
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 325

Fig. 12. The full set of B-splines of order k = 8 relative to the equidistant knot sequence (represented by full circles, with an indication of their numeration).

to the subject they describe. For the given module name to be accessible, the other program units must reference the module in a use
name statement. All arrays in the modules in the present implementation are dynamically allocated for the case under consideration.
The module structure used in BSPLINE simplifies the use of separate subroutines in other applications.
The modules and routines in the description below have been grouped together for convenience. This includes the general
routines and the atomic-orbital-related routines. The latter concern first of all the evaluation of a two-electron integral in the B-
spline basis, and the generation of the B-spline representation for the Hamiltonian matrix. Along with the description of the routines,
each section provides a brief description of the basic theory.

12.1. Description of B-splines

B-splines are functions designed to generalize polynomials used for approximating arbitrary functions in some finite interval.
Consider the interval [a, b] divided into subintervals. The endpoints of these subintervals are given by the knot sequence [ti ], i =
1, 2, . . . , n + k. The B-splines of order k, Bi,k (r), on this knot sequence are defined recursively by the relations

1, ti  r  ti+1 ,
Bi,1 = (12.1)
0, otherwise
and
r − ti ti+k − r
Bi,k (r) = Bi,k−1 (r) + Bi+1,k−1 (r). (12.2)
ti+k−1 − ti ti+k − ti+1
Since k and t are usually fixed, we will indicate only one index, Bi . Each Bi is defined over the interval [ti , ti+k ], which contains
k + 1 consecutive knots, and is indexed by the knot where it starts. An example of a full set of B-splines is shown in Fig. 12. In this
example, the interval [0, 5] is divided into equidistant subintervals, with h = 0.5. Each Bi starts at the knot ti , i = 1, . . . , n, and ends
k knots later. In order to provide the same number of B-splines in each interval, the multiplicity of knots at the endpoints is usually
chosen as the maximum possible value, which is equal to the B-spline order k. The most common choice for the multiplicity at the
inner knots is unity, corresponding to the maximum continuity of B-spline functions inside the interval. With this choice, employed
in the BSPLINE library, the number of subintervals is related to the number of B-splines as
nint = n + 1 − k.
Some general properties of B-splines are the following:

• Over each subinterval [ti , ti + 1], exactly kB-splines are nonzero

Bj (x) = 0 for j = i − k + 1, . . . , i,
the first being Bi−k+1 , which ends at ti + 1, and the last being Bi , which starts at ti . This means that
Bi (x) · Bj (x) = 0 for |i − j |
k.
• In the expansion of an arbitrary function


n

i
f (x) = cj Bj (x) = cj Bj (x) for x ∈ [ti , ti+1 ]
j =1 j =i−k+1
326 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

only k terms contribute. Consequently, the B-spline representation for different operators will exhibit a band structure, simpli-
fying to a great extent the solution
 of the relevant matrix equations.
• B-splines are normalized as i Bi (x) = 1. Since the B-splines are non-negative with minimal support, the expansion coeffi-
cients of an arbitrary function are close to the function values at the knots. This means that wild oscillations in the coefficients
are avoided, cancellation errors are minimal, and the numerical stability is maximized.
• B-splines satisfy the recursion relations (12.2). Together with the definition of B-splines of order k = 1 (12.1), this provides a
simple algorithm for their practical evaluation: given a point x, one generates by recursion the values of all kB-splines that are
nonzero at x. The derivative of a B-spline of order k can also be expressed as a linear combination of B-splines of order (k − 1):
k−1 k−1
DBi,k (r) = Bi,k−1 (r) − Bi+1,k−1 (r). (12.3)
ti+k−1 − ti ti+k − ti+1
The first B1 -spline is the only spline with nonzero value at the left end of interval [a, b] while all other Bi -splines have
(r − a)i−1 dependence. A similar behavior applies at the point b. This makes it very easy to implement boundary conditions
at the endpoints. For instance, f (a) = 0 is satisfied by deleting B1 from the set, or Df (b) = 0 may be satisfied by combining
Bn−1 and Bn as B̃n−1 = c1 Bn−1 + c2 Bn so that D B̃n−1 (b) = 0, and so on.
• With a minimum multiplicity of knots, the Bi,k -spline functions belong to the functions of class C k−2 , i.e. they are all continuous
functions together with their derivatives up to order k − 2. In each subinterval [ti , ti+k ], B-spline functions can be represented
by a polynomial of order k (maximum degree k − 1). That allows, in principle, for analytical manipulations with B-spline
functions in symbolic packages.

In quantum mechanical applications we deal mostly with matrix elements, i.e. definite integrals involving B-splines and their
derivatives. These are best computed to machine accuracy by employing Gauss–Legendre integration of appropriate order over each
subinterval [ti , ti+1 ]. Recall that the Gauss–Legendre method of order n (n points) is exact for a polynomial integrand of order 2n.
This means that a k-points Gauss–Legendre method integrates exactly the product Bi Bj , and also provides accurate results even
for non-polynomial integrands, provided no singularity is present in the given interval. At the same time, in the first subinterval,
(1/r)Bi Bj is also evaluated exactly if i or j > 1. For most applications, therefore, we only need the B-spline values at the k Gauss
points. This also spares us from directly dealing with the vicinity of r = 0, where the potentials may have singularities.
Another important problem in quantum mechanical applications is the solution of the radial Schrödinger equation. In order to
demonstrate the use of B-splines for this problem, consider the one-electron radial equation

d2 l(l + 1)
− + + V (r) unl (r) = Enl unl (r). (12.4)
2dr 2 2r 2
An approximation solution is expressed as a linear combination of B-splines over an interval [0, rmax ], say {Bi (r), i = 1, . . . , ns }.
Then


ns
unl (r) = cinl Bi (r). (12.5)
i=1

Usually the approximate (or trial) functions satisfy the boundary conditions. When the B-spline basis is generated using a grid with
multiple knots at the ends of the interval, the zero boundary conditions, for example, are fulfilled by removing the first and the last
B-splines from the basis set. We shall refer to the column vector c = {ci } as a spline expansion vector of u(r). The solutions of
(12.4) are calculated by solving the system of ns linear equations obtained when substituting (12.5) into (12.4) and projecting onto
Bj (r). Written in matrix form, this procedure is equivalent to solving the generalized eigensystem

H · c = ES · c, (12.6)
where
1 d2 l(l + 1) 1
Hij = − Bi | 2 |Bj  + Bi | 2 |Bj  + Bi |V |Bj , (12.7)
2 dr 2 r
Sij = Bi |Bj . (12.8)
The overlap matrix S originates from the fact that the B-splines do not form an orthogonal set of basis functions. We will refer to
the Hij matrix as the B-spline representation of the Hamiltonian. Note that all matrix elements in (12.6) may be easily computed, as
discussed above, to machine accuracy using Gauss–Legendre integration, applied to each subinterval. Another key property related
to the use of a B-spline basis follows from the fact that a B-spline of order k differs from zero only on k successive subintervals.
Therefore, any local operator expressed in the B-spline basis appears as a banded matrix with nonvanishing values only on a
diagonal band of width 2k − 1. This simplifies considerably the solution of the generalized eigenproblem (12.6).
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 327

Provided the wavefunctions fit into the cavity (often called the box), in which we solve the problem (12.6), the first few eigen-
values and eigenvectors agree precisely with the corresponding bound-state eigenvalues and eigenvectors obtained by numerically
integrating the radial Schrödinger equation, but as the principal quantum number increases, the cavity spectrum departs from the
real spectrum and forms the so-called pseudospectrum, which can be used instead of the real spectrum in different physical prob-
lems. Each of these eigenvectors defines an atomic orbital. By using a B-spline set, it is possible to account for the influence of
an entire Rydberg series without introducing a cut-off. For example, the second-order energy correction in the framework of the
Rayleigh–Schrödinger perturbation theory

0|V |nn|V |0

E (2) = (12.9)
En − E0
n=0

requires the direct evaluation of the sum in terms of an infinite sum over the bound states and an integration over the continuum.
In the B-spline pseudospectrum, only a few states are bound while the remaining states are unbound, and the infinite sum (12.9) is
reduced to a finite sum. As shown in many numerical examples [65], this sum very accurately reproduces the infinite sum (12.9) with
a small number of B-splines, i.e. the B-spline basis is effectively complete in this sense. This effective completeness of the basis set
means that the summation is independent, or almost independent, of the box size and the number of splines if these two quantities
are chosen sufficiently large. We can check numerically if the summation is independent of the box size and the number of splines.
Completeness of the B-spline basis is a very valuable property, simply because it is possible to replace an infinite summation over
bound states and an integration over continuum states with a summation over a small number of L2 -integrable functions. This
property also found wide application in the investigation of the strong interaction between a perturber and a Rydberg series [66], or
in the investigation of the influence of the continuum spectrum in two-electron atoms [67].

12.2. General routines

B-splines can be considered as elementary functions like sin(x) or jl (x). One becomes familiar with them by understanding
their qualitative behavior and how to use them, i.e. how to obtain values of the functions, their derivatives, or integrals. This is the
purpose of the routines described in the next section.

List of general routines

bav computes y = Bx or y = xB, where B is symmetric or non-symmetric band matrix, y and x are vectors.
bvalue calculates at given radius the value and derivatives of a function given in a B-spline representation.
bvmv returns the matrix element x|A|y, where x, y are vectors, A is a band array in symmetric, non-symmetric,
or almost symmetric form.
bxv computes y = Bx, where B is a symmetric band matrix while y and x are vectors.
define_grid gets input data for the grid and sets up the knots for spline.
define_spline initializes the values of the spline and its derivatives, and evaluates the basic elementary operators in the
spline basis.
dinty computes the matrix elements Bi |y(r)|Bj  between individual B-splines for an input function y(r) given
at the Gaussian points of each interval.
facsb factorizes the overlap matrix Bi |Bj .
gauss provides the Gaussian coordinates and weights over the interval [0, 1].
minty computes the matrix elements Bi |y(r)|Bj  between individual B-splines for an input function y(r) given
at the Gaussian points of each interval.
mkgrid sets up the knot sequence.
mrm generates the B-spline representation of the operator r m .
spline_param module, contains the main spline parameters.
spline_grid module, defines the values of the splines and their most general functions at the Gaussian points needed for
integrations on the grid.
328 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

spline_galerkin module, contains common arrays used in the application of splines in the Galerkin method.
splin3 provides coefficients for a cubic interpolating spline.
vbsplvd calculates the values of the B-splines and their derivatives for one point per each interval.
vinty computes the vector elements Bi |y(r), where the function y(r) is given at the Gaussian points of each
interval.
yval computes the values of r m f (r), r m f  (r) or r m f  (r) at the Gaussian points of each interval, where f (r) is
defined by the spline expansion vector.

Internal representation of band matrices

The BSPLINE library uses different storage modes for banded matrices depending on the symmetry of the matrix and on the
operator represented by this matrix. Let A(n, n) be a matrix with a diagonal band of width 2k − 1. Then the possible storage modes,
matrices B, can be described schematically as:

• Symmetric lower-column storage mode

   
1   1
   
• 1   • 1
   
A(n, n) → Aij =  k • 1 , B(n, k) → Bi,j  = k • 1,
   
 k • 1  k • 1
   
k • 1 k • 1
 
A(i, j ) = B(i, j ) with j = j + i − k.

The summation over elements: Do j’= 1,k

Do i = k+1-j’,n
• Symmetric upper-column storage mode
   
1 • k  1 • k
   
 1 • k  1 • k
   
A(n, n) → Aij =  1 • k, B(n, k) → Bi,j  = 1 • k,
   
 1 • 1 • 
   
1 1
A(i, j ) = B(i, j  ) with j = j  + i − 1.

The summation over elements: Do j’= 1,k

Do i = 1,n-j’+1
• Non-symmetric column storage mode
   
1 • k   1 • k
   
• 1 • k   • 1 • k
   
A(n, n) → Aij =  k • 1 • k , B(n, 2k − 1) → Bi,j  =k • 1 • k,
   
 k • 1 • k • 1 • 
   
k • 1 k • 1
A(i, j ) = B(i, j  ) with j = j  + i − k.

The summation over elements: Do j’= 1,k+k-1

Do i = max(1,1+k-j’),min(n,n+k-j’)
• Almost symmetric lower-column storage mode
There is only one non-symmetric element: A(n, n − 1) = A(n − 1, n), the latter is saved in B(1, 1) for low-column storage
mode, or in B(n, n) for upper-column storage mode.

MODULE spline_param

ks order of B-splines.
ns number of B-splines.
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 329

nv number of intervals.
ml number of distinct points from 0 to 1.
me number of points in the “exponential” region.
h initial step in knot sequence.
hmax maximum step.
rmax border radius.

MODULE spline_grid

t(1:ns+ks) knot sequence.

bsp(1:nv+1,1:ks,1:ks) bsp(i,m,ith) contains value of the i+ith-1 B-spline in interval i at Gaussian
point m.
bspd(1:nv+1,1:ks,1:ks,1|2) values of first and second derivatives of B-splines at the same points.
bsp(nv+1,1,.) and bspd(nv+1,1,.,.) contain the corresponding values at
the last knot point, i.e. at rmax .
bsq(1:nv+1,1:ks,1:ks) values of the expression B(r) − 1/rB  (r) at the Gaussian points.
gr(1:nv,1:ks) gr(i,m) contains Gaussian point m in interval i.
grm(1:nv,1:ks) reciprocals of the values of gr(i,m).
grw(1:nv,1:ks) Gaussian weights at the corresponding points.
Internal routines:
Allocate_grid allocates space for arrays bsp, bspd, gr, grm, and grw according to parameters
from module spline_param.

MODULE spline_galerkin

sb(1:ns,1:ks) sb(i,j) -> <Bi |Bi+j-k >

r1(1:ns,1:ks) r1(i,j) -> <Bi | r|Bi+j-k >
rm1(1:ns,1:ks) rm1(i,j) -> <Bi |(1/r)|Bi+j-k >
rm2(1:ns,1:ks) rm2(i,j) -> <Bi |(1/r2 )|Bi+j-k >
rm3(1:ns,1:ks) rm2(i,j) -> <Bi |(1/r3 )|Bi+j-k >
db1(1:ns,1:ks) db1(i,j) -> <Bi |B’i+j-k >
db2(1:ns,1:ks) db2(i,j) -> <Bi |B”i+j-k >
bs(1:ks,1:ns) factorization of array sb
Internal routines:
Allocate_galerkin allocates space for all arrays according to the parameters from module
spline_param. The arrays are stored in the symmetric lower-column storage mode,
except for db2 – almost symmetric mode, and db1 – antisymmetric array.

Remarks. All arrays and parameters in the above modules should be initialized before the B-spline applications. This can be done
with the routines define_grid and define_spline, which should be called at the beginning of each BSR program.
The define_grid subroutine generates the knot sequence. This is a very important first step in the BSR calculations, because the
B-splines are fully defined by the given knot sequence. Hence the grid, together with the spline order, k, determines the accuracy
of the approximation. Though there is complete freedom in choosing the mesh of knots, the optimal choice will depend on the
type of result we are particularly interested in. With splines, we can create any composite grid if needed but too dense a grid can
lead to a rapid saturation of the computational resources. Experience shows that the most appropriate choice is a logarithmic grid
that reflects the exponential behavior of atomic orbitals. On the other hand, in the continuum calculations, the wavelength of the
scattering particle cannot be smaller than the grid step; otherwise the B-spline basis hardly describes the oscillating behavior of the
wavefunction. In the present implementation, we use a mixed sequence which can easily be adjusted to the task under consideration.
330 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

We introduce four main parameters: the step size at the origin, h = 2−m , where m is the integer, the maximum step size hmax , the
maximum range rmax , and the order of splines k. Let Z be the nuclear charge of the atom. The grid points are now defined through
the array ti = Zri such that
ti = 0 for i = 1, . . . , k,
ti+1 = ti + h for i = k, . . . , k + m,
ti+1 = ti (1 + h) for 1 < ti+1 − ti < hmax , (12.10)
ti+1 = ti + hmax for the rest of ti ,
ti = rmax /Z for i = n, . . . , n + k.
This grid is linear near the origin, at Zr < 1, which guarantees an accurate representation of the radial orbitals at small radii.
Experience shows that values Zh = 0.25 or 0.125 are good enough for most applications. In the middle-radii region, the grid is
exponential. This gives some advantage in computational time, because for B-splines entirely in the outer region, relations such as

Bi , r p Bi  = (1 + h)p+1 Bi−1 , r p Bi−1  (12.11)

hold (see, for example, [39]). Similar rules can also be derived for two-electron integrals commonly occurring in atomic structure
calculations (see Section 12.4). Thus, only a few integrals need to be evaluated while a scaling law applies for others.
Another grid parameter defines the maximum range of the solution. Since the range of a 1s orbital is much less than that of a
10s orbital, the software may automatically adjust the range for each orbital, but in a general purpose software, a sufficient range
is needed to allow orbitals to reach their full extent for the desired accuracy. In BSR calculations, rmax is usually chosen from the
condition that |P (rmax )| < 10−5 for all bound radial functions, and this value coincides with the R-matrix boundary, a.
For bound-state calculations we can use the exponential grid in the outer region, whereas in scattering calculations the maximum
size of the subintervals should be restricted according to the maximum energy of the scattering electron. Experience shows that the
oscillating behavior of the scattering radial functions is reproduced accurately with at least three knot points for a half wavelength.
Hence, the user should choose the parameter hmax  1/kmax , where kmax is the maximum linear momentum of the scattering elec-
tron expected in the subsequent calculations. This restriction is primarily important in direct scattering calculations [42], whereas
8 ==> order of splines (ks)
101 ==> number of splines (ns)
6.00000 ==> nuclear charge (z)
0.25000 ==> step size from 0 to 1 (h for z*r, = 1/2\m)
0.70000 ==> maximum step size (hmax for r)
50.00000 ==> maximum r (rmax)
***
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.04167 0.08333
0.12500 0.16667 0.20833 0.26042 0.32552
0.40690 0.50863 0.63578 0.79473 0.99341
1.24176 1.55220 1.94026 2.42532 3.03165
3.63798 4.24431 4.85064 5.45697 6.06330
6.66963 7.27596 7.88229 8.48862 9.09495
9.70128 10.30761 10.91394 11.52027 12.12660
12.73293 13.33926 13.94559 14.55192 15.15825
15.76457 16.37090 16.97723 17.58356 18.18989
18.79622 19.40255 20.00888 20.61521 21.22154
21.82787 22.43420 23.04053 23.64686 24.25319
24.85952 25.46585 26.07218 26.67851 27.28484
27.89117 28.49750 29.10383 29.71016 30.31649
30.92282 31.52915 32.13548 32.74181 33.34814
33.95447 34.56080 35.16713 35.77346 36.37979
36.98612 37.59245 38.19878 38.80511 39.41144
40.01777 40.62410 41.23043 41.83676 42.44309
43.04942 43.65575 44.26208 44.86841 45.47474
46.08106 46.68739 47.29372 47.90005 48.50638
49.11271 50.00000 50.00000 50.00000 50.00000
50.00000 50.00000 50.00000 50.00000
***
Fig. 13. Example of a [Link] file in the case of the e–C scattering.
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 331

in the R-matrix calculations was found that this restriction can be reduced at least three times without noticeable influence on the
numerical accuracy.
The grid parameters are provided in the BSR calculations in the file [Link], and these parameters should hold the same values
during the entire set of calculations. An example of a [Link] file is given in Fig. 13. The routine define_grid first reads the main
parameters from [Link], and then generates the knot sequence in the routine mkgrid. The number of B-splines is an output
parameter in this procedure and defined at this stage. The resulting knot sequence and the relevant parameters are then rewritten
to the [Link] file. If the program cannot find the [Link] file, the execution is halted. The user is strongly advised to first try
different sets of parameters to produce the desired knot sequence with a minimum number of B-splines, because this parameter
strongly affects the computation time for all subsequent calculations.

12.3. Atomic orbital routines

This block of routines deals with the set of atomic radial functions and their B-spline representations, denoted below as vec-
tors pi . The operator L is defined as
d2 l(l + 1) z
L= 2
− + .
dr r2 r

List of routines

azl provides Bl+2 /r l+1 at r = 0, which defines the leading term for orbital nl.
bcore computes the energy of the common closed shells (core).
bhl evaluates the matrix element pi |L|pj  for atomic orbitals i, j .
bhwf computes the spline expansion coefficients for hydrogen radial functions.
coulomb builds the interaction matrix for the Coulomb problem.
grad gives the matrix element pi |d/dr + l/r|pj > for atomic orbitals i, j .
hlc computes the matrix element pi |L|pj  modified by the interaction with closed shells.
hlm sets up the matrix of the L operator in the B-spline basis.
quadr evaluates the matrix element pi |rk |pj  for atomic orbitals i, j .
mvcv computes the matrix elements for the mass-velocity correction in the B-spline basis.
r_bwfn reads the atomic orbitals from the bsw-files.
spline_atomic module, contains the main atomic parameters.
spline_hl module, contains the L matrix and routines for evaluation of relativistic shifts.
spline_orbital module, contains the description of the atomic orbitals.
zeta evaluates the spin–orbit integrals (including the interaction with closed shells).

MODULE spline_orbitals

nbf number of atomic orbitals.

nbs(1:nbf) principal quantum numbers.
lbs(1:nbf) orbital angular momenta of the orbitals.
kbs(1:nbf) set indexes of orbitals.
ebs(1:nbf) spectroscopic labels for orbitals.
mbs(1:nbf) length of the B-spline expansion.
iech(1:nbf) Additional pointer, used to distinguish the continuum orbitals.
pbs(1:ns,1:nbf) B-spline expansion coefficients for atomic orbitals.
332 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

qbs(1:ns,1:nbf) vectors B · p, where B is the B-spline overlap matrix.

obs(1:nbf,1:nbf) one-electron overlaps.
iort(1:nbf,1:nbf) orthogonal constraints, imposed on the one-electron orbitals.

Internal routines:
allocate_bsorb allocates space for all arrays for a given value of nbf.
Ifind_bsorb defines the position of the orbital n, l, k in the common list.
Iadd_bsorb adds the orbital n, l, k to the common list.

MODULE spline_atomic

z charge of nucleus.
kclosd number of closed shells.
rel indicates a relativistic calculation.
fine fine-structure constant.
EC energy of core.
ioo if ioo > 0, the core energy also includes the orbit–orbit interaction.

MODULE spline_hl

hl(1:ns,1:ks) matrix of the L operator in the B-spline basis (in almost symmetric lower-column mode).
lh the corresponding value of the orbital angular momentum.
vc(1:ns,1:ks) matrix of mass-velocity correction in the B-spline basis.

Internal routines:

Allocate_hl allocates space of for the hl and vc arrays according to parameters from module spline_param.

Format of the bsw-files. The atomic orbitals in the B-spline representation are supposed to be recorded in the unformatted bsw-
files, which are similar to the w-files in the MCHF package [68]. In the bsw-files each orbital is recorded sequentially with the
following information:

1. ebs, z, h, hmax, rmax, ksb, nsb, mbs

2. pbs(1:mbs)

where the notation corresponds to the definitions in the module spline_orbitals. Along with the usual information about the quantum
numbers and the B-spline expansions, the bsw-file also contains the B-spline parameters, which enables us to reconstruct the
corresponding knot sequence.

12.4. Evaluation of two-electron integrals

The BSPLINE library contains three sets of routines for the evaluation of two-electron integrals, which are based on the different
numerical methods. They include the usual direct differential method, the B-spline differential method, and the direct summation
over the B-spline cells. The different methods are suitable for different kind of calculations, have different accuracy and speed, and
different representations of the final results. The direct differential method provides only the total values for the integrals, has the
lowest accuracy but the highest speed. This method can be used in large-scale CI calculations. The B-spline differential method
provides the B-spline representation for different two-electron integrals, that can be used in generating the Hamiltonian matrix,
including the Breit–Pauli terms. The evaluation of the total values of integrals, however, is slower, approximately by one order of
magnitude, than in the direct differential method. In the direct summation over B-spline cells, we first evaluate the two-electron
integrals in the different (r1 , r2 ) cells, defined by the B-splines subintervals [ti , ti+1 ]. The integration is performed by using the
Gauss–Legendre quadratures. This provides the highest accuracy, which is comparable to machine accuracy. This method can be
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 333

used efficiently for both the B-spline representation and the evaluation of total integral values, although the speed is still slower that
in the direct differential method.
Below we outline the relevant theory and the different numerical methods used in BSPLINE to deal with the two-electron
integrals. Then we give a short description of the relevant routines and modules. These routines have been developed on the basis
of algorithms and recommendations, previously presented by Brage and Froese Fischer [49] and by Qui and Froese Fischer [64].
We follow closely their notation.

12.4.1. Evaluation of two-electron integrals by direct differential equation method

The two-electron radial integrals in atomic structure calculations are of the form
∞∞ k
  r<
I k,p
(ρ, σ ; ρ , σ ) = fρ (r1 )fσ (r2 ) f  (r )f  (r ) dr1 dr2 ,
k+p ρ 1 σ 2
(12.12)
r>
0 0

where r< and r> are the smaller and larger of the two coordinates r1 and r2 , respectively. The f (r) functions are linked to the radial
orbitals, p is a positive integer, and k can assume a range of different integer values. If we define a function Y k,p (σ, σ  ; r) so that
∞
Y k,p (σ σ  ; r)
I k (ρ, σ ; ρ  , σ  ) = fρ (r) fρ  (r) dr.
rp
0

It follows that
r ∞
sk r k+p
Y k,p (σ σ  , r) = fσ (s) fσ  (s) ds + fσ (s) fσ  (s) ds. (12.13)
rk s k+p
0 r

It can be shown that this function satisfies a differential equation of the form
 
d k+p 2k + p k
− Y k,p (σ σ  ; r) = − Z (σ σ  ; r),
dr r r
Y k (σ σ  ; r) → Z k (σ σ  ; r) as r → ∞. (12.14)
Here we have introduced a second function
r
 sk
Z (σ σ ; r) =
k
fσ (s) fσ  (s), (12.15)
rk
0

which satisfies the differential equation

 
d k
− Z k (σ σ  ; r) = fσ (r)fσ  (r); Z k (σ σ  ; 0) = 0. (12.16)
dr r
The form of the Z k function is quite general when defined in this way, since it is independent of p.
In Breit–Pauli calculations, all two-electron integrals can be reduced to four main types: R k , N k , V k , and T k integrals.

R k integrals
∞∞ k
r<
R k (ρ, σ ; ρ  , σ  ) = Pρ (r1 )Pσ (r2 ) k+1
Pρ  (r1 )Pσ  (r2 ) dr1 dr2 . (12.17)
r>
0 0

This is the special case of Eq. (12.12) with p = 1, and the f functions being the radial orbitals. The R k integrals are highly
symmetric

R k (ρ, σ ; ρ  , σ  ) = R k (ρ  , σ ; ρ, σ  ) = R k (ρ, σ  ; ρ  , σ ) = R k (σ, ρ; σ  , ρ  ),

i.e. they are unchanged as long as the ρ and ρ  orbitals, and thereby the σ and σ  , are linked together with one coordinate.
The two-step procedure described above can be avoided, since the differential equation for the Y k functions is
 2 
d (k + p)2 − (k + p) k (2k + p)(p − 1) k 2k + p
2
− 2
Y (σ σ  ; r) = − 2
Z (σ σ  ; r) − fσ (r)fσ  (r). (12.18)
dr r r r
334 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

In the special case of p = 1, the coefficient of the Z k function disappears, and we are left with a single, though second-order,
differential equation for Y k
 2 
d k(k + 1) k 2k + 1
− Y (σ σ  ; r) = − fσ (r)fσ  (r) (12.19)
dr 2 r2 r
with the boundary conditions
Yσkσ  (0) = 0,
rmax
Yσkσ  (rmax ) = fσ (r)fσ  (r) dr if k = 0, (12.20)
0
d k k
Y  (rmax ) − Yσkσ  (rmax ) = fσ (rmax )fσ  (rmax ) if k > 0.
dr σ σ r
The boundary condition at rmax can be obtained from (12.13), assuming that the orbitals disappear at r > rmax .

N k integrals
∞∞
α2 r2k
N k (ρ, σ ; ρ  , σ  ) = Pρ (r1 )Pσ (r2 ) k+3
ε(r1 − r2 )Pρ  (r1 )Pσ  (r2 ) dr1 dr2
4 r1
0 0
∞
α2 1 k
= Pρ (r)Pρ  (r) Z (Pσ , Pσ  ; r) dr (12.21)
4 r3
0

where k
−1 and ε(x) = 0 for x  0, ε(x) = 1 for x > 0. The ε(x) step function reduces the Y k function to Z k , and the N k
integral is the special case of Eq. (12.12) with p = 3, Y k = Z k and the f functions as the radial orbitals. Hence, in this case we
need to solve only one first-order differential equation (12.16). But an additional complication occurs when dealing with the case
N −1 , i.e. k < 0. In this case the homogeneous solution to Eq. (12.16) is

Z k (σ σ  , r) = const × r −k
which leads to an ambiguous result when r → 0, since k > 0. This is solved by imposing a stronger boundary condition
d k
Z (σ σ  , 0) = 0 if k < 0. (12.22)
dr
The N k integrals have fewer symmetries than the R k integrals:

N k (ρ, σ ; ρ  , σ  ) = N k (ρ  , σ ; ρ, σ  ) = N k (ρ, σ  ; ρ  , σ ).
In applications, we also use the integrals
1 k
M k (ρ, σ ; ρ  , σ  ) = N̄ k (ρ, σ ; ρ  , σ  ) = N (ρ, σ ; ρ  , σ  ) + N k (σ, ρ; σ  , ρ  ) (12.23)
2
which are as symmetric as the R k integral and differ from the latter only by p = 3.

V k integrals
∞∞ k  
 α2
 r< d 1
V (ρ, σ ; ρ , σ ) =
k
Pρ (r1 )Pσ (r2 ) k+3
− Pρ  (r1 )r2 Pσ  (r2 ) dr1 dr2 . (12.24)
4 r> dr1 r1
0 0
This is the special case of Eq. (12.12) with p = 3 and
 
d 1
fρ (r1 ) = Pρ (r1 ); fρ  (r1 ) = − Pρ  (r1 );
dr1 r1
fσ (r2 ) = Pσ (r2 ); fσ  (r2 ) = r2 Pσ  (r2 );
The V k integral allows only the two-step procedure (12.14)–(12.16), and it has the most restricted symmetry relations:

V k (ρ, σ ; ρ  , σ  ) = V k (ρ, σ  ; ρ  , σ ).
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 335

In applications, we also use more the symmetric integrals [25]:

W k (ρ, σ ; ρ  , σ  ) = V̄ k (ρ, σ ; ρ  , σ  )
= 2V k (ρ, σ ; ρ  , σ  ) − (k + 1)N k+1 (ρ, σ ; ρ  , σ  ) + (k + 2)N k−1 (σ, ρ; σ  , ρ  ), (12.25)
Qk (ρ, σ ; ρ  , σ  ) = V̄¯ k (ρ, σ ; ρ  , σ  ) = V̄ k (ρ, σ ; ρ  , σ  ) − V̄ k (σ, ρ; , σ  ρ  ). (12.26)

T k integrals
∞∞ k    
 α2 r< ∂ 1 ∂ 1
T (ρ, σ ; ρ , σ ) =
k
Pρ (r1 )Pσ (r2 ) + Pρ  (r1 ) + Pσ  (r2 ) dr1 dr2 (12.27)
4(2k + 1) k+1
r> ∂r1 r1 ∂r2 r2
0 0

with p = 1 and
 
d 1
fρ (r1 ) = Pρ (r1 ); fρ  (r1 ) = − Pρ  (r1 );
dr1 r1
 
d 1
fσ (r2 ) = Pσ (r2 ); fσ  (r2 ) = − Pσ  (r2 ).
dr2 r2
Here we can use the same procedure as for the R k integrals, but the symmetry relations are stricter:

T k (ρ, σ ; ρ  , σ  ) = T k (σ, ρ; σ  ; ρ  ).
In the Breit–Pauli approximation, one also frequently uses the U k integrals [25]. In the present implementation we do not use any
special procedure for U k , but instead the relationship

U k (i, j ; i  , j  ) = T k (i, j ; i  , j  ) + T k (i  , j ; i, j  ). (12.28)

List of routines

spline_slater module, contains common arrays used in the evaluation of two-electron integrals by the direct differential
equation method.
bzk computes the spline solution for the Z k function.
facdyk sets up and factorizes the matrix of the Y k operator d2 /dr 2 − k(k + 1)/r 2 in the B-spline basis.
facdzk sets up and factorizes the matrix of the Z k operator d/dr + k/r.
rky evaluates the R k integral by the direct differential method.
nky evaluates the N k integral by the direct differential method.
mky evaluates the M k integral by the direct differential method.
tky evaluates the T k integral by the direct differential method.
vky evaluates the V k integral by the direct differential method.
ykf computes the spline solution for the Y k function.

MODULE spline_slater

kz multipole index for the Z k function.

dzk(1:ktx,ns) factorization of the matrix for the Z k operator d/dr + k/r in the B-spline basis.
ipvtz(1:ns)
ktx = 3*ks-2
ky multipole index for the Y k function.
dyk(1:ktx,1:ns) factorization of the matrix for the Y k operator d2 /dr 2 − k(k + 1)/r 2 in the B-spline basis.
ipvtd(1:ns)
yk(1:ns) B-spline representation of the Y k (or Z k ) function.
336 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

fyk(1:nv,1:ks) values of Y k (or Z k ) function at the Gaussian points.

fc1(1:nv,1:ks), ic1 values at the Gaussian points and pointers for four one-electron
fc2(1:nv,1:ks), ic2 radial functions in the two-electron integral under consideration.
fy1(1:nv,1:ks), iy1
fy2(1:nv,1:ks), iy2
fc(1:nv,1:ks) work array.

Internal routines:

Allocate_slater allocates space for all arrays according to parameters from the module spline_param.
Dealloc_slater deallocates space in the module spline_slater.

Remarks.

• The routines ykf, bzk, ytk, and yvk differ only by their way of calculating the Y k functions. The specific routines ytk and yvk
are incorporated in the calling routines.
• The Y k and T k integrals are defined through the solution of second-order differential equations. The N k integral needs only the
Z k solution of the first-order equation while the V k integral is evaluated by the sequential solution of two first-order differential
equations.
• The ykf routine includes the relativistic shift modification of the Y k function if the parameter rel = ‘TRUE’.
• The quadr routine is used to determine the appropriate boundary conditions (12.20) when k = 0. It is important in the case of
non-orthogonal orbitals.
• These routines are very fast but not as accurate as the cell integration. They can be used for evaluating two-electron integrals,
but not for generating the interaction matrix in the B-spline basis.

Calling sequences

rky nky
--------------- ----------
/ \ / \
ykf yval bzk yval
| |
------------------------- -------------
| | | | \ / | \
facdyk yval vinty quadr dgbsl facdzk yval dgbsl
| |
dgbfa dgbfa

tky vky
---------- ------------
// \ // \
ytk yval yvk yval
| |
-------------------------- -------------------------
| | | | \ | | | | \

facdyk yval vinty quadr dgbsl bzk yval vinty facdzk dgbsl
| | |
dgbfa -------------- dgbfa
facdzk yval dgbsl
|
dgbfa
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 337

12.4.2. Evaluation of two-electron integrals in the B-spline basis

In the B-spline representation of radial atomic orbitals (12.5), the two-electron integrals (12.12) become





I k (a, b; c, d) = ai bi ci  dj  I k (i, j ; i  , j  ), (12.29)
i j i j

where the matrix elements I k (i, j ; i  , j  ) are

∞∞ k
  r<
R (i, j ; i , j ) =
k
Bi (r1 )Bj (r2 ) k+1
Bi  (r1 )Bj  (r2 ) dr1 dr2 ,
r>
0 0
∞ r1
  1 r2k
N (i, j ; i , j ) =
k
dr1 Bi  (r1 )Bi (r1 ) dr2 Bj (r2 ) Bj  (r2 ),
r13 r1k (12.30)
0 0
∞∞ k
r<
T k (i, j ; i  , j  ) = Bi (r1 )Bj (r2 ) k+1
B̄i  (r1 )B̄j  (r2 ) dr1 dr2 ,
r>
0 0
∞∞ k
r<
V k (i, j ; i  , j  ) = Bi (r1 )Bj (r2 ) k+3
B̄i  (r1 )Bj  (r2 )r2 dr1 dr2 ,
r>
0 0

with B̄(r) = B  (r) − (1/r)B(r). The matrix elements (12.30) can be obtained by using the same procedures described above in
Section 12.4.1 after replacing the f (x) functions by the corresponding B-splines.
Therefore, the two-electron integral is given as a four-fold summation, but considerable simplification is possible. Note that

I k (i, j ; i  j  ) = 0 if either |i  − i|
ks or |j  − j |
ks . (12.31)

Furthermore, the integrals may be symmetric relative to the (i, i  ) or (j, j  ) indexes. As a result, the B-spline integrals can be stored
as a matrix with band character relative to the (i, i  ) or (j, j  ) indexes, namely:

R k (ns , ns , ks , ks ),
N k (ns , ns , ks , ks ),
(12.32)
V k (ns , ns , 2ks − 1, ks ),
T k (ns , ns , 2ks − 1, 2ks − 1).

In particular, the array element R k (i, j ; ip , jp ) refers to the matrix element in which i  = i + ip − 1 and j  = j + jp − 1, where
ip , jp = 1, . . . , ks (the symmetric upper-column mode). However, the array element T k (i, j ; ip , jp ) refers to the matrix element in
which i  = i + ip − ks and j  = j + jp − ks , where ip , jp = 1, . . . , 2ks − 1 (the non-symmetric column mode).
Thus, it is convenient to regroup the summation as follows:





 
R (a, b; c, d) =
k
ai ci  bj dj  R (i, j ; i , j ) .
k
(12.33)
i i j j

Let us consider the innermost double summation. Only the matrix elements with j 
j are stored, but the summation in Eq. (12.33)
runs over all values. However, the summation may be written as






bj dj  R k (i, j ; i  , j  ) = bj dj R k (i, j ; i  , j ) + (bj dj  + bj  dj )R k (i, j ; i  , j  ). (12.34)
j j j j j  >j

Let us define a direct “density matrix”


b j dj if jp = 1,
B(j, jp ) = (12.35)
bj dj +jp −1 + dj bj +jp −1 if jp > 1 and j + jp − 1  ns .
With a similar definition for the another pair of indexes,

ai ci if ip = 1,
A(i, ip ) = (12.36)
ai ci+ip −1 + ci ai+ip −1 if ip > 1 and i + ip − 1  ns ,
338 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

the Slater integral is then defined as

ns

ks

ns

ks
R k (a, b; c, d) = A(i, ip ) B(j, jp )R k (i, j ; ip , jp ). (12.37)
i=1 ip =1 j =1 jp =1

Similarly, we can obtain

ns 2ks −1


ns 2ks −1
T (a, b; c, d) =
k
Ā(i, ip ) B̄(j, jp )T k (i, j ; ip , jp ), (12.38)
i=1 ip =1 j =1 jp =1

where
Ā(i, ip ) = ai ci+ip −ks , 1 < i + ip − ks < ns ,
(12.39)
B̄(j, jp ) = bj dj +jp −ks , 1 < j + jp − ks < ns .

The N k integrals are also defined through the ‘symmetric’ summation (12.37), whereas the V k integrals are the mixed case

ns 2ks −1

ns

ks
V k (a, b; c, d) = Ā(i, ip ) B(j, jp )V k (i, j ; ip , jp ). (12.40)
i=1 ip =1 j =1 jp =1

List of routines

spline_integrals module, contains the B-spline representation of the given two-electron integral.
convol convolutes the four-dimensional array of B-spline two-electron spline integrals with a two-dimensional
density matrix.
density provides the density matrices (12.35) or (12.39).
mrk_diff sets up the matrix for the B-spline representation of the R k -integral using the differential equation method.
mnk_diff sets up the matrix for the B-spline representation of the N k -integral using the differential equation method.
mmk_diff sets up the matrix for the B-spline representation of the M k -integral using the differential equation method.
mtk_diff sets up the matrix for the B-spline representation of the T k -integral using the differential equation method.
mvk_diff sets up the matrix for the B-spline representation of the V k -integral using the differential equation method.
rk evaluates the R k integral by assembling the B-spline integrals.
nk evaluates the N k integral by assembling the B-spline integrals.
mk evaluates the M k integral by assembling the B-spline integrals.
tk evaluates the T k integral by assembling the B-spline integrals.
vk evaluates the V k integral by assembling the B-spline integrals.

MODULE spline_integrals

krk multipole index for the integral under consideration.

itype character (rk, nk, mk, tk, or vk) which indicates the type of integral.
rkb the B-spline representation of the given two-electron integral, stored in column mode with dimensions
(1:ns ,1:ns ,1:2ks -1,1:2ks -1).
Internal routines:
Allocate_integrals allocates space for the rkb array according to parameters from module spline_param.
Dealloc_integrals deallocates space in module spline_integral.
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 339

Calling sequences

rk nk

| |

mrk_diff mnk_diff

| |
----------------- -----------------
/ || | \ / || | \

facdyk bspyk yval minty facdzk bspzk yval minty

| | | |
dgbfa dgbsl dgbfa dgbsl

tk vk
| |
mtk_diff mvk_diff
| |
----------------- -----------------
/ || | \ / || | \

facdyk bsptyk yval dinty facdzk bspvyk yval dinty

| | | |
dgbfa dgbsl dgbfa ----------
/ | \
yval vinty dgbsl

Remarks.

• The module spline_integrals contains only one array for all types of integrals. If requirements on memory are not very strict,
it is possible to introduce specific arrays (with different shape) for each type of integrals.
• The rk, nk, mk, vk, and tk routines evaluate the two-electron integrals by assembling the matrix elements in the B-spline
basis with the corresponding orbital expansion coefficients (see Eqs. (12.33)–(12.40)). The calculations of matrix elements in
the B-spline basis are performed by the routines mrk_diff, mnk_diff, mmk_diff, mvk_diff, and mtk_diff by the differential
equation method, or with the routines mrk_cell, mnk_cell, mmk_cell, mvk_cell, and mtk_cell by the cell integration method
(see next section). All these routines first check the parameters krk and itype in the module spline_integrals. Only if these
parameters differ from what is required, the calculation of the corresponding B-spline representation for the given integrals is
initialized.
• The above routines differ mainly by their way of calculating the corresponding Z k and Y k functions between the individual
B-spline functions in the bspyk, bspzk, bsptyk, and bspvyk routines, respectively. These specific routines are incorporated in
the calling routines, and the algorithms are similar to the algorithms in routines ykf, bzk, ytk, and yvk, based upon the direct
integration of two-electron integrals (see Section 12.3.2).
• When the parameter rel = ‘TRUE’, the bspyk routine includes the relativistic shift modification for the Y k function due to
the two-electron Darvin and spin–spin contact terms (see Section 2.6). The fine-structure multiplicative factor α 2 /4 is applied
directly to the rkb array, if it contains a B-spline representation for the relativistic integral.

12.4.3. Cell algorithm for the evaluation of two-electron integrals

In this algorithm, the two dimensional region (r1 , r2 ) is divided into a number of rectangular cells according to the chosen
grid. The direct two-dimensional integration using Gaussian quadrature is performed for the evaluation of integrals within each
cell. Over the off-diagonal cells, the integration is separable, so the two-dimensional cell integral is reduced to a product of two
one-dimensional integrals. Furthermore, the scaling invariance of the B-splines in the logarithmic region of the chosen grid is fully
exploited, and only a restricted amount of cell integrations need to be performed. The values of the required integrals can then
be obtained by assembling the cell integrals. This algorithm significantly improves the efficiency and accuracy of the traditional
methods that rely on the solution of the differential equations.

Integration over the off-diagonal cells. In this case the integrand in the matrix elements is separable and the two-dimensional
integrals are reduced to a product of two one-dimensional integrals. Consider, for example, the Slater integral and assume iv < jv ,
340 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

where iv and jv enumerate the knot subintervals [riv , riv +1 ] and [rjv , rjv +1 ]. Then
rjv +1 riv +1
k
r<
 
R (i, j ; i , j ; iv , jv ) =
k
k+1
Bi (r1 )Bj (r2 )Bi  (r1 )Bj  (r2 ) dr1 dr2
r>
rjv riv
riv +1 rjv +1
1
= r1k Bi (r1 )Bi  (r1 ) dr1 Bj (r2 )Bj  (r2 ) dr2
r2k+1
riv rjv
 −(k+1)
= r (i, i ; iv ) × r
k
(j, j  ; jv ), (12.41)
where
riv +1

r (i, i , iv ) =
m
r m Bi (r)Bi  (r) dr. (12.42)
riv

The r m (i, i  ; iv ) quantities are called the integral moments. Although there are nv (nv − 1) off-diagonal cells, we only need to
calculate and store the nv two-dimensional moment arrays r m (i, i  ; iv ). Recall that each B-spline is a positive polynomial over a
finite range. Hence, the moment integrals can be effectively evaluated using Gaussian quadrature with ks Gaussian points. For a
given order ks , the spline Bi is nonzero only in the range from knot i to knot i + ks , i.e. for a given cell iv , the i and j indexes only
run over the ks values iv , . . . , iv + ks − 1. Hence, the moment arrays r m (i, i  ; iv ) have dimension (1:ks ;1:ks ,1:nv ). Furthermore,
since r m (i, i  ; iv ) is symmetric with respect to interchanging the indexes i, i  , only r m (i, i  ; iv ) with i  i  needs to be evaluated.
The evaluation can be further reduced by using their scaling properties in the logarithmic grid region. Suppose that the splines lie
entirely within the exponential region. Let the left-most knot defining Bi (r) be ti and let r = ti + s. Then
   
Bi (ti + s) = Bi+1 (1 + h)(ti + s) = Bi+1 (ti+1 + s(1 + h)) . (12.43)
Based on this displacement invariance, the following scaling law holds for the moment integrals:
r m (i + 1, i  + 1; iv + 1) = (1 + h)1+m r m (i, i  ; iv ). (12.44)
Expressions similar to (12.41) can also be written for the Nk, V k, and Tk integrals:
N k (i, j ; i  , j  ; iv , jv ) = r −(k+3) (i, i  ; iv ) × r k (j, j  ; jv ), iv > jv ,
T k (i, j ; i  , j  ; iv , jv ) = d −(k+1) (i, i  ; iv ) × d k (j, j  ; jv ), iv > jv ,
T k (i, j ; i  , j  ; iv , jv ) = d k (i, i  ; iv ) × d −(k+1) (j, j  ; jv ), jv > iv , (12.45)
V k (i, j ; i  , j  ; iv , jv ) = d −(k+3) (i, i  ; iv ) × r k+1 (j, j  ; jv ), iv > jv ,
V k (i, j ; i  , j  ; iv , jv ) = d k (i, i  ; iv ) × r −(k+2) (j, j  ; jv ), jv > iv ,
where we introduced the moment integrals d m (i, i  , iv ) containing derivatives:
riv +1  
  1
d (i, i , iv ) =
m
r Bi (r) Bi  (r) − Bi  (r) dr.
m
(12.46)
r
riv

Integration over the diagonal cells. In the diagonal cells, the two coordinates in the integrand of the two-electron integrals are
coupled. Moreover, the integrand is not continuous across the cell diagonal where r1 = r2 . We therefore separate the rectangular
cell into two triangles, so that the integrand inside each triangle is continuous. The integration over the rectangular cell is then a
summation of integrations over the two triangles, which are symmetric with respect to index exchange:
R k (i, j ; i  , j  ; iv ) = R
k
(i, j ; i  , j  ; iv ) + R
k
(j ; i; j  , i  ; iv ), (12.47)
where
riv +1 r1
1
k
R (i, j ; i  , j  ; iv ) = Bi (r1 )Bi  (r1 ) dr1 r2k Bj (r2 )Bj  (r2 ) dr2 . (12.48)
r1k+1
riv riv

We use again Gaussian quadrature for the integration over the triangular cells. Here, however, we should generate new B-splines at
the Gaussian points for the interval [riv , r1 ]. To minimize the number of new B-spline values needed to perform the two-dimensional
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 341

integration, we apply the original Gaussian points used in the evaluation of the one-dimensional integrals to the integration with
respect to the coordinate r1 (i.e., the first integral in (12.48)). In the integration over the coordinate r2 , we also use the ks Gaussian
points for each interval [riv , r1 ], dynamically calculating the required B-spline values. With this two-dimensional grid, the difficulty
of the cell integrals near the origin due to the singularity of the integrand is minimized and uniformly accurate results for all Slater
integrals can be obtained. The full storage of the R k (i, j ; i  , j  ; i ) integrals takes n × k 4 locations ignoring the symmetry over the
v v s
indexes. The computational time for these integrals can also be reduced by using the scaling law in the logarithmic grid region [64]:
k
R (i + 1, j + 1; i  + 1, j  + 1; iv + 1) = (1 + h) × R
k
(i, j ; i  , j  ; iv ). (12.49)
For the Nk integrals we have
N k (i, j ; i  , j  ; iv ) = N
k
(i, j ; i  , j  ; iv ),
where
riv +1 r1
1
k
N (i, j ; i  , j  ; iv ) = Bi (r1 )Bi  (r1 ) dr1 r2k Bj (r2 )Bj  (r2 ) dr2 .
r1k+3
riv riv

Then for the Tk integrals:

T k (i, j ; i  , j  ; iv ) = Tk (i, j ; i  , j  ; iv ) + Tk (j ; i; j  , i  ; iv ),
where
riv +1 r1
1
Tk (i, j ; i  , j  ; iv ) = Bi (r1 )Bi  (r1 ) dr1 r2k B̄j (r2 )B̄j  (r2 ) dr2 ,
r1k+1
riv riv

and, finally, the V k integrals:

V k (i, j ; i  , j  ; iv ) = Vk (i, j ; i  , j  ; iv ),
where
riv +1 r1
1
Vk (i, j ; i  , j  ; iv ) = Bi (r1 )B̄i  (r1 ) dr1 r2k+1 Bj (r2 )Bj  (r2 ) dr2
r1k+3
riv riv
riv +1 r1
1
+ Bj (r1 )B (r1 ) dr1
j r2k Bi (r2 )B̄i  (r2 ) dr2 .
r1k+2
riv riv

For all relativistic integrals, the scaling law in the exponential region is the same, but it differs from that for the R k case:
Ik (i + 1, j + 1; i  + 1, j  + 1; iv + 1) = Ik (i, j ; i  , j  ; iv )/(1 + h). (12.50)
Now one can obtain the B-spline representation for any integral by summing the corresponding cell integrals over the relevant cells.

Direct summation of cell integrals. Consider the cell (iv , jv ). Its contribution to the total R k (a, b; c, d) integral is
iv +k

s −1 jv +k

s −1 iv +k

s −1 jv +k

s −1

R k (a, b; c, d; iv , jv ) = ai bj ci  dj  R k (i, j ; i  , j  ; iv , jv ). (12.51)

i=iv j =jv i  =i v j  =j v

Using the formula for the cell integrals, we obtain

' i +k −1 i +k −1 ('jv +ks −1 jv +ks −1 (
v
s v
s


 −(k+1) 
R (a, b; c, d; iv < jv ) =
k
ai k
ci  r (i, i , iv ) bj dj  r (j, j ; jv )
i=iv i  =iv j =jv j  =jv

= X k (a, c; iv ) × X −(k+1) (b, d; jv ),

' i +k −1 i +k −1 ('jv +ks −1 jv +ks −1 (
v
s v
s


−(k+1)  
R (a, b; c, d; iv > jv ) =
k
ai ci  r (i, i , iv ) bj dj  r (j, j ; jv )
k

i=iv i  =iv j =jv j  =jv

−(k+1)
=X (a, c; iv ) × X (b, d; jv ),
k
342 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

where
iv +k

s −1 iv +k

s −1
X m (x, y; iv ) = xi yi  r m (i, i  , iv ) = x × r m (iv ) × y. (12.52)
i=iv i  =i v

Hence, the off-diagonal cell integrals for the atomic orbitals can also be represented as a product of two quantities, which depend
only on the cell index. Introducing the four vectors

A(iv ) = X k (a, c; iv ),
B(jv ) = X −(k+1) (b, d; jv ),
(12.53)
C(iv ) = X −(k+1) (a, c; iv ),
D(jv ) = X k (b, d; jv ),

the cell representation of the total R k integral has the form

 
 .. 
 DI AG1 C1 D 2 C1 D 3 . C1 Dnv 
 
 .. 
 A2 B1 DI AG2 C2 D 3 . C2 Dnv 
 
R k (a, b; c, d) =  . . (12.54)
 A3 B1 A3 B2 DI AG3 .. C3 Dnv 
 
 ... ... ... ... . . . 
 ..
 
Anv B1 Anv B2 Anv B3 . DI AGnv
Therefore, the R k integral can be obtained by directly summing the above quantities, without generating the corresponding B-spline
representation for the given integral described in the previous sections. Such procedure is much more efficient.

List of routines

spline_moments module, contains the off-diagonals moments and diagonal-cell integrals for a given type of integral.
moments evaluates the moments Bi |r λ |Bj  or Bi |r λ |Bj  for a given power index.
rk_moments generates all off-diagonals moments and diagonal-cell integrals for the R k -integral.
nk_moments generates all off-diagonals moments and diagonal-cell integrals for the N k -integral.
mk_moments generates all off-diagonals moments and diagonal-cell integrals for the M k -integral.
tk_moments generates all off-diagonals moments and diagonal-cell integrals for the T k -integral.
vk_moments generates all off-diagonals moments and diagonal-cell integrals for the V k -integral.
mrk_cell generates the B-spline representation of the R k -integral using the cell-integration method.
mnk_cell generates the B-spline representation of the N k -integral using the cell-integration method.
mmk_cell generates the B-spline representation of the M k -integral using the cell-integration method.
mtk_cell generates the B-spline representation of the T k -integral using the cell-integration method.
mvk_cell generates the B-spline representation of the V k -integral using the cell-integration method.
rkc evaluates the R k integral by direct assembling the moments.
nkc evaluates the N k integral by direct assembling the moments.
mkc evaluates the M k integral by direct assembling the moments.
tkc evaluates the T k integral by direct assembling the moments.
vkc evaluates the V k integral by direct assembling the moments.

MODULE spline_moments. This module contains moments, and the two-electron integral over the diagonal cell, Eq. (12.49).
These moments are used for calculating the two-electron integrals according to the cell algorithm.
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 343

kmk multipole index for the integral under consideration.

mtype character (rk, nk, mk, tk, or vk), which indicates the type of integral.
rkd (1: k s × k s ,1: k s × k s ,1: nv ) stores the two-electron integrals between individual B-splines over the diagonal cells.
rkd1 (1: k s ,1: k s ,1: nv ) one-electron moments r m (i, i  ; iv ) and d m (i, i  ; iv ) as defined in Eqs. (12.44) and (12.46),
rkd2 (1: k s ,1: k s ,1: nv ) which define the integration over the off-diagonal cells for a given type of integrals.
rkd3 (1: k s ,1: k s ,1: nv )
rkd4 (1: k s ,1: k s ,1: nv )

Internal routines:
Allocate_monents allocates space for arrays according to parameters from module spline_param.
Dealloc_moments deallocates space in module spline_moments.

Calling sequences

rk_moments nk_moments
| |
---- ----
/ \\ / \\
moments rk_pdiag moments nk_pdiag
|| ||
rk_triang nk_triang

tk_moments vk_moments
| |
---- ----
/ \\ / \\
moments tk_pdiag moments vk_pdiag
|| ||
tk_triang vk_triang

Remarks.

• All the above routines for the calculations of the moments have the same structure and only differ by their way of calculating
the diagonal cell contributions with the triang-routines. The pdiag-routines control the scaling properties.
• The mrk_cell, mnk_cell, mmk_cell, mtk_cell and mvk_cell routines generate the B-spline representation for the given two-
electron integral and place the corresponding information in the array rbk in the module spline_integral.
• The direct cell-integration routines rkc, nkc, mkc, tkc, and vkc use the direct summation of elements in (12.54). They are the
most effective and accurate way to calculate individual atomic integrals. However, mrk_cell and similar routines should be
used to generate the interaction matrices.

12.5. The library ZCONF

This library contains subroutines and functions that deal with the configuration and close-coupling expansions. They were de-
signed to unify the representation of different expansions and are used in almost all programs of the BSR package. The main
modules in ZCONF are configs, channels and target, with handle data related to the configuration expansion, the scattering chan-
nels, or the description of the target states, respectively. The description of the atomic states in the present package is based upon
the configuration model, and the corresponding expansions for each state should be represented in the so-called c-files, which are
used in the MCHF atomic structure package [52]. For completeness of the description, we recall below the main features of the
c-files.
A general configuration consists of groups of equivalent electrons


m
(n1 l1 )w1 (n2 l2 )w2 . . . (nm lm )wm , N= wi . (12.55)
i=1
344 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

C 3P -37.7603504
1s
2s( 2) 2p( 2) 0.97657139
1S0 3P2 3P
2p( 4) -.13532612
3P2
2s( 1) 2p( 2) 3d( 1) 0.10709232
2S1 1D2 2D1 2D 3P
2s( 1) 2p( 2) 3d( 1) 0.06814467
2S1 3P2 2D1 4P 3P
2s( 1) 2p( 2) 3d( 1) -.06450278
2S1 3P2 2D1 2P 3P
2s( 1) 2p( 1) 3s( 1) 3p( 1) 0.05052374
2S1 2P1 2S1 2P1 3P 2P 3P
2s( 2) 3p( 2) -.03579417
1S0 3P2 3P
2s( 2) 3d( 2) 0.03468992
1S0 3P2 3P
2p( 2) 3d( 2) -.03455937
3P2 1S0 3P
2s( 1) 2p( 1) 3s( 1) 3p( 1) -.02901590
2S1 2P1 2S1 2P1 3P 2P 3P
2p( 2) 3s( 2) -.02546691
3P2 1S0 3P
*
Fig. 14. Structure of a c-file (representing the ground state of C as an example).

Each subshell is additionally characterized by its term αi Li Si . The total orbital angular momentum and the total spin are obtained
by coupling the momenta of all occupied subshells according to the given scheme, completely specified by the intermediate angular
momentum quantum numbers Li , Si . In this package, the coupling is from left to right, for which the notation

γ LS = (n1 l1 )w1 αL1 S1 ; (n2 l2 )w1 αL2 S2 , L2 S2 ; (n3 l3 )w1 αL3 S3 , L3 S3 ; . . . (nm lm )wm αLm Sm , Lm Sm

; LS (12.56)
can be used. Here |L1 − L2 |  L2  |L1 + L2 |, |L3 − L2 |  L3  |L3 + L2 |, etc. Other coupling schemes may prove to be more
advantageous in special circumstances. Therefore, each configuration state function, CSF, is defined by its configuration (12.55)
and the coupling scheme (12.56). The c-file for a given atomic state contains the list of the CSFs involved, with the corresponding
expansion coefficients as an option when the given c-file describes the specific atomic state. An example of a c-file is given in
Fig. 14. All orbital angular momenta in the c-file are given in spectroscopic notation. The configuration states are specified by first
listing the subshells and their occupation in a linear form where, for example, 2p 3 becomes 2p(3), followed by a separate line
with multiplicity (2S + 1), the coupled orbital angular momentum (L), and the seniority, first for each subshell, and then for the
resultant when coupling these subshells from left to right. For the latter, the seniority is no longer relevant. Notice that for the same
configuration there are several possible CSFs, which differ by the intermediate angular quantum numbers. For fine-structure LSJ
atomic states, the c-file may also contain CSFs with different total terms. The expansion coefficients are given at the end of the
configuration line. The first line of the c-file is reserved for the title (optional) and the energy of the given state. In the present
implementation, the user can also define the J -value (if relevant) by including the statement 2J = value at any position in the
first line. The second line in the c-file contains the list of fully occupied, or closed, subshells common to all CSFs. This saves
space and is also justified by the fact that the treatment of closed shells at the computation of matrix elements can be simplified by
redefining the radial integrals for the outer electrons [68].
The spectroscopic notation adopted in the c-files is very convenient for the user, but it imposes several restrictions to the possible
CSFs. First, the spectroscopic symbols for the angular momenta are restricted by the character map used in the computer. For
values up to l = 21 we adopt the usual alphabetical characters s, p, d, f, . . . , z, whereas for higher values (which can be important
in scattering calculations), we use the symbols from the computer character map with indexes from 101 and higher. Note that the
one-electron spectroscopic notation, in addition to the nl-values, may also contain the so-called ‘set index’, which indicates the
non-orthogonal orbitals. For example, we can use the notation 3p1, 3p2, . . . for the term-dependent radial functions for the 3p
orbital in the various 2p5 3p LS states, or we can use different set indexes for orbitals obtained in independent calculations. The c-
file provides only a four-character space for each orbital, and this imposes additional restrictions on the maximum value of possible
principal quantum numbers and set indexes (for details, see the routines al, elf, new_index). However, in practical calculations
these restrictions are usually not critical.
Note also that the c-file format was changed from its original presentation in the MCHF program [52,68]. Previously, it was
assumed that there will be at most five subshells (filled or unfilled) outside the closed subshells common to all CSFs. This turned
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 345

out to be a serious restriction, especially in scattering calculations, where we have an additional continuum orbital, i.e. an additional
shell in the CSF. The format of the c-file was changed, so that we can now use CSFs with up to eight subshells outside the closed
subshells.
In addition to the c-files, the present BSR package can also use the l- and j -files provided by the MCHF complex. The name.l
file represents the list of solutions obtained in the frame of the configuration-interaction method on the basis of configurations in
the file name.c. The name.j file corresponds to the list of solutions obtained in the frame of LSJ Breit–Pauli calculations.
Below we provide only the alphabetical list of modules and routines, with a brief description of the main points. The detailed
description of the input/output parameters can be found at the beginning of the corresponding source codes.

Modules

configs contains different information used for the description of the configuration expansions under
consideration.
channel defines the scattering channels for one partial wave.
channels defines the scattering channels for all partial waves.
target contains the description of the target states.

Routines

allocations provides different array allocations in the modules.

decode_c decodes the spectroscopic representation of the CFSs in the c-file to the ‘integer’ representation used by
the programs.
det_exp provides the determinant expansion for the given subshell.
gen_conf generates all CFSs that can be built from a given configuration by adding the additional electron nl.
idef_jot defines the total electronic angular momentum J for the configuration expansion in the given c-file.
idef_label generates the list of labels for the CSFs in the given c-file.
idef_ncfg defines the number of configurations in the given c-file.
idef_ne defines the number of electrons for the CSFs in the given c-file.
idef_prt defines the parity for the CSFs in the given c-file.
idef_st defines the number of states in the given l- or j -file.
idef_term defines the total term LS for the configuration expansion in the given c-file.
incode_c incodes the ‘integer’ configuration representation used by the programs into spectroscopic c-file format.
new_index assigns a new set index for the orbital on the base of the given orthogonal conditions.
pre_iort prepares the orthogonality conditions for the one-electron radial functions according to their set indexes
and additionally imposed conditions.
pri_conf prints the given configuration in spectroscopic format.
r_closd reads the list of common closed shells from the given c-file.
r_confc reads the configuration list from the given c-file.
r_expn reads the expansion coefficients from the given c-file.
r_orbit defines the list of one-electron orbitals in the given c-file.
r_target reads the information about the target states and the close-coupling expansion from the file target.
r_term defines the list of different terms in the given c-file.
sym_conf defines the configuration and angular symmetries for the given CSF.
test_a checks the “AFTER” relations between the orbitals in the given configuration.
test_c checks the configurations (number of electrons, parity, couplings).
346 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

12.6. The library ZCOM

This library contains subroutines and functions that are part of the BSR package, but they also provide a variety of auxiliary
computations common to different atomic-structure calculations. This includes, for example, the calculation of 3j-symbols or var-
ious sorting procedures. The routines in ZCOM do not use any COMMON or USE Fortran statements, and all input and output
parameters are transferred through the formal arguments only. No call is made from the library to any routine outside the library,
with the exception of interface subroutines used for calling the LAPACK routines.
Below we provide the alphabetical list of routines, with only a brief description of the main points. The input/output parameters
and the method used by each routine are fully commented at the beginning of the corresponding source code. If a user needs to use
some routine in other applications than a BSR calculation, or wants to replace it with a more effective procedure, we advise to first
read the relevant source code for a detailed description of formal arguments.

al provides the spectroscopic symbol (s, p, d, . . .) for an orbital angular momentum l.

clebsh determines a Clebsch–Gordan coefficient via the corresponding 3j-symbol.
conv_au provides the transformation factor from a.u. to cm−1 or eV.
det calculates the determinant value in the Gauss method.
dj_fact defines the J -dependent factor for the reduced matrix elements of an electric multipole transition operator.
djm_fact defines the J -dependent factor for the reduced matrix elements of a magnetic multipole transition operator.
el_nlk decodes the spectroscopic notation for the electron orbital (n, l, k).
elf provides the spectroscopic notation for the electron orbital (n, l, k).
f_gauss gives the value of a given Gaussian at the point x.
gauss looks up the values of the Gaussian coordinates and weights for a k-point Gaussian quadrature over the
interval [0, 1].
gen provides the fractional parentage coefficients for p q - and d q -shells.
hwf routines in this module provide hydrogenic bound radial functions and average values of
r 1 , r 2 , r 3 , r −1 , r −2 , r −3 for a given nuclear charge Z.
iglq total statistical weight of l q -subshell.
interv a reformatted version of the de Boor routine [60] with the same name.
inv calculates the inverse matrix for the square input matrix A(n,n).
ipointer determines the position of an element i in the input integer array, with output value 0 if that element is not
a member of the set.
iterm provides information about possible atomic shell terms.
itri checks the triangular relations for three coupled angular momenta.
la gives the value of the orbital momentum l for a given spectroscopic symbol.
lap_call contains interface routines for calling the most ‘popular’ LAPACK procedures, with the possibility of
automatic allocation of work arrays if needed.
num finds the integer representation for an arbitrary real number r, i.e. it finds such i1 and i2 that abs(r) =
sqrt(i1 /i2 ) and sign(r) = sign(i1 ). It is used for the output of angular coefficients.
read_par a set of routines for reading the parameter values, presented in the input file as name = value.
read_arg a set of routines for reading the input arguments presented in the command line as name = value.
rttc an interface program for timing routines. The user may choose the proper Fortran time routine that is most
appropriate for a given computational platform.
t_ls_jk provides the recoupling coefficients from LS- to jj-coupling.
t_ls_jk provides the recoupling coefficients from LS- to jK-coupling.
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 347

xlagr contains k points of a Lagrange interpolation procedure.

z_3j calculates a 3j-symbol without directly using factorials.
z_6j calculates a 6j-symbol without directly using factorials.
z_9j calculates a 9j-symbol as a sum of 6j-symbols.
zcb provides the square of a 3j-symbol with zero magnetic numbers, using a simplified algorithm.
zclkl calculates the reduced matrix elements of spherical harmonics, l C k l  , without directly using factorials.

13. Additional programs-utilities

In addition to the programs that were described separately, a number of short utility programs are available that assist in the
processing and managing of the data.

BOUND_BSW

Description: converts the close-coupling B-spline expansions to the c- and bsw-files

Input files: [Link], bound_bsw.inp
Output files: c- and bsw-files for given states
Call: bound_bsw < bound_bsw.inp

The results of the bound-state B-spline calculations are recorded in the [Link] files in the form of the corresponding close-
coupling expansions (7.3). They are similar to the l- or j -files in the MCHF complex, and contain descriptions of all states for a
given total term. The BOUND_BSW utility collects the information for a given bound state indicated in the file bound_bsw.inp,
and records it as pair name.c and [Link]. These files can be used as input target files for the BSR scattering calculations.

BOUND_TAB

Description: produces the total list of resulting bound states

Input files: [Link]
Output files: [Link]
Call: bound_tab

This program sorts the energy levels in different [Link] files and prints them in atomic units, eV or cm−1 , relative to the lowest
state. The user can define the range of partial waves and the energy range to restrict the output.

BSW_W

Description: converts the bsw-file to the w-file

Input files: [Link]
Output files: name.w
Call: bsw_w [Link]

CFILE

Description: extracts the given state from the l- or j -file in separate c-files
Input files: name.l or name.j
Output files: result.c
Call: cfile name.l(j) nn 2J result.c eps_c

In the BSR complex, each target state is described by a pair name.c and [Link]. The CI program from the MCHF complex
produces the results in the l(j )-files as a list of states. The CFILE utility allows one to extract the expansion coefficients for a
348 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

given state and create the corresponding c-file. The user should indicate in the command line the name of the input l- or j -file, the
pointer for the given state, nn, the value of 2J if it is a j -file, the name for the output c-file, and the tolerance for the expansion
coefficients, eps_c.

PHOTO_TAB

Description: produces the tables for channel or total cross sections and asymmetry parameters
Input files: bsr_phot.nnn
Output files: indicated by user
Call: photo_tab with interactive response

This utility serves for processing the data recorded during photoionization calculations with the program BSR_PHOT. The program
BSR_PHOT can be run for different energy intervals and different partial waves separately. The results are accumulated in files
bsr_phot.nnn. Then final tables for channel or total cross sections and asymmetry parameters may be generated with the utility
photo_tab.

SLATER_W

Description: transfers Slater-type orbitals to MCHF format

Input files: [Link]
Output files: name.w
Call: slater_w [Link]

The Slater orbitals are supposed to be given as in the RMATRX1 code, i.e.


NCO
Pnl (r) = C(j ) ∗ r IRAD(j ) ∗ e−ZE(j )∗r .
j =1

The formatted file [Link] contains the Slater-type orbital parameters in the form:

1. Z, atom, term nuclear charge, and optional identification for atom and term
2. nwf number of radial wave functions
3. n,l,k quantum numbers for given orbital
4. NCO
5. IRAD(1:NCO)
6. ZE(1:NCO)
7. C(1:NCO)

Repeat records 3–7 nwf times.

SUM_HH

Description: merges the set of [Link] files to the final [Link]

Input files: [Link]
Output files: [Link]
Call: sum_hh klsp1 klsp2

W_BSW

Description: converts the w-file to the bsw-file format

Input files: name.w
Output files: [Link]
Call: bsw_w name.w
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 349

14. Compilation and test runs

This section describes the main features concerning the distribution, installation and compilation, and running of the system.
Several test runs are presented in order to verify that the system works as expected.

14.1. Compilation

The description below and the script files are prepared for a UNIX/LINUX user in the c-shell. Other user can take these files as
a guide for compilation in a specific environment.
The source-code files for the BSR complex are distributed in the compressed bsr_cpc.[Link] file. In order to uncompress the
source files, the user may employ the system routines

gzip -d bsr_cpc.[Link]
tar -xvf bsr_cpc.tar.

All source files are unpacked into the directory BSR_CPC. Each main program, libraries or utilities are presented by a single
source file. The instructions for compiling the entire system are given in the make_BSR script file. Before running make_BSR, the
user should check the assignment of two additional environmental variables in this file. First the user should define the name of the
FORTRAN compiler and optimization level in the variable FC. For example:

setenv FC ’fort -fast’.

Then the user should indicate in the variable LAPACK the directory where the linker can find the LAPACK and BLAS libraries.
Recall that these popular libraries were chosen to provide the linear-algebra calculations in the BSR complex, and the user should
install the LAPACK and BLAS libraries before compiling the BSR complex. The LAPACK and BLAS libraries can be obtained
without charge from the web-site [Link] Many FORTRAN compilers also provide these libraries
as an additional option. If the user copies the liblapack.a and libblas.a files to the BSR_CPC directory, then the variable
LAPACK is defined as

setenv LAPACK ’-L. -llapack -lblas’.

Otherwise the corresponding path should be indicated after the parameter -L.
After setting the environment variables, the make_BSR file provides the recommended sequence of commands for compilation
of the libraries and other BSR programs. The user can run this file as a command script file. As a result, all executable files for
the main programs and utilities will be placed in the directory BSR_CPC. In order to use these programs in the test runs or in
the production calculations, the user should make the corresponding modification in the shell variable path. For example, as a
command

set path = ($HOME/BSR_CPC $path)

or by putting

setenv PATH ”$HOME/BSR_CPC: $PATH”

into the shell script file.

14.2. Test runs

We provide four test runs to illustrate the main modes of operation of the BSR complex, namely, LS-coupling scattering calcula-
tions, JK-coupling scattering calculations, photoionization, and bound-state calculations. In order to provide examples of realistic
calculations and choices of appropriate physical models, the test runs were taken from published calculations for bound states [69]
and electron scattering [70] from carbon atoms. Since these calculations were rather extensive, the test runs are given here in a
simplified version. They only serve as a benchmark to ensure that the programs are running correctly, and as a template for users to
construct their own runs.
Each scattering calculation begins with the choice and generation of the target states. In the framework of the BSR complex, this
is the responsibility of the user. In principle, the user can employ any atomic-structure program, and then convert the results to the
form acceptable in the BSR programs. The BSR programs themselves are adapted to use the target states from the MCHF complex
of Froese Fischer [68]. Below we assume that all target states are generated in the framework of the MCHF approach, and we will
show the run sequences in the different modes.
350 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

In the following test-run description we only discuss the input data files and the sequence of required commands. A complete
set of input files and files with the final results is provided in the separate tar-file bsr_test.[Link], which contains the main directory
BSR_TEST and four subdirectories C_LS, C_JK, C_PHOT, C_BOUND for each test run described below. Each subdirectory
contains the file run_script, with a short description of the test and a required sequence of commands. Note that each target state is
represented by a pair of c- and frm-files. The radial functions are given in the ASCII frm-file instead of unformatted w- or bsw-files,
in order to simplify the test runs on different computer platforms, which may use a different default structure for unformatted files.
As a consequence, the run_script file contains first a rather long sequence of preliminary commands for transformation of ASCII
frm-files to the unformatted bsw-files required by the BSR complex. Each test-run subdirectory contains the file with the correct
final results. These files have the extension .test. At the end of each run_script file we give information what files with final results
should be compared in order to verify that the system works correctly.

14.3. LS-coupling test run

This test corresponds to low-energy elastic and inelastic electron scattering from carbon atoms. The close-coupling expansion
includes the seven lowest states of neutral carbon with configurations 1s2 2s2 2p2 , 1s2 2s2p3 , and 1s2 2s2 2p3s. These states exhibit a
strong term-dependence of the one-electron orbitals, and it is difficult to obtain accurate wavefunctions with only one orthogonal
set of orbitals. In order to illustrate the advantage of the present complex, the atomic wavefunctions for odd-parity and even-
parity states were generated using two different non-orthogonal sets of orbitals. This allows us to obtain fairly accurate target
wavefunctions with relatively small configuration expansions. The results of the MCHF calculations are presented in the c- and
w-files, for the configuration expansions and the one-electron radial functions, respectively. These files are the input for the present
BSR calculations. The c-file for the carbon ground state is given in Section 12.5.
The first step of the BSR calculations concerns the choice of the B-spline knot sequence. The user should prepare the [Link]
file, indicating the R-matrix radius and the maximum size of subintervals, as discussed in Section 12.2 (see also Fig. 13). The
w-files are then transformed to bsw-files with the program-utility w_bsr (see Section 13). At this stage the user can check the
R-matrix radius by considering the value of the one-electron orbitals at the boundary. It is recommended that all orbitals should
at least have fallen to magnitudes < 10−3 . Recall that the value of the R-matrix radius also automatically defines the number of
B-splines. Since the time of the calculations is approximately proportional to the square of the number of B-splines, the user should
choose the minimal possible R-matrix radius with a minimum number of B-splines.
The next step involves the preparation of the file target, with the specification of the scattering model, as shown below.

e + C
--------------------------------------------------------------------------
coupling = LS ! LS coupling
nz = 6 ! nuclear charge
nelc = 6 ! number of electrons
--------------------------------------------------------------------------
ntarg = 7 ! number of target states
--------------------------------------------------------------------------
2p2_3P.c
2p2_1D.c
2p2_1S.c
2p3_5So.c
2p3s_3Po.c
2p3s_1Po.c
2p3_3Do.c
--------------------------------------------------------------------------
nlsp = 12 ! number of partial waves
--------------------------------------------------------------------------
2Se 0 2 1 no
4Se 0 4 1 no
2So 0 2 -1 no
4So 0 4 -1 p_4So.c
2Pe 1 2 1 no
4Pe 1 4 1 no
2Po 1 2 -1 p_2Po.c
4Po 1 4 -1 no
2De 2 2 1 no
4De 2 4 1 no
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 351

2Do 2 2 -1 p_2Do.c
4Do 2 4 -1 no
--------------------------------------------------------------------------

The test run contains 12 partial waves, and for three partial waves there are additional c-files. These provide an example of
(N + 1)-electron configurations, which may be added to the close-coupling expansion to improve the description of short-range
correlations. In the present example, the (N + 1)-electron configurations were added for a more accurate description of the lowest
resonance states with configuration 2p3 . The orbitals in these files were optimized in separate MCHF calculations for each given
term of the 2p3 configuration, and they are not orthogonal to the target orbitals. These additional multi-configuration expansions for
the 2p3 states may be treated in two ways. First, they may be completely included in the close-coupling expansion, in which case
we need to impose the orthogonality constraint for the corresponding continuum orbital according to kp|3p = 0. Alternatively, we
may remove from these perturber files the main configuration 2p3 , leaving only the correlation configurations. We chose this way
in the test run.
The last step in preparing the input data concerns the bsr_par file. Its parameters will control the specific run. This test run has
only a few parameters that differ from the default values, namely:

# bsr_conf

< 2s| ks>=0

# bsr_hd:

itype = 0
idiag = 2
Emax = -0.1

The above orthogonality condition requires that all ks continuum orbitals should be orthogonal to the 2s orbital. By default,
the continuum orbitals are not orthogonal to the bound target orbitals, but it is very important to keep orthogonality to the closed
subshell as discussed in Section 5.1. Orthogonality to the closed core shells (here the 1s orbital) is imposed automatically. The
orthogonality of the ks orbital to the 2s orbitals in the channel 2s2 2p2 ks is also imposed automatically by the BSR_CONF program.
However, it is also important to keep the ks orbital orthogonal in all other channels (such as 2s2p3 ks); otherwise different channels
may generate the same (N + 1)-electron states (2s2 2p3 states in the present example). Note that no additional (N + 1)-electron
states are used in this case to compensate for the orthogonality conditions.
In this test we also chose the two-step diagonalization procedure (idiag=2, see Section 8.1). It is more effective, because it
reduces the size of the interaction matrix to be diagonalized.
After preparing the target and input files, the run of the BSR programs is described by the following script file:

bsr_prep
bsr_conf
bsr_breit klsp1=1 klsp2=12
bsr_mat klsp1=1 klsp2=12
bsr_hd klsp1=1 klsp2=12
sum_hh 1 12

The result of this run will be the [Link] file, which may be used in subsequent cross-section calculations with some external-region
program. If additional partial waves are needed (or we want to recalculate some partial waves), we may run additional calculations
only for specific partial waves and then again the sum_hh utility to obtain the combined [Link] file. The user may control the
calculations using the information from the BSR_MAT and BSR_HD programs in the corresponding log-files. For example, the
BSR_MAT program compares the input and the calculated target energies. If they do not agree within some tolerance, it means that
the input c- or w-files contain errors. If the BSR_HD program fails to diagonalize the interaction matrix, this means that the total
overlap matrix is not positively definite. This may indicate that not all important orthogonal constraints are imposed, or that the
input target states are not orthogonal or not normalized. In order to verify that the system works correctly, user should only compare
the final R-matrix eigenvalues given by the system in the file sum_hh.log and provided in file sum_hh.test for comparison. If these
final eigenvalues do not agree with each other, we suggest to check all log-files for possible warnings.

14.4. jK-coupling test run

This test corresponds to low-energy elastic and inelastic scattering of electron by atomic carbon in the same approximation as
the above LS-coupling test run, but with the inclusion of relativistic corrections. The target file now has the following structure:
352 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

e + C
-----------------------------------------------------------------------
coupling = JK ! coupling
nz = 6 ! nuclear charge
nelc = 6 ! number of electrons
-----------------------------------------------------------------------
ntarg = 13 ! number of target states
-----------------------------------------------------------------------
2p2_1D2.c
2p2_1S0.c
2p2_3P0.c
2p2_3P1.c
2p2_3P2.c
2p3s_1Po1.c
2p3s_3Po0.c
2p3s_3Po1.c
2p3s_3Po2.c
2p3_3Do1.c
2p3_3Do2.c
2p3_3Do3.c
2p3_5So2.c
-----------------------------------------------------------------------
nlsp = 12 ! number of partial waves
-----------------------------------------------------------------------
1+ 1 0 1 no
1- 1 0 -1 p_1o.c
3+ 3 0 1 no
3- 3 0 -1 p_3o.c
5+ 5 0 1 no
5- 5 0 -1 p_5o.c
7+ 7 0 1 no
7- 7 0 -1 no
9+ 9 0 1 no
9- 9 0 -1 no
11+ 11 0 1 no
11- 11 0 -1 no
-----------------------------------------------------------------------
The main difference from the LS-case is that target c-files are now provided for each J -level. In the MCHF package, the Breit–
Pauli corrections are included in the framework of the CI calculations, and the solutions are recorded in the j -files, which contain
the configuration expansions for all J -values. In order to obtain the c-file for a specific J -level, the user may use the utility cfile
(see Section 13). For example, the command

cfile even.j 1 4 2p2_3P2.c 0.0000001

will generate the c-file for the state 2s2 2p2 3 P2 , provided that CI results are recorded in the file even.j. The bsr_par file contains
only one additional parameter, mrel = 2, indicating that the Hamiltonian now includes also all one-electron Breit–Pauli terms. The
sequence of programs in the script file is the same as in the LS test run, with one correction for the BSR_BREIT program, namely:

bsr_breit
oper=1111000 klsp1=1 klsp2=12

The additional parameter oper indicates that the calculations of the angular coefficients also include the spin–orbit interaction.
Recall that the angular coefficients are recorded in the databank files int_bnk.nnn, which can be used in repeated or corrected
calculations.

14.5. Photoionization calculations

As an example of a photoionization calculation, we will consider the photodetachment of the negative carbon ion from its ground
state 2p3 4 So . The final term can only be 4 P, and hence the target file only contains expansions for two partial waves, corresponding
to the initial bound and the final continuum channels:
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 353

hv + C-
--------------------------------------------------------------------
coupling = LS ! LS coupling
nz = 6 ! nuclear charge
nelc = 6 ! number of electrons
--------------------------------------------------------------------
ntarg = 7 ! number of target states
--------------------------------------------------------------------
2p2_3P.c
2p2_1D.c
2p2_1S.c
2p3_5So.c
2p3s_3Po.c
2p3s_1Po.c
2p3_3Do.c
--------------------------------------------------------------------
nlsp = 2 ! number of partial waves
--------------------------------------------------------------------
4So 0 4 -1 p_4So.c
4Pe 1 4 1 no
--------------------------------------------------------------------
Note that the description of the initial state through the close-coupling expansion is optional, because the BSR complex also
allows to represent the initial state as a real bound state obtained separately from independent MCHF calculations. A typical run of
the BSR program for the case of photoionization calculations is described by the following script file:

bsr_prep
bsr_conf
bsr_breit klsp1=1
bsr_mat klsp1=1
bsr_hd klsp1=1 itype=-1 msol=1
bsr_breit klsp1=2
bsr_mat klsp1=2
bsr_hd klsp1=2 itype=1 idiag=1
mult cfg.001 cfg.002 E1
bsr_dmat cfg.001 cfg.002 b p
bsr_phot

We generate the bound state in the first partial wave, and the results are stored in the file bound.001. Then we generate the h.002 and
rsol.002 files in the second partial wave, with the latter file containing all R-matrix solutions in the internal region. In comparison
to the scattering calculations described in Section 14.1 above, we have two additional runs in order to generate the dipole matrix.
The MULT program generates the databank for the dipole angular coefficients, and the BSR_DMAT program generates the dipole
matrix file d.002 itself. The h.002 and d.002 files are then used by the BSR_PHOT program for the photoionization calculations.
The parameters and input energies are provided in the file bsr_phot.inp (see Section 11). The user can run the BSR_PHOT program
several times for different energy intervals and with different energy steps. The results are accumulated in the files photo.002 and
bsr_photo.002. The file photo.002 contains the total photoionization cross sections and can be used for a quick check of the
calculations. The more detailed data recorded in the file bsr_photo.002 may then be used for the calculation of channel cross
sections or asymmetry parameters (see the description of the program utility photo_tab in Section 13). If we have several final
continuum states with different total terms, the calculations are supposed to be carried out for each partial wave separately, and the
total cross sections or asymmetry parameters may be generated again with the utility photo_tab. In order to verify that the system
works correctly, the user should only compare the final photoionization cross sections in the file photo.002 with the data provided
in the file [Link].

14.6. Bound-state calculations

This test is a simplified version of the calculation of oscillator strengths in carbon, presented in detail in [69]. The task was to
obtain the oscillator strengths for transitions from the ground state to the Rydberg states 2s2 2pns and 2s2 2pnd for intermediate
354 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

values of n = 4–8. In this range traditional methods, such as MCHF or the quantum defect method, are difficult to apply. The
Rydberg series 2s2 2pns and 2s2 2pnd strongly interact with each other through both electrostatic and spin–orbit interaction. Besides,
there are also strong corrections from core-valence and inner-core correlations. The test presents the simplest model of describing
the 2s2 2pns and 2s2 2pn series though the close-coupling expansion with five target states: 2s2 2p 2 Po , 2s2p2 2 S, 2 P, 2 D, and 2p3
2 Po . The 2s2p2 states are included to simulate the core-valence correlation due to the strong 2s-2p transition, whereas the inclusion

of the 2p3 state accounts, to some extent, for the inner-core correlation. The target file for this case has the following form:

C - bound
------------------------------------------------------------------------
coupling = LS ! LS coupling
nz = 6 ! nuclear charge
nelc = 5 ! number of electrons
------------------------------------------------------------------------
ntarg = 5 ! number of target states
------------------------------------------------------------------------
2s2_2p_2Po.c
2s_2p2_2D.c
2s_2p2_2P.c
2s_2p2_2S.c
2p3_2Po.c
------------------------------------------------------------------------
nlsp = 4 ! number of partial waves
------------------------------------------------------------------------
0- 0 0 -1 p0.c
1- 2 0 -1 p1.c
2- 4 0 -1 p2.c
3- 6 0 -1 p3.c
------------------------------------------------------------------------
All target states in this example are represented by simple single-configuration wavefunctions, but with term-dependent 2s and
2p one-electron orbitals. We are interested in the final J -level splitting, so the partial waves are constructed for total J -values from
0 to 3. Note that the target states here are represented in LS coupling. Hence this test is an example of using the LSJ coupling in the
framework of the present BSR complex. The bsr_par input file for this test has the form:
# bsr_conf
< 2s| ks>=0
< 2p| kp>=0
max_LT = 7
max_ST = 3
max_ll = 3
# bsr_mat:
mrel = 2
# bsr_hd:
itype = -1
jmvc = 2
idiag = 2
Emax = -0.1
msol = 30
n_max= 12
qd_0 = 1.5
qd_1 = 0.6
qd_2 = 0.2
qd_3 = 0.1

The most important difference from the previous examples is the orthogonality condition imposed on both the 2s and 2p bound
orbitals. It is due to the fact that in this case the same 2s2p3 bound states can be generated in different channels. This, in turn,
may lead to an ‘overloading’ of the total overlap matrix as discussed in Section 5. The corresponding 2s2p3 states are included in
the close-coupling expansion through the additional perturber expansion p0.c, p1.c, etc. For a more accurate representation of the
2s2p3 states, their wavefunctions were also obtained with a full re-optimization of the 2s and 2p orbitals. The mrel parameter in
the above list indicates that the interaction matrix will include all one-electron Breit–Pauli terms, and itype = −1 means that it is a
bound-state calculation with zero boundary condition. In order to obtain the states with high principal quantum numbers we use a
border radius of 300 a.u. However, this does not lead to an increase in the number of B-splines, because we can use the exponential
 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356 355

radial grid (see Section 12.2) for a bound-state calculation. The input parameters also contain the values for the expected quantum
defects in the 2s2 2pnl series. These parameters have no effect on the actual calculations; they are used only for a more appropriate
labeling of the resulting states.
The run of the BSR programs in this test is described by the following script file:

bsr_prep
bsr_conf
bsr_breit klsp1=1 klsp2=4 oper=1111000
bsr_mat klsp1=1 klsp2=4
bsr_hd klsp1=1 klsp2=4
bound_tab
mult 2s2_2p2_3P.c cfg.001 E1
mult 2s2_2p2_3P.c cfg.002 E1
mult 2s2_2p2_3P.c cfg.003 E1
mult 2s2_2p2_3P.c cfg.004 E1
bsr_dmat 2s2_2p2_3P0.c cfg.002 c b
bsr_dmat 2s2_2p2_3P1.c cfg.001 c b
bsr_dmat 2s2_2p2_3P1.c cfg.002 c b
bsr_dmat 2s2_2p2_3P1.c cfg.003 c b
bsr_dmat 2s2_2p2_3P2.c cfg.002 c b
bsr_dmat 2s2_2p2_3P2.c cfg.003 c b
bsr_dmat 2s2_2p2_3P2.c cfg.004 c b

After running the BSR_HD program we obtained the set of [Link] files with solutions for a given J -value. The final list of
bound states in the [Link] files may be obtained with the program utility bound_tab. The initial states 2s2 2p2 3 P0,1,2 in this
test are supposed to have been obtained in separate configuration-interaction calculations, with a list of configurations being placed
in the file 2s2_2p2_3P.c and the J -dependent expansion coefficients in the files 2s2_2p2_3Pj.c. First, we run the MULT program
to obtain all the required dipole angular coefficients in the databank mult_bnk. Then we run the program BSR_DMAT to calculate
the oscillator strengths for each combination of initial and final J -values. Finally, the oscillator strengths from different runs are
accumulated in the file zf_res. The file zf_res.test contains the correct results provided for comparison.

Acknowledgements

I am very grateful to Prof. Charlotte Froese Fischer, who proposed and initiated this work and also provided her B-spline library
and suggested many useful improvements to the code. I am grateful to Dr. S. Tayal and Prof. T. Gorczyca for assisting with the
first implementation of the code in scattering calculations and for helpful discussions at different stages of this work. I would like
to thank Prof. K. Bartschat for many useful discussions, extensive testing of the code, and careful reading of the manuscript. This
work was supported by the United States National Science Foundation under grant PHY-0311161.

References

[1] P.G. Burke, A. Hibbert, W.D. Robb, J. Phys. B 4 (1971) 153.

[2] P.G. Burke, W.D. Robb, Adv. At. Mol. Phys. 11 (1975) 143.
[3] P.G. Burke, K.A. Berrington, Atomic and Molecular Processes: An R-Matrix Approach, IOP Publishing, Bristol, 1993.
[4] K.A. Berrington, W.B. Eissner, P.H. Norrington, Comput. Phys. Comm. 92 (1995) 290.
[5] K.A. Berrington, P.G. Burke, J.J. Chang, A.T. Chivers, W.D. Robb, K.T. Taylor, Comput. Phys. Comm. 8 (1974) 149.
[6] K.A. Berrington, P.G. Burke, M. LeDourneuf, W.D. Robb, K.T. Taylor, Vo Ky Lan, Comput. Phys. Comm. 14 (1978) 367.
[7] N.S. Scott, K.T. Taylor, Comput. Phys. Comm. 25 (1982) 347.
[8] A.G. Sunderland, C.J. Noble, V.M. Burke, P.G. Burke, Comput. Phys. Comm. 143 (2002) 311.
[9] K.G. Dyall, I.P. Grant, C.T. Johnson, F.A. Parpia, E.P. Plummer, Comput. Phys. Comm. 55 (1989) 425.
[10] P.H. Norrington, I.P. Grant, J. Phys. B 20 (1987) 4869; see also [Link]
[11] K. Bartschat, Comput. Phys. Comm. 75 (1993) 219.
[12] P.G. Burke, M.P. Scott, Electron and Positron Scattering from Atoms and Ions, in: K. Bartschat (Ed.), Computational Atomic Physics, Springer, Heidelberg,
1996.
[13] P.G. Burke, C.J. Noble, M.P. Scott, Proc. Roy. Soc. A 410 (1987) 911.
[14] U. Fano, C.M. Lee, Phys. Rev. Lett. 31 (1973) 1573.
[15] M. Plummer, C.J. Noble, J. Phys. B 32 (1999) L345.
[16] P.J.A. Buttle, Phys. Rev. 160 (1967) 719.
[17] K. Bartschat, E.T. Hudson, M.P. Scott, P.G. Burke, V.M. Burke, J. Phys. B 29 (1996) 115.
356 O. Zatsarinny / Computer Physics Communications 174 (2006) 273–356

[18] T.W. Gorczyca, N.R. Badnell, J. Phys. B 30 (1997) 3897.

[19] T.W. Gorczyca, F. Robicheaux, M.S. Pindzola, D.C. Grifin, N.R. Badnell, Phys. Rev. A 52 (1995) 3877.
[20] K. Bartschat, Comput. Phys. Comm. 114 (1998) 168.
[21] A. Hibbert, Comput. Phys. Comm. 1 (1970) 359.
[22] A. Hibbert, C. Froese Fischer, Comput. Phys. Comm. 64 (1991) 417.
[23] A. Hibbert, R. Glass, C. Froese Fischer, Comput. Phys. Comm. 64 (1991) 455.
[24] O. Zatsarinny, Comput. Phys. Comm. 98 (1996) 235.
[25] O. Zatsarinny, C. Froese Fischer, Comput. Phys. Comm. 124 (2000) 247.
[26] H. Bachau, E. Cormier, P. Decleva, J.E. Hansen, F. Martin, Rep. Prog. Phys. 64 (2001) 1815.
[27] H.W. van der Hart, J. Phys. B 30 (1997) 453.
[28] V.M. Burke, C.J. Noble, Comput. Phys. Comm. 85 (1995) 471.
[29] O. Zatsarinny, C. Froese Fischer, J. Phys. B 33 (2000) 313.
[30] O. Zatsarinny, S.S. Tayal, J. Phys. B 34 (2001) 3383.
[31] O. Zatsarinny, S.S. Tayal, J. Phys. B 35 (2002) 241.
[32] O. Zatsarinny, T.W. Gorczyca, C. Froese Fischer, J. Phys. B 35 (2002) 4161.
[33] O. Zatsarinny, T.W. Gorczyca, Abstracts of XIII ICPEAC, We026, Th006, 2003.
[34] O. Zatsarinny, K. Bartschat, J. Phys. B 37 (2004) 2173.
[35] O. Zatsarinny, K. Bartschat, J. Phys. B 37 (2004) 4693.
[36] M.A. Crees, M.J. Seaton, P.M.H. Wilson, Comput. Phys. Comm. 15 (1978) 23.
[37] R.J.W. Henry, S.P. Rountree, E.R. Smith, Comput. Phys. Comm. 23 (1981) 233.
[38] I. Bray, A.T. Stelbovics, Comput. Phys. Comm. 85 (1995) 1.
[39] C. Froese Fischer, M. Idrees, Comput. Phys. 3 (1989) 53.
[40] T. Brage, C. Froese Fischer, G. Miecznik, J. Phys. B 25 (1992) 5289.
[41] M. Venuti, P. Decleva, J. Phys. B 30 (1997) 4839.
[42] J. Xi, C. Froese Fischer, Phys. Rev. A 59 (1999) 307.
[43] M.A. Crees, Comput. Phys. Comm. 19 (1980) 103.
[44] M.J. Seaton, J. Phys. B 5 (1972) L91.
[45] H.E. Saraph, J. Phys. B 13 (1980) 3129.
[46] C. Mendoza, J. Phys. B 14 (1981) 397.
[47] M.J. Seaton, J. Phys. B 18 (1985) 2111.
[48] K.A. Berrington, M.J. Seaton, J. Phys. B 18 (1985) 2587.
[49] T. Brage, C. Froese Fischer, J. Phys. B 27 (1994) 5467.
[50] T. Brage, C. Froese Fischer, Phys. Scripta 49 (1994) 651.
[51] A.I. Akhiezer, V.B. Berestetsky, Quantum Electrodynamics, Interscience, New York, 1960.
[52] C. Froese Fischer, T. Brage, P. Jönsson, Computational Atomic Structure. An MCHF Approach, IOP Publishing, Bristol, 1997.
[53] A. Hibbert, Comput. Phys. Comm. 9 (1975) 141.
[54] W.B. Eissner, M. Jones, H. Nussbaumer, Comput. Phys. Comm. 8 (1974) 270.
[55] M. Bentley, J. Phys. B 27 (1994) 637.
[56] C. Froese Fischer, M.R. Godefroid, A. Hibbert, Comput. Phys. Comm. 64 (1991) 486.
[57] P.O. Löwdin, Phys. Rev. 97 (1955) 1474.
[58] D.A. Varshalovich, A.N. Moskalev, V.K. Khersonskij, Quantum Theory of Angular Momentum, World Scientific, Singapore, 1988.
[59] N. R Badnell, J. Phys. B 32 (1999) 5583; see also [Link]
[60] C. de Boor, A Practical Guide to Splines, Springer, New York, 1978.
[61] C. Froese Fischer, Int. J. Supercomp. Appl. 5 (1991) 5.
[62] C. Froese Fischer, in: T.N. Chang (Ed.), Many-Body Theory of Atomic Structure, Photoionization, World Science, Singapore, 1992.
[63] C. Froese Fischer, W. Guo, Z. Shen, Int. J. Quant. Chem. 42 (1992) 849.
[64] Y. Qiu, C. Froese Fischer, J. Comput. Phys. 156 (1999) 257.
[65] J. Sapirstein, W.R. Johnson, J. Phys. B 29 (1996) 5213.
[66] J.E. Hansen, M. Bently, H.W. van der Hart, M. Landtman, G.M.S. Lister, Y.-T. Shen, N. Vaeck, Phys. Scripta T 47 (1993) 7.
[67] T.N. Chang, T.K. Fang, Phys. Rev. A 52 (1995) 2638.
[68] C. Froese Fischer, Comput. Phys. Comm. 64 (1991) 369.
[69] O. Zatsarinny, C. Froese Fischer, J. Phys. B 35 (2002) 4669.
[70] O. Zatsarinny, K. Bartschat, L. Bandurina, V. Gedeon, Phys. Rev. A 71 (2005) 042702.