GRASPC - GRASP package adapted for the generation of continuum orbitals wave functions

This repository, GRASPC, is a fork of The General-purpose Relativistic Atomic Structure Package (GRASP) - a set of Fortran 90 programs for performing fully relativistic electron structure calculations of atoms.

Link to GRASP repository

About GRASPC

GRASPC is a fork of GRASP, intended to incorporate the scattered (continuum) electron calculations in elastic (single-channel) and inelastic (multi-channel) processes.

Please note: Currently, only elastic processes are implemented.

Link to GRASPC paper in Computer Physics Communications

The main idea is to use as many computational apparatus as they are implemented in GRASP, by adapting them to calculate continuum orbitals wave functions.

GRASPC is entirely transparent for usual calculations in GRASP (bound states and their properties). Only when calculations involving continuum orbital are requested (during the execution of the rmcdhf program), the default flow is changed, and different outputs are produced. The only differences are listed below:

• the default number of points in the radial (NNNP) is increased to 5000 in the src/lib/libmod/parameter_def_M.f90 file, since continuum orbitals have to be calculated on the far-reaching computational grid; for larger systems, it is recommended to magnify the grid size even several times
• rmcdhf: in addition to the other properties, $\langle r^3\rangle$, $\langle r^5\rangle$ and $\langle r^7\rangle$ are calculated for each orbital, since they are part of the model polarization potential
• rwfnplot: generation of the input file for the gnuplot plotting program has been added
• CMakeList.txt: The file has been reorganized so that the additional user settings in the CMakeList.user file are correctly taken into account
• CMakeList.user: file with -fallow-argument-mismatch flag, required for recent versions of the gfortran compiler; this file is not part of the original GRASP repository

All modifications in the original source files are marked in the following way:

