Module Namespace Reference

Detailed Description

Modules are user-defined calculations run after the wavefunction has been solved.

A module is a self-contained C++ function that takes an already-computed atomic wavefunction (Hartree-Fock, MBPT, etc.) and performs some downstream calculation: matrix elements, polarisabilities, hyperfine constants, PNC amplitudes, lifetimes, and so on. Any number of modules can be run from a single input file by adding Module::ModuleName{ ... } blocks; the ampsci driver looks up each name in the registry below and dispatches to the corresponding function.

See Writing custom modules for a more detailed tutorial.

Adding a new module

Adding a module requires editing exactly one new .cpp file. No central list needs to be updated; the build system picks up new .cpp files automatically, and the module self-registers at program startup via a static initialiser.

The minimal template:

// src/Modules/myModule.cpp
#include "IO/InputBlock.hpp"
#include "Modules/Modules.hpp"
#include "Wavefunction/Wavefunction.hpp"
 
namespace Module {
 
// Declare, register, then define below.
void myModule(const IO::InputBlock &input, const Wavefunction &wf);
namespace {
const Register r_myModule{"myModule", "One-line description shown in -m",
                          &myModule};
} // namespace
 
void myModule(const IO::InputBlock &input, const Wavefunction &wf) {
 
  input.check({{"option1", "Description of option1 [default1]"},
               {"option2", "Description of option2 [default2]"}});
  if (input.has_option("help"))
    return;
 
  const auto option1 = input.get("option1", default1);
  // ... physics ...
}
 
} // namespace Module

The first string passed to Register is the name used in input files: Module::myModule{ ... }. The Register variable itself is never referenced by name – it exists purely so its constructor fires before main() and adds the entry to the singleton Registry. The anonymous namespace gives it internal linkage, so different module files can each use the same identifier (r_myModule, or whatever you like) without colliding at link time.

A complete worked example is provided in src/Modules/exampleModule.cpp, intended as a template for new users to copy.

Running a module

In the input file:

Module::myModule{
  option1 = 3.14;
  option2 = true;
}

From the command line, to list available modules or query a module's options:

./ampsci -m                 # list all modules
./ampsci -m myModule        # show options for `myModule`

Classes
struct	ModuleEntry
	One entry in the module registry. More...
struct	Register
	Helper struct: constructing a Register adds a module to the Registry. More...
class	Registry
	Singleton registry of all compiled-in modules. More...

Typedefs
using	ModuleFn = void(*)(const IO::InputBlock &, const Wavefunction &)
	Function-pointer signature shared by every module.

Functions
void	Kionisation (const IO::InputBlock &input, const Wavefunction &wf)
void	photo (const IO::InputBlock &input, const Wavefunction &wf)
void	formFactors (const IO::InputBlock &input, const Wavefunction &wf)
void	tests (const IO::InputBlock &input, const Wavefunction &wf)
	Runs basic numerical tests on HF orbitals and basis/spectrum.
void	writeOrbitals (const IO::InputBlock &input, const Wavefunction &wf)
	Writes all orbital sets to disk for plotting.
void	continuum (const IO::InputBlock &input, const Wavefunction &wf)
	Computes continuum orbitals and tests orthogonality.
void	Breit (const IO::InputBlock &input, const Wavefunction &wf)
	Breit corrections to HF energies and matrix elements.
void	dcp (const IO::InputBlock &input, const Wavefunction &wf)
void	exampleModule (const IO::InputBlock &input, const Wavefunction &wf)
void	HFAnomaly (const IO::InputBlock &input, const Wavefunction &wf)
	Bohr-Weisskopf effect and hyperfine anomaly.
void	b_plot (const IO::InputBlock &input, const Wavefunction &wf)
	b moments KS(R) and KL(R) as a function of nuclear rms radius.
void	BW_effect (const std::vector< DiracSpinor > &valence, const DiracOperator::TensorOperator *const h0, const ExternalField::DiagramRPA *const rpa0, const DiracOperator::TensorOperator *const h, const ExternalField::DiagramRPA *const rpa)
void	tune_Rmag (const DiracSpinor &Fv, const double eps_target, const std::optional< IO::InputBlock > &hfs_options, const Wavefunction &wf, const DiracOperator::TensorOperator *const h0, const ExternalField::DiagramRPA *const rpa0)
void	BW_screening_factor (const Wavefunction &wf, const DiracOperator::TensorOperator *const h0, const ExternalField::DiagramRPA *const rpa0, const DiracOperator::TensorOperator *const h, const ExternalField::DiagramRPA *const rpa)
std::array< double, 3 >	fit (std::vector< double > &xd, std::vector< double > &yd, std::size_t terms)
std::pair< std::array< double, 3 >, std::array< double, 3 > >	b_moments (const std::string &iso, const DiracSpinor &v, double R0_fm, int max_power)
void	fieldShift (const IO::InputBlock &input, const Wavefunction &wf)
	Isotope field shift using TDHF/RPA and MBPT.
void	fieldShift_direct (const IO::InputBlock &input, const Wavefunction &wf)
	Isotope field shift by direct Hartree-Fock recalculation.
