SplineBasis Namespace Reference

Detailed Description

Constructs spinor/orbital basis using B-splines.

Uses Bsplines to form a set of B-spline orbitals (using method from [1] "Derevianko", [2] "Shabaev", [3] "Johnson", or [4] "Fischer", or the "Hybrid" of [1] and [3]; see SplineType). Diagonalises B-splines over Hamiltonian to produce a set of basis orbitals.

[1] K. Beloy, A. Derevianko, Comput. Phys. Commun. 179, 310 (2008). [2] V. M. Shabaev, I. I. Tupitsyn, V. A. Yerokhin, G. Plunien, and G. Soff, Phys. Rev. Lett. 93, 130405 (2004). [3] W. Johnson, S. Blundell, J. Sapirstein, Phys. Rev. A 37, 307 (1988). [4] C. F. Fischer and F. A. Parpia, Phys. Lett. A 179, 198 (1993). See also: Bachau et al., Reports Prog. Phys. 64, 1815 (2001).

By default, uses the Derevianko Dual Kinetic Balance (DKB) basis.

If \(\{|i\rangle\}\) are the set of spline orbitals (of a given angular symmetry; \(2N\) of them, except Fischer with \(\kappa>0\), see spline_params), and

\[ H_{ij} = \langle{S_i}|\hat H_{\rm HF}|{S_j}\rangle \, , \qquad S_{ij} = \langle{S_i|S_j}\rangle. \]

The eigenvalue problem:

\[ H_{ij}p_j = \epsilon S_{ij}p_j, \]

is solved, yielding one eigenvalue \(\epsilon\) (with eigenvector \(p\)) for each spline orbital: roughly half are positive-energy solutions, and half negative-energy solutions.

Note: form_basis() does not store the eigenvectors, instead, it expands the basis orbitals and stores them on the regular grid (coordinate space). i.e., for each eigenvalue, n, the corresponding basis orbital is:

\[ |{n}\rangle = \sum_i p_i |i\rangle\, \]

Classes
struct	Parameters
	Input parameters for B-spline basis construction (see form_basis). More...
struct	SplineParams
	Type-dependent spline selection for a given kappa (see spline_params). More...
struct	SplineSpinors
	The spline spinors for one kappa (see form_spline_basis) More...
struct	SplineSpinorValues
	Values at r of the two spline spinors built from one B-spline B_i. More...

Enumerations
enum class	SplineType {   Derevianko , Shabaev , Johnson , Fischer ,   Hybrid }
	B-spline construction method (see details). More...

Functions
std::vector< double >	default_r0 (SplineType type)
	Default first internal knot, r0, for each l, for the given type.
std::vector< DiracSpinor >	form_basis (const Parameters &params, const Wavefunction &wf, const bool correlationsQ=false)
	Forms and returns the basis orbitals (expanded in terms of splines).
double	check (const std::vector< DiracSpinor > &basis, const std::vector< DiracSpinor > &orbs, bool print_warning=true)
	Checks basis against orbitals: normality, orthogonality, and energies.
SplineParams	spline_params (SplineType type, int kappa, std::size_t n_states)
	Returns SplineParams for the given type, kappa, and number of states.
SplineSpinorValues	spline_spinor_values (double B, double dB, double d2B, int kappa, double r, double kb)
	Values at r of the f-like and g-like spline spinors built from B_i.
SplineSpinors	form_spline_basis (const int kappa, const std::size_t n_states, const std::size_t k_spl, const double r0_spl, const double rmax_spl, std::shared_ptr< const Grid > rgrid, const double alpha, SplineType itype=SplineType::Derevianko)
	Forms the underlying spline basis (which is not kept).
std::pair< LinAlg::Matrix< double >, LinAlg::Matrix< double > >	fill_Hamiltonian_matrix (const SplineSpinors &splines, const Wavefunction &wf, const bool correlationsQ=false, SplineType itype=SplineType::Derevianko)
	Calculates and returns the Hamiltonian matrix H_ij and overlap S_ij.
void	add_JohnsonBoundary (LinAlg::Matrix< double > *Aij, const int kappa, const double alpha, const std::vector< double > &f_rmax, const std::vector< double > &g_rmax)
	Applies Johnson/Notre-Dame boundary conditions to Hamiltonian matrix.
