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	photoionisation (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.

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

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

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

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

void	exampleModule (const IO::InputBlock &input, const Wavefunction &wf)
	Template module: copy this file to write your own module.

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	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	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	normalisation (const IO::InputBlock &input, const Wavefunction &wf)

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

void	polarisability (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.

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

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

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.

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. Used to test Breit, compare sensitivities of different systems to Breit effects, and to compare energy corrections between works. Should not be used for final/accurate breit calculations; see note below.

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.

Note: The interplay between correlations and the Breit effect is very important, and not captured in this perturbative approach. To accurately inlcude Breit corrections to energies and wavefunction, Breit should be added to the wavefunction in the HartreeFock{} input block.

◆ exampleModule()

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

Template module: copy this file to write your own module.

To create a new module:

Copy this file and rename it (e.g., myModule.cpp).
Rename exampleModule to your module name throughout (function declaration, definition, and the Register entry).
Either drop the file anywhere under src/ and recompile (internal), or keep it outside ampsci and add it to your Makefile (external): EXTERNAL_MODULES = path/to/myModule.cpp
Add the module to the Register (see example) - this is what makes your module visible to ampsci and allows it to be run from the input options.

No other files need to be edited: the Register object below adds the module to the list ampsci knows about at compile time.

Every module is a free function with this exact signature, taking the user input block and a pre-solved Wavefunction. Inside the module, use the precomputed wavefunction to calculate any property you require. See example for how to parse input options using IO::InputBlock, and a few basic examples of things to do with the wavefunctions.

◆ 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 ('each').

Frequency-dependent operators are evaluated at each pair's transition frequency (their physical frequency); omega_operator pins them at a fixed frequency instead (rarely meaningful).

Calculations are performed by Amplitudes::matrix_elements; this module only parses input and prints.

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). The single-particle basis is taken from the stored CI integrals (see Wavefunction::CI_integrals), so it need not be re-listed here.

◆ 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 (Flambaum-Ginges 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.

It should be used to test sensitivity of energies to parts of the QED potential, and can be used to compare QED results between works.

For accurate calculations, correlation corrections are required, which are not accounted for in this module. For those cases, add QED into the wavefunction itself using the options in HartreeFock{} input block.

This module also does not calculate QED corrections to matrix elements. You can do that manually using the usual MatrixElements module, and running with/without QED.

Note: Assumes QED is NOT already included in the input wavefunction.