! PS      (or # PS in CMakeList.* )
  ... code modified/added
! PS END  (or # PS END in CMakeList.* )

Theoretical background

Generation of the continuum orbital wave function

The basics of the relativistic multiconfiguration Dirac-Hartree-Fock method (RMCDHF) applied to scattering are described in [1]. In short, the scattering system is constructed as $N+1$ electron system, where $N$ electrons are bound (with discrete, negative energy levels), and one electron is from the continuum spectrum, with given (positive) energy $\epsilon$ and quantum number $\kappa$.

To obtain the large $P_{\kappa\epsilon}(r)$ and small $Q_{\kappa\epsilon}(r)$ components of the continuum orbitals, the Dirac-Hartree-Fock equations

$\big(\frac{\textrm{d}}{\textrm{d}r}+\frac{\kappa}{r}\big)P_{\kappa\epsilon}=\big(\frac{2}{\alpha}+\alpha(\epsilon-V-V_{pol})\big)Q_{\kappa\epsilon}-X^{(Q)}$

$\big(\frac{\textrm{d}}{\textrm{d}r}-\frac{\kappa}{r}\big)Q_{\kappa\epsilon}=-\alpha\big(\epsilon-V-V_{pol}\big)P_{\kappa\epsilon}+X^{(P)}$

are solved using the outward integration method implemented in GRASP. In the above equations, $\alpha$ is the fine structure constant, $\epsilon$ is the energy of the scattered particle, $V = V(r)$ is the Coulomb potential, $V_{pol} = V_{pol}(r)$ is a polarization potential, and $X^{(P)}$ and $X^{(Q)}$ are the exchange terms [2].

Please note: From the perspective of the GRASP source code, bound orbital energies are positive. Thus, the continuum electron energy is considered negative. This does not affect the user's actions, for him the continuum electron energy is always positive (or zero in a special case).

If the calculated wave function is to be coupled with some other function (e.g., the bound state), it should be normalized first. The currently implemented per energy continuum wave function normalization is described in [3].

Please note: For calculations of phase shifts and scattering lengths, normalization is usually not required.

Polarization potential

At greater distances, the polarization potential coincides with dipole, quadrupole, and octupole polarization but is limited near the nucleus. It is currently modeled and implemented in the following form:

$V_{pol}\left(r\right)=-\frac{1}{2}\frac{\alpha_d r^2}{\left(r^3+\langle r_0^3\rangle\right)^2} -\frac{1}{2}\frac{\alpha_q r^4}{\left(r^5+\langle r_0^5\rangle\right)^2} -\frac{1}{2}\frac{\alpha_o r^6}{\left(r^7+\langle r_0^7\rangle\right)^2}$,

where $\alpha_d$, $\alpha_q$ and $\alpha_o$ represent the static dipole, quadrupole and octupole polarizabilities, respectively; $\langle r_0^3\rangle$, $\langle r_0^5\rangle$ and $\langle r_0^7\rangle$ are the cut-off parameters.

For atoms, a great source for static dipole polarizabilities is [4] (except for the Livermorium atom, atomic number 116). There is no such a similar source for static quadrupole and octupole polarizabilities, unfortunately.

For atoms, cut-offs $\langle r_0^3\rangle$, $\langle r_0^5\rangle$ and $\langle r_0^7\rangle$ can be taken from bound states calculations, assuming that $\langle r_0\rangle$ is the radius of the outermost orbital of the target atom.

Polarization potential can also be provided in numerical form. It has to be provided in a text file named vpol located in the working directory, containing space-separated pairs $r\ \ V_{pol}(r)$ in rows (in Bohr radii $a_0$, and hartree units, respectively), ordered by increasing $r$ values, e.g.

1.00000E-05 -4.10660E-07
1.05127E-05 -4.31715E-07
1.10517E-05 -4.76698E-07
(...)

Phase shift

Relativistic phase shift $\delta_l$ is obtained by comparing the numerical solution with the asymptotic one at large $r$ [1], determined to the nearest $2\pi$:

$P_\kappa(r)/r = j_l(kr)\cos(\delta_l) - n_l(kr)\sin(\delta_l)$,

where $j_l(kr)$ and $n_l(kr)$ are the spherical Bessel functions of the first and the second kind, respectively, and $k=\sqrt{2\epsilon+\alpha^2\epsilon^2}$ is the wavenumber. Orbital angular momentum quantum number $l$ is related to $\kappa$ such that $\kappa = l$ or $\kappa = -l-1$, whereby $\kappa \neq 0$. Therefore, for the $l$-th partial wave (except the $l=0$ case) we obtain two phase shifts, sometimes referred to $\delta_l^+$ and $\delta_l^-$.

Please note: To correctly determine the phase shift, the continuum orbital has to be calculated far enough from the origin, to ensure that atom-electron interaction can be neglected.

Please note: The phase shift is determined to the nearest $2\pi$, and in many scientific publications is given in the $[−\pi, \pi]$ range. Thus, the calculated value is shifted by $2\pi$ if necessary, giving the second, equivalent value in the output

Scattering length

Scattering length is one of the most useful parameters for describing low-energy electron-atom collisions. It is defined as the radius of a rigid sphere in the zero-energy total cross-section. The sign of a scattering length represents the type of interaction: positive for repulsion and negative for attraction.

The typical method for determining the scattering length $a$ involves analyzing the asymptotic behavior of the wave function:

$a = -\lim_{k \to 0}\frac{\tan(\delta_0)}{k}$,

where $\delta_0$ is the phase shift for $l=0$ (s-wave, $\kappa=-1$), calculated for (energy-dependent) wavenumber tending towards zero in the limit.

In addition to the above method (which is subject to the inaccuracy associated with the non-zero energy), the calculation of the scattering length is implemented as the intersection of the asymptote of the zero-energy ($\epsilon=0$) wave function with the $r$-axis. The details of that approach are described in [5].

Please note: Zero energy wave function, used for accurate scattering length calculations is not a typical oscillating wave function.

More detailed explanations of the theoretical background and computational procedure can be found in [6].

References

1. P. Syty and J.E. Sienkiewicz, Relativistic Multiconfiguration Dirac-Hartree-Fock in the scattering of electrons from argon atoms, J. Phys. B: At. Mol. Opt. Phys. 38 (2005) 2859. https://doi.org/10.1088/0953-4075/38/16/001
2. I.P. Grant, B.J. McKenzie, P.H. Norrington, D.F. Mayers, and N.C. Pyper. An atomic multiconfigurational Dirac-Fock package. Comput. Phys. Commun. 21 (1980) 207–231. https://doi.org/10.1016/0010-4655(80)90041-7
3. R.D. Cowan, The Theory of Atomic Structure and Spectra, University of California Press, Berkeley (1981) pp. 516–521. https://www.ucpress.edu/book/9780520038219/the-theory-of-atomic-structure-and-spectra
4. P. Schwerdtfeger and J.K. Nagle, 2018 Table of static dipole polarizabilities of the neutral elements in the periodic table, Molecular Physics 117 (2019) 9-12. https://doi.org/10.1080/00268976.2018.1535143
5. P. Syty, M.P. Piłat, J.E. Sienkiewicz, Calculation of electron scattering lengths on Ar, Kr, Xe, Rn and Og atoms, J. Phys. B: At. Mol. Opt. Phys. 57 (2024) 175202. https://doi.org/10.1088/1361-6455/ad4fd1
6. P. Syty, M. Piłat, J.E. Sienkiewicz, GRASPC – GRASP package adapted for the generation of continuum orbitals wave functions, Comput. Phys. Commun. 315 (2025) 109691. https://doi.org/10.1016/j.cpc.2025.109691

Current status

Implemented and tested

• Calculations of the continuum wave function of electrons elastically scattered from atoms and ions, with the model or numerical polarization term
• Phase shift calculations (with grid control to ensure correct results)
• Calculations of electronic scattering lengths using the zero energy wave function
• Normalization of the calculated continuum orbital

Changelog

• 2025-03-13

  • co_add_polarization_potential Added: Octupole term in the model polarization potential (this entailed a slight change in the input reception)
  • co_scattering_length Added: Calculation of relative error of straight line determination at the tail, in scattering length calculation (in the zero-energy case)
  • Several minor code and test-case improvements

• 2024-07-18

  • rmcdhf Added: phase shift calculation
  • rwfnestimate Added: Dedicated method for initial estimate of the radial wave function for continuum electron
  • grasptest/continuum Added: new example for electron-argon scattering (d-wave calculation)
  • README Updated: most sections, the most important updates for User guide and Theoretical background
  • rmcdhf Now, the energy of the continuum electron should be entered as a positive value (or zero)
  • Several minor code and test-case improvements

• 2024-06-25 (Initial release)

  • Continuum orbitals generator with scattering length calculations, together with examples for e-Ar and e-Sr scattering

TO-DO list

• Positron scattering (work in progress)
• Validation of user-input parameters (work in progress)
• Inelastic scattering (work at the final concept phase)

User guide

1. Bound states

   Optimize bound states of the selected target (atom/ion) in a usual way (rnucleus => rcsfgenerate => rwfnestimate => rangular => rmcdhf => rsave in the simplest case), or take nuclear data (isodata) and radial wave functions (rwfn.out / .w) files from previous calculations.

   GRASPC may also be used for that calculation, since it works exactly as the original GRASP for bound states.

   Please note: Since only radial orbitals are used (and not the expansion coefficients), Calling the rci code does not improve anything.

2. CSF list

   By invoking rcsfgenerate, create a special rcsf.inp file with only one CSF, where 'core' is the configuration of the target atom or ion, and 'peel' consists of one additional inactive electron. This electron will be treated as a scattered one; its principal quantum number will be ignored, and its quantum number $\kappa$ will be determined from the subshell designation and final $J$ value, resulting from coupling with the core.

   Examples for argon-electron scattering:

   • Core subshells: 1s 2s 2p- 2p 3s 3p- 3p,
   • Additional subshell:
     • 4s(1,i), coupled to $J=1/2$ means s-wave electron of $\kappa=-1$
     • 4p(1,i), coupled to $J=1/2$ means p-wave electron of $\kappa=1$
     • 4p(1,i), coupled to $J=3/2$ means p-wave electron of $\kappa=-2$
     • 4d(1,i), coupled to $J=3/2$ means d-wave electron of $\kappa=2$
     • 4d(1,i), coupled to $J=5/2$ means d-wave electron of $\kappa=-3$, etc.

   Remember to provide the $J$ value as $2J, 2J$ range (e.g. 1,1 for $J=1/2$), and to provide 0 as the number of excitations.

   Example rcsf.inp file for electron-argon scattering (s-wave):

Core subshells:
1s   2s   2p-  2p   3s   3p-  3p
Peel subshells:
4s
CSF(s):
4s ( 1)
    1/2
     1/2+

3. Angular coefficients

   Run rangular as usual (with Full interaction option enabled).

4. Initial estimations of the radial wave functions

   Run rwfnestimate using the previously calculated radial wave functions as initial estimation for the 'core' orbitals (option 1 -- GRASP92 File). To estimate the additional electron designed to be a continuum one, use the new option 5 -- Continuum orbital. Any other method (options 2 - 4) should also work (try them if you encounter convergence problems during the actual calculations in the next step).

5. SCF procedure

   Invoke rmcdhf, then

   • answer n when asked about Default settings
   • answer y when asked about Perform continuum wave function calculations?
   • provide continuum electron energy in hartree (should be positive or zero for accurate scattering length calculation)
   • decide if polarization potential should be included.
     • 0 - do not include polarization potential
     • 1 - include the dipole term of the model potential with default parameters: $\alpha_d$ taken from [2], and cut-off $\langle r_0^3\rangle$ taken from bound state calculations as the size of the outermost orbital; here, the quadrupole and octupole terms are omitted
     • 2 - include model potential with all parameters provided manually by the user; in turn: $\alpha_d$, $\langle r_0^3\rangle$, $\alpha_q$, $\langle r_0^5\rangle$, $\alpha_o$ and $\langle r_0^7\rangle$
     • 3 - include numerical potential from the file named vpol
   • decide (y/n) if the calculated continuum wave function should be normalized
   • answer y when asked Change default speed of light or radial grid parameters?
   • answer y when asked Revise default radial grid parameters?
   • enter new RNT and H values (firstly, the defaults might be kept)
   • enter new HP; this value is the step size ath the tail; use non-zero value to force the linearly-logarithmic grid, which ensures adequate grid density far from the scattering center; 1.0 or less is a good choice for a first-try
   • enter new N; in general, use as a big number as possible to ensure as long grid as possible; 5000 is the default, but even tens of thousands of points might be required in some cases, which would require code modification (NNNP value in src/lib/libmod/parameter_def_M.f90) and recompilation

   Please note: If calculations do not converge (Maximal iterations exceeded), experiment with the other grid parameters, and/or try the other method to estimate the radial wave function for the continuum electron.

   Please note: Only the rmcdhf code has been adapted for continuum orbital calculations, not rmcdhf_mem, rmcdhf_mem_mpi or rmcdhf_mpi.

6. Output

   The calculated continuum orbital wave function will be stored in the rwfn.out file (together with the bound orbitals), and also in a text-formatted file continuum.csp.

   If the grid is long enough and electron energy is not zero, calculated phase shift and scattering length are written to the screen and the rmcdhf.sum summary file. The accuracy of the scattering length strongly depends on the electron energy. (smaller energy means better accuracy since it should be calculated in the $k\to0$ limit).

   If the electron energy is set to zero, only the scattering length is calculated in a more accurate approach, as the intersection of the asymptote of the zero-energy wave function with the r-axis. In that case, additional information is given, which may be helpful for error estimation. The first one is the relative error of straight line determination at the tail. The second one is the percentage difference between the scattering length calculated from the last two points on the grid and the penultimate and the one before the penultimate point. In both cases, a lower value means better accuracy.

Examples

Elastic scattering of electrons from argon atoms

Bound states estimation, continuum orbital wave functions generation (for two different partial waves), phase shift calculations, scattering length calculations using the zero energy approach; calculations with default, numerical and manually entered parameters of polarization potential. The energies of the incident electron were chosen to compare the results with those in Y. Cheng, S. Liu, S. B. Zhang, and Y.-B. Tang, Phys. Rev. A 102, 012824 (2020), https://doi.org/10.1103/PhysRevA.102.012824 easily.

Files in /grasptest/continuum/argon-electron_scattering directory:

• 1_Ar_continuum_electron_s-wave_function - script calculating continuum orbital of $\kappa = -1$ (s-wave, $J=1/2$) and electron energy $\epsilon=0.00124997$ hartree, using the dipole term of the polarization potential with default parameters
• 2_Ar_continuum_electron_d-wave_function- script calculating normalized continuum orbital of $\kappa = 2$ (d-wave, $J=3/2$) and electron energy $\epsilon=0.04499880$ hartree, using polarization potential read from a file
• 3_Ar_electronic_scattering_length - script calculating electronic scattering length using the zero-energy wave function
• Ar_bound_DF - auxiliary script creating the nuclear data and calculating the bound states of argon atom (very simple case without correlations); there's no need to call it by hand; it is invoked automatically by all of the above scripts
• clean - auxiliary script, deleting all of the output files created during program execution
• vpol - auxiliary file with the polarization potential given in numerical form; it is used by the second script above

Scripts realizing the main calculations (beginning with a number) can be run independently and in any order.

The results stay in very good agreement with the results from Cheng when it comes to phase shifts in the first two examples ($\delta_0 = 0.03432$ vs. $0.03999$, and $\delta_2^+ = 0.03241$ vs. $0.03253$), and electronic scattering length in the third example ($a_0 = -1.320$ vs. $-1.394$).

Electronic scattering length of strontium

Calculation of scattering length with default parameters of polarization potential, for a given .w file with optimized bound states.

Files in the /grasptest/continuum/strontium_electronic_scattering_length directory:

• Sr_scattering_length - script calculating electronic scattering length using the zero-energy wave function
• isodata, Sr.w - nuclear data and previously calculated bound states of strontium
• clean - script cleaning the output files

Contributors

Code development and testing; preparing of the documentation; preparing and scripting the test cases; maintaining the repository; co-authoring the publication

• Paweł Syty, Gdańsk University of Technology, pawel.syty@pg.edu.pl

Giving the ideas and theoretical background; proposing and improving the test cases; improving the documentation;

co-authoring the publication

• Józef E. Sienkiewicz, Gdańsk University of Technology
• Michał Piłat, Gdańsk University of Technology

How to cite

P. Syty, M. Piłat, J.E. Sienkiewicz, GRASPC – GRASP package adapted for the generation of continuum orbitals wave functions, Comput. Phys. Commun. 315 (2025) 109691. https://doi.org/10.1016/j.cpc.2025.109691

Installation

Please note: All the installation instructions for the original GRASP are also valid for GRASPC, with the repository address changed to: https://github.com/sylaspg/grasp-continuum.git. Further in this section are original instructions adapted for the GRASPC repository. All the credits goes to the GRASP contributors, who prepared the original installation procedure.

To compile and install GRASPC, first clone this Git repository:

git clone https://github.com/sylaspg/grasp-continuum.git

There are two ways to build GRASPC: either via CMake or via the Makefiles in the source tree. Either works, and you will end up with the GRASPC binaries in the bin/ directory.

CMake is the recommended way to build GRASPC. The Makefile-based workflow is still there to make it smoother to transition from Makefiles to a modern build system.

CMake-based build

The first step with CMake is to create a separate out-of-source build directory. The configure.sh script can do that for you:

cd grasp-continuum/ && ./configure.sh

This will create a build/ directory with the default Release build configuration. However, configure.sh is just a simple wrapper around a cmake call, and if you need more control over the build, you can always invoke cmake yourself (see CMake documentation for more information).

To then compile GRASPC, you need to go into the out-of-source build directory and call make:

cd build/ && make install

Remarks:

• Running make install instructs CMake to actually install the resulting binaries into the conventional bin/ directory at the root of the repository.

  When you run make, the resulting binaries will end up under the build/ directory (specifically in build/bin/). This is useful when developing and debugging, as it allows you to compile many versions of the binaries from the same source tree with different compilation options (e.g. build with debug symbols enabled) by using several out of source build directories.

• With CMake, GRASPC also supports parallel builds, which can be enabled by passing the -j option to make (e.g., make -j4 install to build with four processes).

• The CMake-based build allows running the (non-comprehensive) test suite by calling ctest in the build/ directory. The configuration and source files for the tests are under test//

Makefile-based build

The legacy Makefile-based build can be performed by simply calling the make in the top level directory:

make

In this case, each library and program is compiled in their respective directory under src/, and the built artifacts are stored in the source tree. The resulting binaries and libraries will directly get installed under the bin/ and lib/ directories.

To build a specific library or binary, you can pass the path to the source directory as the Make target:

# build libmod
make src/lib/libmod
# build the rci_mpi binary
make src/appl/rci90_mpi

Note that any necessary library dependencies will also get built automatically.

WARNING: The Makefiles' do not know about the dependencies between the source files, so parallel builds (i.e., calling makewith the-j` option) do not work.

Customizing the build

By default, the Makefile is designed to use gfortran. The variables affecting GRASPC builds are defined and documented at the beginning of the Makefile.

It should never be necessary for the user to modify the Makefile itself. Rather, a Make.user file can be created next to the main Makefile where the build variables can be overridden. E.g., to use the Intel Fortran compiler instead, you may want to create the following Make.user file:

export FC = ifort
export FC_FLAGS = -O3 -save -mkl=sequential
export FC_LD =
export FC_MPI = mpiifort
export OMPI_FC=${FC}

where -mkl=sequential should be set depending on what version of ifort you can access.

Alternatively, to customize the GNU gfortran build to e.g., use a specific version of the compiler, you can create a Make.user file such as

export FC = gfortran-9
export FC_FLAGS = -O3 -fno-automatic
export FC_LD =
export FC_MPI= mpifort
export OMPI_FC=${FC}

To set up a linker search path for the BLAS or LAPACK libraries, you can set FC_LD as follows:

export FC_LD = -L /path/to/blas

The repository also contains the Make.user.gfortran and Make.user.ifort files, which can be used as templates for your own Make.user file.