void	ladder (const IO::InputBlock &input, const Wavefunction &wf)
void	check_L_symmetry (const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, const std::vector< DiracSpinor > &valence, const Coulomb::QkTable &qk, bool include_L4, const Angular::SixJTable &sj, const Coulomb::LkTable *const lk=nullptr)
void	matrixElements (const IO::InputBlock &input, const Wavefunction &wf)
	Matrix elements of any operator for HF/Brueckner valence states.
void	CI_matrixElements (const IO::InputBlock &input, const Wavefunction &wf)
	Matrix elements of any operator between CI many-body states.
void	structureRad (const IO::InputBlock &input, const Wavefunction &wf)
	Structure radiation and normalisation corrections to matrix elements.
void	normalisation (const IO::InputBlock &input, const Wavefunction &wf)
	Normalisation correction via derivative of correlation potential.
void	runModules (const IO::InputBlock &input, const Wavefunction &wf)
	Iterate over the input blocks and run any that are modules.
void	runModule (const IO::InputBlock &input, const Wavefunction &wf)
	Run a single module.
void	list_modules ()
	Print the list of compiled-in modules (name + description).
void	pnc (const IO::InputBlock &input, const Wavefunction &wf)
void	polarisability (const IO::InputBlock &input, const Wavefunction &wf)
void	dynamicPolarisability (const IO::InputBlock &input, const Wavefunction &wf)
void	transitionPolarisability (const IO::InputBlock &input, const Wavefunction &wf)
void	QED (const IO::InputBlock &input, const Wavefunction &wf)
	QED corrections to HF energies.

---

Class Documentation

Module::ModuleEntry

struct Module::ModuleEntry

Class Members
string	name
string	description
ModuleFn	function

Typedef Documentation

ModuleFn

using Module::ModuleFn = typedef void (*)(const IO::InputBlock &, const Wavefunction &)

Function-pointer signature shared by every module.

Function Documentation

tests()

void Module::tests	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Runs basic numerical tests on HF orbitals and basis/spectrum.

Tests orthonormality, Hamiltonian expectation values, radial boundaries, and basis/spectrum comparisons (energies, E1, HFS, sum rules). Basis/spectrum tests are skipped if the corresponding sets are empty.

writeOrbitals()

void Module::writeOrbitals	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Writes all orbital sets to disk for plotting.

continuum()

void Module::continuum	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Computes continuum orbitals and tests orthogonality.

Solves for continuum states at each specified energy for l = 0..max_l. Optionally computes matrix elements with a chosen operator, with or without RPA. Orbitals can be written to file for plotting.

Breit()

void Module::Breit	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Breit corrections to HF energies and matrix elements.

Calculates first-order (one-body) Breit corrections to core and valence energies, with and without HF relaxation. Decomposes into Gaunt, retardation, and frequency-dependent contributions.

Optionally computes second-order Breit corrections using the Johnson formula (requires a basis). The second-order term includes both HF one-body (de(B1,2)) and Sigma two-body (de(B2,2)) contributions.

Note

Assumes Breit is NOT already included in the input wavefunction.

Warning

Not compatible with Brueckner (correlation potential) orbitals.

To calculate Breit corrections including MBPT, just use ampsci manually (i.e., calculate once with and once without Breit). This is because Breit must be included into the basis/Green's function when forming correlation potential.

HFAnomaly()

void Module::HFAnomaly	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Bohr-Weisskopf effect and hyperfine anomaly.

Computes the Bohr-Weisskopf (BW) effect: the fractional change in the hyperfine constant due to a finite nuclear magnetisation distribution, defined:

\[ A \equiv A_0(1 + \varepsilon_{\rm BW}) \]

where \( A_0 \) is the pointlike result and \( A \) uses a finite distribution (Ball/Fermi/Sp model, set via hfs_options). Note: Some places define \( \varepsilon \) with opposite sign.

Optionally:

• Includes RPA corrections (diagram method only, required a basis).
• Computes BW screening factors \( x \) and \( \eta_{sp} \) (relative to hydrogenic values), as in PhysRevA.105.052802.
• Computes b moments in the \( K_S(R) \) and \( K_L(R) \) polynomial expansions (see b_power).
• Tunes the magnetic radius to reproduce an experimental BW effect (eps_target).
• Computes the differential hyperfine anomaly \({}^1\Delta^2\) between two isotopes (A2), optionally tuning both magnetic radii to match a target (1D2_target).

Note

Requires valence states. For RPA, a basis is also needed.

Since this module runs automatically, calculations should be checked/confirmed using manual calculation with the Module::MatrixElements module

b_plot()

void Module::b_plot	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

b moments KS(R) and KL(R) as a function of nuclear rms radius.

Computes the b-moment polynomial coefficients in the expansions

\[ K_S(R_m) = \sum_n k_{2n}^{(S)} R_m^{2n}, \quad K_L(R_m) = \sum_n k_{2n}^{(L)} R_m^{2n}, \]

across a range of nuclear rms radii \( [r_1, r_2] \). At each radius, HF is re-solved and b_moments fits the radial integrals up to b_power. Numerical uncertainties are estimated by varying grid density, inner grid boundary, and fit radius.