void	add_rmaxBoundary (LinAlg::Matrix< double > *Aij, const double alpha, const std::vector< double > &f_rmax, const std::vector< double > &g_rmax)
	Applies the rmax (MIT bag) boundary term only: f(R) = -g(R).
void	expand_basis_orbitals (std::vector< DiracSpinor > *basis, std::vector< DiracSpinor > *basis_positron, const std::vector< DiracSpinor > &spl_basis, const int kappa, const int max_n, int max_n_positron, const LinAlg::Vector< double > &e_values, const LinAlg::Matrix< double > &e_vectors, const Wavefunction &wf)
	Expands basis orbitals in terms of spline orbitals.
std::vector< double >	sumrule_TKR (const std::vector< DiracSpinor > &basis, const std::vector< double > &r, bool print=false)
	TKR sum rule (basis test); should be zero (needs -ve energy states).
std::vector< double >	sumrule_DG (int nDG, const std::vector< DiracSpinor > &basis, const Grid &gr, double alpha, bool print)
	Drake-Goldman sum rules (basis test); should be zero.
std::pair< double, double >	r_completeness (const DiracSpinor &Fv, const std::vector< DiracSpinor > &basis, const Grid &gr, bool print=false)
	Measures radial completeness of the basis for orbital Fv.
std::string_view	parseSplineType (SplineType type)
	Returns the canonical name string for a SplineType.
auto	parseSplineType (std::string_view type)
	Parses a string to SplineType (case-insensitive).
void	r_completeness_header ()
	Prints the title and column headers for r_completeness output.

---

Class Documentation

SplineBasis::SplineSpinors

struct SplineBasis::SplineSpinors

Class Members
vector< DiracSpinor >	basis {}	Spline spinors on the radial grid: f-like, then g-like.
vector< DiracSpinor >	d_basis {}	Radial derivatives (d/dr of both f and g) of the spline spinors.
vector< double >	f_rmax {}	f and g of each spline spinor exactly at rmax (for boundary terms)
vector< double >	g_rmax {}

SplineBasis::SplineSpinorValues

struct SplineBasis::SplineSpinorValues

Class Members
double	f_f	f-like spinor (f = B_i): f, g, and their radial derivatives
double	g_f
double	df_f
double	dg_f
double	f_g	g-like spinor (g = B_i): f, g, and their radial derivatives
double	g_g
double	df_g
double	dg_g

Enumeration Type Documentation

SplineType

enum class SplineBasis::SplineType

strong

B-spline construction method (see details).