arXiv:2603.20090

Output is a table: r_rms, KS coefficients and uncertainties, KL coefficients and uncertainties.

fieldShift()

void Module::fieldShift	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Isotope field shift using TDHF/RPA and MBPT.

Calculates the field shift constants F, G2, G4 by varying the nuclear charge radius and fitting the resulting energy shifts:

\[ \delta E = F \, \delta\langle r^2 \rangle + G_2 (\delta\langle r^2 \rangle)^2 + G_4 \, \delta\langle r^4 \rangle \]

Uses TDHF/RPA to include core polarisation. Fits are performed via GSL.

Warning

Old module, not well tested. Use with caution. G definitions may not match literature or expectations.

fieldShift_direct()

void Module::fieldShift_direct	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Isotope field shift by direct Hartree-Fock recalculation.

Calculates the field shift constant F by directly re-solving Hartree-Fock for a range of nuclear charge radii and fitting:

\[ \delta E = F \, \delta\langle r^2 \rangle \]

Core relaxation (equivalent to RPA) is optionally included. Does not recompute the correlation potential Sigma.

Warning

Old module, not well tested. Use with caution.

matrixElements()

void Module::matrixElements	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Matrix elements of any operator for HF/Brueckner valence states.

Computes matrix elements \( \redmatel{a}{h}{b} \) of any registered operator operator between valence states, with optional RPA corrections. Output can be reduced matrix elements, stretched-state matrix elements (with \( m = j \)), or hyperfine constants \( A, B, \ldots \).

Optionally:

• Includes RPA corrections (TDHF, diagram, or basis method).
• Includes core-state matrix elements.
• Uses the spectrum instead of HF valence states.
• Solves RPA at a fixed frequency or at each transition frequency (for frequency-dependent operators or 'each').

Note

For hyperfine operators (hfs, MLVP), the default output is HFS constants rather than reduced matrix elements.

CI_matrixElements()

void Module::CI_matrixElements	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Matrix elements of any operator between CI many-body states.

Computes matrix elements between CI (configuration interaction) many-body wavefunctions. Loops over angular momentum \( J \) and parity, applying selection rules. Optionally includes:

• RPA corrections (diagram method).
• Structure radiation and normalisation corrections (via StructureRadiation{}).
• Frequency-dependent operators, solved at each transition frequency or at a fixed value.

Note

Requires CI wavefunctions to be computed (via the CI{} module).

structureRad()

void Module::structureRad	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Structure radiation and normalisation corrections to matrix elements.

Computes second-order MBPT corrections to matrix elements using the structure radiation (SR) framework. Includes:

• SR "top", "bottom", and "centre" diagrams.
• Normalisation of states.
• Brueckner orbital (BO) corrections (if no Brueckner orbitals present).

The total corrected matrix element is

\[ t_{ab}^{\rm tot} = t_{ab}^{(0)} + \delta V_{ab} + \delta t_{ab}^{\rm SR} + \delta t_{ab}^{\rm Norm} + \delta t_{ab}^{\rm BO}. \]

Optionally reads/writes Coulomb \( Q^k \) integrals to file for speed.

Note

Requires a basis (B-splines). For large bases, the QkTable option gives ~10x speedup at the cost of memory.

normalisation()

void Module::normalisation	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Normalisation correction via derivative of correlation potential.

Computes the normalisation factor for each valence state from the energy derivative of the correlation potential,

\[ \delta t_{ab}^{\rm Norm} = \frac{1}{2} \left( \langle v | \frac{\partial \Sigma}{\partial \varepsilon} | v \rangle_a + \langle v | \frac{\partial \Sigma}{\partial \varepsilon} | v \rangle_b \right) t_{ab}, \]

evaluated numerically via a finite difference \( \delta \varepsilon \). This is an alternative to the perturbative SR normalisation.

Note

Requires Brueckner orbitals (correlation potential Sigma must be present).

runModules()

void Module::runModules	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Iterate over the input blocks and run any that are modules.

Scans input for blocks whose name matches Module::* and dispatches each to runModule in turn. Called once by the ampsci driver after the wavefunction has been solved.

Parameters

input	Top-level input block.
wf	Solved wavefunction passed to each module.

runModule()

void Module::runModule	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Run a single module.

Looks up the module by name in the Registry and calls its function. If the name is not found, prints an error message together with a nearest-match suggestion and a list of available modules.

Parameters

input	Input block for this module (typically Module::Name{}).
wf	Solved wavefunction.

list_modules()

void Module::list_modules

(

)

Print the list of compiled-in modules (name + description).

Invoked by ./ampsci -m. Iterates over Registry::entries() and prints each module's name and description (word-wrapped).

QED()

void Module::QED	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

QED corrections to HF energies.

Calculates first-order QED corrections (Ginges-Flambaum radiative potential) to core and valence energies, with and without HF relaxation. Decomposes into Uehling (vacuum polarisation), self-energy (high- and low-frequency), and magnetic form-factor contributions.

Note

Assumes QED is NOT already included in the input wavefunction.