• Derevianko (Reno): dual kinetic balance (DKB) splines; the first |kappa| and the last spline are dropped, so the B-spline parts vanish at both ends of the cavity (the kinetic-balance parts need not). K. Beloy, A. Derevianko, Comput. Phys. Commun. 179, 310 (2008).
• Shabaev: DKB splines; the first and last spline are dropped, for all kappa, with no boundary terms (as published: dropping more at r = 0 for |kappa| >= 2 was found to make practically no difference). V. M. Shabaev, I. I. Tupitsyn, V. A. Yerokhin, G. Plunien, and G. Soff, Phys. Rev. Lett. 93, 130405 (2004).
• Johnson (Notre Dame): plain splines for f and g (no kinetic balance), all kept; the boundary conditions f(0) = 0 and f(R) = -g(R) (MIT bag) are imposed through boundary terms in the Hamiltonian matrix. W. R. Johnson, S. A. Blundell, J. Sapirstein, Phys. Rev. A 37, 307 (1988).
• Fischer: plain splines (no kinetic balance); the first and last spline are dropped (f = g = 0 at both ends), with no boundary terms. For kappa > 0, the second spline is also dropped from f (f'(0) = 0), which removes the spurious state found by Fischer and Parpia for kappa = 1. C. F. Fischer and F. A. Parpia, Phys. Lett. A 179, 198 (1993).
• Hybrid: DKB splines; the first |kappa| are dropped (as for Derevianko), but the last is kept, and the Johnson rmax (MIT bag) boundary term f(R) = -g(R) is imposed (see add_rmaxBoundary). Differs from Derevianko only in the boundary condition at rmax.

See parseSplineType for the input strings.

Warning

Fischer: bound states are accurate, but the pseudo-continuum (cavity) states are not reliable: they disagree with those of the other types, and for kappa > 0 some lie below the corresponding free-atom energies. May not be suitable for MBPT; intended for testing only.

Function Documentation

default_r0()

std::vector< double > SplineBasis::default_r0

(

SplineType

type

)

Default first internal knot, r0, for each l, for the given type.

Element l is the first knot for orbital angular momentum l; the last element is used for all higher l (see form_basis).

• Dual-kinetic-balance types (Derevianko, Shabaev, Hybrid): 1e-4 for all l
• Types without kinetic balance (Johnson, Fischer): 3e-4 for s, 1e-4 for p, 1e-3 for d and f, 1e-2 for higher l. As in Johnson et al. (Table III), the first knot is moved out for higher l, since these orbitals are negligible close to the nucleus, which improves the knot distribution. The s and p values give the most accurate hyperfine constants and Drake-Goldman sum rules for Cs (see the basis tests); in particular, a first s knot at 1e-4 degrades the s-state hyperfine constants.

form_basis()

std::vector< DiracSpinor > SplineBasis::form_basis	(	const Parameters &	params,
		const Wavefunction &	wf,
		const bool	correlationsQ = false
	)

Forms and returns the basis orbitals (expanded in terms of splines).

Members of params:

• states = which states to keep e.g., "25spd10f" (up to n=25 for s,p,d-states, and up to n=10 for f states)
• n: Number of splines (the underlying spline set may be larger; see spline_params)
• k: k order of the B-splines. NB: must have \(k\geq l_{\rm max}+3\) [1]
• r0: first internal knot, for each l: r0[l] is used for orbital angular momentum l, and the last entry for all higher l (so a single value applies to all l). Empty: default_r0(type). Values below the first grid point are replaced by it.
• rmax: last knot (basis orbitals only non-zero before this point)
• positron: basis string giving the number of negative-energy states to keep for each kappa (e.g., "9spdf" keeps 9 of each). Those closest to -2mc^2 are kept. They are labelled with n continuing on from the positive-energy solutions (so {n,kappa} stays unique), are flagged (DiracSpinor::negativeEnergyStateQ), and are appended to the end of the basis std::vector. Empty string discards them. Must be a subset of states. A warning is also printed if any label is too large for DiracSpinor::nk_index (n > 256; e.g., many splines with many negative-energy states kept).
• type: B-spline construction method (see SplineType)
• verbose: print construction details to screen
• orthogonalise: not used here (see Wavefunction::formBasis)

Parameters

params	Basis parameters (see above)
wf	Wavefunction object: needed to form Hartree-Fock Hamiltonian
correlationsQ	If true, includes the correlation potential (Sigma) of wf in the Hamiltonian (if it exists)

Returns

Basis orbitals: electron states for each kappa (in the order of params.states, increasing n), followed by any negative-energy states.

Note

This function calls the below functions, they rarely need to be called explicitly, unless you are trying to do something different to usual.

check()

double SplineBasis::check	(	const std::vector< DiracSpinor > &	basis,
		const std::vector< DiracSpinor > &	orbs,
		bool	print_warning = true
	)

Checks basis against orbitals: normality, orthogonality, and energies.

For each orbital in orbs with a matching basis state (same n and kappa), finds the worst \( ||\langle a|b\rangle| - 1| \) and relative energy difference, \( |(\varepsilon_b - \varepsilon_a)/\varepsilon_b| \). Also finds the worst orthonormality between orbs and basis states of the same kappa (see DiracSpinor::check_ortho). Prints each, and lists basis states with opposite sign to the matching orbital.

Parameters

basis	Basis orbitals
orbs	Reference orbitals (e.g., core or valence)
print_warning	Flag printed results worse than 1e-3

Returns

Largest of the three deviations; 0 (prints nothing) if orbs is empty.

spline_params()

SplineParams SplineBasis::spline_params	(	SplineType	type,
		int	kappa,
		std::size_t	n_states
	)

Returns SplineParams for the given type, kappa, and number of states.

Of the n_spl underlying B-splines, n_states are kept for each of f and g (indices [imin_f, imax_f) and [imin_g, imax_g)):

• Johnson : all; no kinetic balance
• Fischer : drops the first and the last; no kinetic balance. For kappa > 0, the second is also dropped from f (f'(0) = 0), which removes the spurious kappa = 1 state found by Fischer and Parpia; f then has n_states - 1 splines
• Hybrid : drops the first |kappa|; kinetic balance
• Shabaev : drops the first and the last; kinetic balance
• Derevianko : drops the first |kappa| and the last; kinetic balance

spline_spinor_values()

SplineSpinorValues SplineBasis::spline_spinor_values	(	double	B,
		double	dB,
		double	d2B,
		int	kappa,
		double	r,
		double	kb
	)

Values at r of the f-like and g-like spline spinors built from B_i.

See form_spline_basis: the f-like spinor has f = B_i and g = kb (B_i' + kappa B_i/r); the g-like spinor has g = B_i and f = kb (B_i' - kappa B_i/r), where kb = lambda alpha/2 is the kinetic balance factor (zero without kinetic balance).

Parameters

B	B_i(r)
dB	B_i'(r)
d2B	B_i''(r)
kappa	Dirac quantum number
r	Radius (> 0)
kb	Kinetic balance factor, lambda alpha/2

form_spline_basis()

SplineSpinors SplineBasis::form_spline_basis	(	const int	kappa,
		const std::size_t	n_states,
		const std::size_t	k_spl,
		const double	r0_spl,
		const double	rmax_spl,
		std::shared_ptr< const Grid >	rgrid,
		const double	alpha,
		SplineType	itype = SplineType::Derevianko
	)

Forms the underlying spline basis (which is not kept).

Forms the spline spinors for given kappa from the B-splines \(B_i\) selected by spline_params: \(N_f\) "f-like" followed by \(N_g\) "g-like" (both equal to n_states, except Fischer with kappa > 0),

\[ S_i = \begin{pmatrix} B_i \\ \frac{\lambda\alpha}{2}\left(B_i' + \frac{\kappa}{r}B_i\right) \end{pmatrix}, \qquad S_{N_f+i} = \begin{pmatrix} \frac{\lambda\alpha}{2}\left(B_i' - \frac{\kappa}{r}B_i\right) \\ B_i \end{pmatrix}, \]

where \(\lambda\) is the kinetic balance prefactor (SplineParams). Knots are distributed log-linearly between r0_spl and rmax_spl, with the end points (0 and rmax) repeated k times: only the first spline is non-zero at r = 0, and only the last at r = rmax.

Parameters

kappa	Dirac quantum number
n_states	Number of splines, N, of each type (f-like and g-like; N-1 f-like for Fischer with kappa > 0)
k_spl	Order of the B-splines
r0_spl	First internal knot
rmax_spl	Last knot (splines are zero beyond this point)
rgrid	Radial grid the splines are stored on
alpha	Fine-structure constant
itype	B-spline construction method

Returns

The N_f + N_g spline spinors, their radial derivatives (d/dr of both f and g), and the values of f and g of each exactly at rmax (which is generally between grid points).

fill_Hamiltonian_matrix()

std::pair< LinAlg::Matrix< double >, LinAlg::Matrix< double > > SplineBasis::fill_Hamiltonian_matrix	(	const SplineSpinors &	splines,
		const Wavefunction &	wf,
		const bool	correlationsQ = false,
		SplineType	itype = SplineType::Derevianko
	)

Calculates and returns the Hamiltonian matrix H_ij and overlap S_ij.

\( H_{ij} = \langle S_i|\hat H|S_j\rangle \) and \( S_{ij} = \langle S_i|S_j\rangle \), over the spline basis. The Hamiltonian is that of wf: the local part (see Wavefunction::H0ab), plus Hartree-Fock exchange (if non-local), Breit (if included in wf), and Sigma (if correlationsQ and wf has Sigma). For Johnson and Hybrid types, also adds the boundary terms (see add_JohnsonBoundary and add_rmaxBoundary).

Parameters

splines	Spline spinors (see form_spline_basis)
wf	Wavefunction: provides the Hamiltonian
correlationsQ	Include the correlation potential (Sigma) of wf
itype	B-spline construction method (sets boundary terms)

Returns

Pair {H_ij, S_ij}.

add_JohnsonBoundary()

void SplineBasis::add_JohnsonBoundary	(	LinAlg::Matrix< double > *	Aij,
		const int	kappa,
		const double	alpha,
		const std::vector< double > &	f_rmax,
		const std::vector< double > &	g_rmax
	)

Applies Johnson/Notre-Dame boundary conditions to Hamiltonian matrix.

Adds the boundary term of the action, S' [Eq. (10) of Johnson et al.], together with the surface term \(\frac{c}{2}[fg]_0^R\) by which the symmetric kinetic matrix (Wavefunction::H0ab) differs from the antisymmetrised form of Johnson et al. [Eq. (7)], so that f(0) = 0 and f(R) = -g(R) [MIT bag, P(R) = Q(R); g = -Q here] become the natural boundary conditions. Only the first spline is non-zero at r = 0 (B_1(0) = 1), and f(R), g(R) are those of the spline spinors at rmax (see add_rmaxBoundary):

\[ S' + \frac{c}{2}[fg]_0^R = \frac{c}{2} f(0)^2 \ [\kappa<0] \ {\rm or}\ c^2 f(0)^2 \ [\kappa>0] + \frac{c}{4}\left[f(R)^2 - g(R)^2\right] + \frac{c}{2} f(R) g(R), \]

(the f(0)g(0) terms cancel). The matrix gets its second derivatives. The cancellation is exact analytically; on the radial grid it holds only to the accuracy of the quadrature at the cavity ends. W. R. Johnson, S. A. Blundell, J. Sapirstein, Phys. Rev. A 37, 307 (1988).

Parameters

Aij	Hamiltonian matrix, 2N x 2N (N f-like splines, then N g-like); updated in place
kappa	Dirac quantum number
alpha	Fine-structure constant
f_rmax	f(rmax) of each spline spinor (see SplineSpinors)
g_rmax	g(rmax) of each spline spinor

add_rmaxBoundary()

void SplineBasis::add_rmaxBoundary	(	LinAlg::Matrix< double > *	Aij,
		const double	alpha,
		const std::vector< double > &	f_rmax,
		const std::vector< double > &	g_rmax
	)

Applies the rmax (MIT bag) boundary term only: f(R) = -g(R).

The r = R part of add_JohnsonBoundary (no r = 0 term); also used alone for the Hybrid type. With \(f(R) = \sum_i p_i u_i\) and \(g(R) = \sum_i p_i w_i\), where \(u_i\), \(w_i\) are f and g of spline spinor i at R, adds the second derivatives of

\[ \frac{c}{4}\left[f(R)^2 - g(R)^2\right] + \frac{c}{2} f(R) g(R): \quad \frac{c}{2}\left(u_iu_j - w_iw_j + u_iw_j + w_iu_j\right). \]

Without kinetic balance (Johnson), only the last f-like (u = 1) and last g-like (w = 1) spinors are non-zero at R. With kinetic balance (Hybrid), the kinetic-balance parts of the last two splines of each set also contribute, so these are included.

Parameters

Aij	Hamiltonian matrix (f-like spinors, then g-like); updated in place
alpha	Fine-structure constant
f_rmax	f(rmax) of each spline spinor (see SplineSpinors)
g_rmax	g(rmax) of each spline spinor

expand_basis_orbitals()

void SplineBasis::expand_basis_orbitals	(	std::vector< DiracSpinor > *	basis,
		std::vector< DiracSpinor > *	basis_positron,
		const std::vector< DiracSpinor > &	spl_basis,
		const int	kappa,
		const int	max_n,
		int	max_n_positron,
		const LinAlg::Vector< double > &	e_values,
		const LinAlg::Matrix< double > &	e_vectors,
		const Wavefunction &	wf
	)

Expands basis orbitals in terms of spline orbitals.

Basis orbitals are formed from the eigenvectors of the Hamiltonian. Keeps the electron states up to n = max_n (labelled n = l+1, l+2, ..., from the bottom of the positive-energy branch up), and the max_n_positron negative-energy states closest to the gap (labelled with n continuing on from the positive-energy solutions, and flagged).

Each orbital is normalised, with sign chosen so that f is positive at small r.

Parameters

basis	Electron states are appended to this
basis_positron	Negative-energy states are appended to this
spl_basis	Spline basis the eigenvectors are expressed in
kappa	Dirac quantum number
max_n	Maximum n of electron states kept
max_n_positron	Number of negative-energy states kept
e_values	Eigenvalues (ascending)
e_vectors	Eigenvectors: e_vectors[i] corresponds to e_values[i]
wf	Wavefunction: provides grid and alpha

Note

States whose norm differs from 1 by more than 1e-4 are treated as spurious: a warning is printed, and the orbital is set to zero (but is still kept).

sumrule_TKR()

std::vector< double > SplineBasis::sumrule_TKR	(	const std::vector< DiracSpinor > &	basis,
		const std::vector< double > &	r,
		bool	print = false
	)

TKR sum rule (basis test); should be zero (needs -ve energy states).

Thomas-Kuhn-Reiche sum rule. With a the first state in basis, for each l up to the maximum l in basis:

\[ \sum_n \frac{w_n}{2l+1} (\varepsilon_n - \varepsilon_a) \langle a|r|n\rangle^2 = 0, \]

summed over basis states with \(\kappa_n = l\) ( \(w_n = l\)) and \(\kappa_n = -l-1\) ( \(w_n = l+1\)), where \(\langle a|r|n\rangle\) is the radial integral.

Parameters

basis	Basis orbitals (including negative-energy states)
r	Radial grid points
print	Print electron and negative-energy contributions, for each l

Returns

The sum, for each l; empty if basis is empty.

sumrule_DG()

std::vector< double > SplineBasis::sumrule_DG	(	int	nDG,
		const std::vector< DiracSpinor > &	basis,
		const Grid &	gr,
		double	alpha,
		bool	print
	)

Drake-Goldman sum rules (basis test); should be zero.

Must include negative-energy states in basis. For each kappa in basis (a: the first basis state with that kappa):

\[ \frac{1}{2j_a+1}\sum_n (\varepsilon_n - \varepsilon_a)^{m} |\langle a||r||n\rangle|^2 \]

where m = nDG, compared with its exact value: \(\langle a|r^2|a\rangle\) (m=0) or 0 (m=1). For m=2, the sum is scaled by \(\alpha^2/3\), and compared with 1. Kappas with l equal to the maximum l in basis are skipped (no l+1 states to sum over).

Parameters

nDG	Power of the energy difference, m: 0, 1, or 2
basis	Basis orbitals (including negative-energy states)
gr	Radial grid
alpha	Fine-structure constant
print	Print sum, exact value, and difference, for each kappa

Returns

Sum minus exact value, for each kappa included. Kappas not present in basis (e.g., p states for "20s5d") are skipped.

r_completeness()

std::pair< double, double > SplineBasis::r_completeness	(	const DiracSpinor &	Fv,
		const std::vector< DiracSpinor > &	basis,
		const Grid &	gr,
		bool	print = false
	)

Measures radial completeness of the basis for orbital Fv.

Sums over basis states n with the same kappa as Fv:

\[ \sum_n \langle v|r|n\rangle\langle n|r^{-1}|v\rangle = 1, \qquad \sum_n \langle v|r|n\rangle\langle n|r|v\rangle = \langle v|r^2|v\rangle, \]

where the matrix elements are radial integrals.

Parameters

Fv	Orbital to test
basis	Basis orbitals
gr	Radial grid
print	Print one line of results (see r_completeness_header)

Returns

Relative deviation of each sum from its exact value.

parseSplineType() [1/2]

std::string_view SplineBasis::parseSplineType ( SplineType  type )

inline

Returns the canonical name string for a SplineType.

parseSplineType() [2/2]

auto SplineBasis::parseSplineType ( std::string_view  type )

inline

Parses a string to SplineType (case-insensitive).

• Derevianko (default; aliases: reno, dkb)
• Shabaev
• Johnson (aliases: nd, notredame, notre-dame)
• Fischer
• Hybrid

See SplineType for a description of each.

Returns

The SplineType; Derevianko if type is empty.

Warning

An unrecognised type prints a warning and returns Derevianko.

r_completeness_header()

void SplineBasis::r_completeness_header ( )

inline

Prints the title and column headers for r_completeness output.