Detailed Description

Many-body perturbation theory.

Namespaces
namespace	Sigma2
	Functions for each Sigma2 diagram; called by Sk_vwxy.

Classes
class	CorrelationPotential
	Correlation potential, Sigma, for valence states. More...

class	Feynman
	Direct correlation potential diagrams via the Feynman technique. More...

struct	FeynmanOptions
	Options for the Feynman method (see Feynman) More...

class	Goldstone
	Second-order correlation potential, Sigma, from Goldstone diagrams. More...

class	RadialMatrix
	Radial matrix: stored on a sub-grid (defined by i0, stride, size). More...

struct	rgrid_params
	Sigma sub-grid parameters: range (aB) and stride (of full grid) More...

struct	SigmaData
	Stored Sigma matrix, with the parameters it was formed with. More...

struct	SigmaLData
	Ladder correlation potential matrix, Sigma_L, for a single state. More...

class	SpinorMatrix
	Radial Dirac spinor matrix. More...

class	StructureRad
	Structure Radiation and Normalisation of states, using diagram method. More...

Typedefs
using	RMatrix = RadialMatrix< double >

using	ComplexRMatrix = RadialMatrix< std::complex< double > >

using	GMatrix = SpinorMatrix< double >

using	ComplexGMatrix = SpinorMatrix< std::complex< double > >

using	ComplexDouble = std::complex< double >

Enumerations
enum class	SigmaMethod { Goldstone , Feynman }
	Method used to form the direct Sigma diagrams. More...

enum class	HoleParticle { exclude , include , include_k0 }
	Options for including the hole-particle interaction. More...

enum class	Screening { exclude , include }
	Options for including Screening. More...

enum class	GreenStates { both , core , excited }
	Which states to include in the Green's function. More...

enum class	SigmaLMethod { basis , ratio , direct }
	Method used to construct the ladder correlation potential, Sigma_L. More...

enum class	Denominators { RS , Fermi , Fermi0 , DFK , BW }
	Type of energy denominators: DFK, BW, RS, Fermi, Fermi0. More...

Functions
double	fit_tail_coefficient (const GMatrix &Sigma, const DiracSpinor &Fv, std::size_t n_fit=10, double f_cut=1.0e-2)
	Fits the coefficient a_eff of the large-r polarisation tail of Sigma.

void	add_tail (DiracSpinor *SigmaFn, const GMatrix &Sigma, const DiracSpinor &Fn, double a_eff)
	Adds the large-r polarisation tail of Sigma to `SigmaFn`.

double	best_omre (const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &valence, bool print=false)
	Best real part of the frequency, omre, for the direct diagram.

SigmaLMethod	parseSigmaLMethod (const std::string &method)
	Converts string (name) to SigmaLMethod enum (case-insensitive).

std::string	parseSigmaLMethod (SigmaLMethod method)
	Converts SigmaLMethod enum to string (name)

double	Lkmnij (int k, const DiracSpinor &m, const DiracSpinor &n, const DiracSpinor &i, const DiracSpinor &j, const Coulomb::QkTable &qk, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, bool include_L4, const Angular::SixJTable &SJ, const Coulomb::LkTable *const Lk=nullptr, std::optional< double > e_i={}, std::optional< double > e_m={})
	Full ladder integral summed over all diagrams.

double	L1 (int k, const DiracSpinor &m, const DiracSpinor &n, const DiracSpinor &i, const DiracSpinor &j, const Coulomb::QkTable &qk, const std::vector< DiracSpinor > &excited, const Angular::SixJTable &SJ, const Coulomb::LkTable *const Lk=nullptr, std::optional< double > e_i={})
	Particle–particle ladder diagram L1.

double	L4 (int k, const DiracSpinor &m, const DiracSpinor &n, const DiracSpinor &i, const DiracSpinor &j, const Coulomb::QkTable &qk, const std::vector< DiracSpinor > &core, const Angular::SixJTable &SJ, const Coulomb::LkTable *const Lk=nullptr, std::optional< double > e_m={})
	Core–core (hole–hole) ladder diagram L4.

double	L2 (int k, const DiracSpinor &m, const DiracSpinor &n, const DiracSpinor &i, const DiracSpinor &j, const Coulomb::QkTable &qk, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, const Angular::SixJTable &SJ, const Coulomb::LkTable *const Lk=nullptr, std::optional< double > e_j={}, std::optional< double > e_m={})
	Particle–hole ladder diagram L2.

void	fill_Lk_mnib (Coulomb::LkTable *lk, const Coulomb::QkTable &qk, const std::vector< DiracSpinor > &excited, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &i_orbs, bool include_L4, const Angular::SixJTable &sjt, int max_k=-1, bool print=true)
	Fills the ladder integral table for all new index combinations.

GMatrix	Sigma_ladder (const DiracSpinor &v, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, const std::vector< DiracSpinor > &projection, const Coulomb::QkTable &qk, const Coulomb::LkTable *lk, const Angular::SixJTable &sjt, bool include_L4=false, double r0=1.0e-4, double rmax=30.0, std::size_t stride=4, bool include_G=false)
	Ladder-diagram correction to the correlation potential, Sigma_L(e_v).

DiracSpinor	Lkv_mnia (int k, const DiracSpinor &v, const DiracSpinor &m, const DiracSpinor &n, const DiracSpinor &a, const Coulomb::QkTable &qk, const Coulomb::YkTable &yk, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, bool include_L4, const Angular::SixJTable &SJ, const Coulomb::LkTable *const Lk=nullptr)
	Vertex (ket) form of ladder integral L^k_mnia over external index i.

DiracSpinor	Lkv_inab (int k, const DiracSpinor &v, const DiracSpinor &n, const DiracSpinor &a, const DiracSpinor &b, const Coulomb::QkTable &qk, const Coulomb::YkTable &yk, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, bool include_L4, const Angular::SixJTable &SJ, const Coulomb::LkTable *const Lk=nullptr)
	Vertex (ket) form of ladder integral L^k_inab over external index i.

GMatrix	Sigma_ladder_direct (const DiracSpinor &v, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, const Coulomb::QkTable &qk, const Coulomb::YkTable &yk, const Coulomb::LkTable *lk, const Angular::SixJTable &sjt, bool include_L4=false, double r0=1.0e-4, double rmax=30.0, std::size_t stride=4, bool include_G=false)
	Ladder correlation potential Sigma_L via the direct method.

void	update_Lk_mnib (Coulomb::LkTable lk, const Coulomb::QkTable &qk, const std::vector< DiracSpinor > &excited, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &update_i, bool include_L4, const Angular::SixJTable &sjt, const Coulomb::LkTable const lk_prev, double a_damp, bool print)
	Updates the ladder integral table with L(Q,Q) -> L(Q,Q+L)

bool	write_SigmaL (const std::string &fname, const std::vector< SigmaLData > &SLs, const Grid &grid)
	Writes Sigma_L (ladder) matrices to binary file.

std::vector< SigmaLData >	read_SigmaL (const std::string &fname, const std::shared_ptr< const Grid > &grid)
	Reads Sigma_L (ladder) matrices from binary file.

std::pair< int, int >	k_minmax_L (const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d)
	Min and max k (multipolarity) allowed for ladder integral L^k_abcd.

std::pair< int, int >	k_minmax_L (int kap_a, int kap_b, int kap_c, int kap_d)
	Min and max k allowed for ladder integral L^k_abcd (kappa version).

double	L3 (int k, const DiracSpinor &m, const DiracSpinor &n, const DiracSpinor &i, const DiracSpinor &j, const Coulomb::QkTable &qk, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, const Angular::SixJTable &SJ, const Coulomb::LkTable *const Lk=nullptr, std::optional< double > e_i={})
	Exchange partner of L2; equals L2 with m,n and i,j swapped.

template<typename Qintegrals , typename QorLintegrals >
double	de_valence (const DiracSpinor &v, const Qintegrals &qk, const QorLintegrals &lk, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited)
	Second-order (or ladder) correction to the valence energy.

template<typename Qintegrals , typename Lintegrals >
double	de_valence_w (const DiracSpinor &v, const Qintegrals &qk, const Lintegrals &lk, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, const Angular::SixJTable *sj=nullptr)
	Ladder (or MBPT2) valence energy, antisymmetrising the FIRST integral.

template<typename Qintegrals , typename QorLintegrals >
double	de_core (const Qintegrals &qk, const QorLintegrals &lk, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited)
	Second-order (or ladder) correction to the core energy.

void	ladder (const IO::InputBlock &input, const Wavefunction &wf)
	Driver for the ladder-diagram calculation: the Ladder{} input block.

template<typename T >
bool	equal (const RadialMatrix< T > &lhs, const RadialMatrix< T > &rhs)
	Checks if two matrices are equal (to within parts in 10^12)

template<typename T >
double	max_element (const RadialMatrix< T > &a)
	returns maximum element (by abs)

template<typename T >
double	max_delta (const RadialMatrix< T > &a, const RadialMatrix< T > &b)
	returns maximum difference (abs) between two matrices

template<typename T >
double	max_epsilon (const RadialMatrix< T > &a, const RadialMatrix< T > &b)
	returns maximum relative difference, \|(aij-bij)/(aij+bij)\|, of two matrices

std::string	parse_Denominators (Denominators d)
	Returns string representation of Denominators enum.

Denominators	parse_Denominators (std::string_view s)
	Parses string to Denominators enum (case-insensitive); DFK if unrecognised.

double	leg_de (Denominators denominators, double et_bar, double et, double ei_bar, double ei, double es, double E0)
	External-leg part of a Sigma_2 energy denominator.

double	e_bar (int kappa_v, const std::vector< DiracSpinor > &excited)
	Returns energy of first excited state matching a given kappa.

bool	Sk_vwxy_SR (int k, const DiracSpinor &v, const DiracSpinor &w, const DiracSpinor &x, const DiracSpinor &y)
	Selection rule for Sk_vwxy.

std::pair< int, int >	k_minmax_S (const DiracSpinor &v, const DiracSpinor &w, const DiracSpinor &x, const DiracSpinor &y)
	Minimum and maximum k allowed by selection rules for Sk_vwxy.

std::pair< int, int >	k_minmax_S (int twoj_v, int twoj_w, int twoj_x, int twoj_y)
	Overload taking 2j values directly.

double	Sk_vwxy (int k, const DiracSpinor &v, const DiracSpinor &w, const DiracSpinor &x, const DiracSpinor &y, const Coulomb::QkTable &qk, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, const Angular::SixJTable &SixJ, Denominators denominators=Denominators::DFK, const std::vector< double > &fk={}, double E0=0.0)
	Reduced two-body Sigma (2nd-order correlation) operator matrix element.

Coulomb::LkTable	calculate_Sk (const std::string &filename, const std::vector< DiracSpinor > &external, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, const Coulomb::QkTable &qk, int max_k, bool exclude_wrong_parity_box, Denominators denominators, bool no_new_integrals=false, const std::vector< double > &fk={}, double E0=0.0)
	Calculates (or reads in) a table of two-body Sigma_2 matrix elements.

std::vector< double >	average_hk (const Coulomb::LkTable &Sk, const Coulomb::QkTable &qk, const std::vector< DiracSpinor > &external, int max_k=-1)
	Average Sigma_2 correction ratios, h_k, for each multipole k.

template<class CoulombIntegral >
double	Sigma_vw (const DiracSpinor &v, const DiracSpinor &w, const CoulombIntegral &qk, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, int max_l_internal=99, std::optional< double > ev=std::nullopt)
	Matrix element of the 1-body Sigma (2nd-order correlation) operator.

template<class CoulombIntegral >
std::pair< double, double >	Sigma_vw_direct_exchange (const DiracSpinor &v, const DiracSpinor &w, const CoulombIntegral &qk, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, int max_l_internal=99, std::optional< double > ev=std::nullopt)
	Direct and exchange parts of Sigma_vw(), returned separately.

template<class CoulombIntegral >
double	dSigma_dE_vw (const DiracSpinor &v, const DiracSpinor &w, const CoulombIntegral &qk, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, double ev, int max_l_internal=99, double delta=0.01)
	Energy derivative of the one-body correlation correction.

template<typename T >
bool	equal (const SpinorMatrix< T > &lhs, const SpinorMatrix< T > &rhs)
	Checks if two matrices are equal (to within parts in 10^12)

template<typename T >
double	max_element (const SpinorMatrix< T > &a)
	returns maximum element (by abs)

template<typename T >
double	max_delta (const SpinorMatrix< T > &a, const SpinorMatrix< T > &b)
	Returns maximum difference (abs) between two matrices.

template<typename T >
double	max_epsilon (const SpinorMatrix< T > &a, const SpinorMatrix< T > &b)
	Returns maximum relative difference (abs) between two matrices.

Variables
constexpr std::size_t	sk_array_size = 32

const auto	vroot = [](auto x) { return std::sqrt(x); }

Class Documentation

◆ MBPT::FeynmanOptions

struct MBPT::FeynmanOptions

Class Members
Screening	screening {Screening::exclude}	Include all-order screening of the Coulomb interaction.
HoleParticle	hole_particle {HoleParticle::exclude}	Include hole-particle interaction in the polarisation operator.
int	max_l_internal {6}	Maximum l for internal lines (at least the maximum l in the core)
double	omre {-0.3}	Real part of the frequency, for integration along Im(w) (see best_omre)
double	w0 {0.01}	First non-zero point of the frequency grid along Im(w). u=0 is always included. The [0, w0] panel is a trapezoid: w0 must be small compared to the pole distances (~\|omre\|) so the integrand is flat across the panel.
double	w_ratio {1.5}	Ratio for the logarithmic frequency grid along Im(w)
bool	complex_green {true}	Method for the Green's function at complex energy. false : solve at Re(en), extend to complex energy via Dyson (resolvent) true : solve the Dirac equation directly at complex energy (Dyson still used at large Im(en); see Feynman::green)

◆ MBPT::rgrid_params

struct MBPT::rgrid_params

Class Members
double	r0 {1.0e-4}
double	rmax {30.0}
size_t	stride {4}

◆ MBPT::SigmaData

struct MBPT::SigmaData

Class Members
int	kappa	Dirac quantum number of the state Sigma was formed for.
double	en	Energy Sigma was formed at.
SpinorMatrix< double >	Sigma	Sigma matrix, on the sub-grid.
int	n {0}	Principal quantum number of the state Sigma was formed for (0: none given)
double	lambda {1.0}	Scaling factor, Sigma -> lambda*Sigma; not applied to the stored matrix.
vector< double >	fk {}	fk screening factors used for this Sigma (empty if none)
double	a_eff {0.0}	Coefficient of large-r tail, Sigma = a_eff/r^4 beyond sub-grid (0: none)

◆ MBPT::SigmaLData

struct MBPT::SigmaLData

Class Members
int	kappa	Dirac quantum number of the state Sigma_L was formed for.
int	n	Principal quantum number of the state Sigma_L was formed for.
double	en	Energy Sigma_L was evaluated at.
GMatrix	SL	Sigma_L matrix, on its own sub-grid.

Enumeration Type Documentation

◆ SigmaMethod

enum class MBPT::SigmaMethod

strong

Method used to form the direct Sigma diagrams.

Goldstone (basis sum) or Feynman (Green's function). Exchange diagrams always use Goldstone.

◆ HoleParticle

enum class MBPT::HoleParticle

strong

Options for including the hole-particle interaction.

exclude : not included
include : included, all k
include_k0 : included, k=0 term only

◆ Screening

enum class MBPT::Screening

strong

Options for including Screening.

◆ GreenStates

enum class MBPT::GreenStates

strong

Which states to include in the Green's function.

both : core and excited (full Hartree-Fock Green's function)
core : core states only
excited : excited states only (orthogonalised to the core)

◆ SigmaLMethod

enum class MBPT::SigmaLMethod

strong

Method used to construct the ladder correlation potential, Sigma_L.

basis : project onto the basis; requires extending Qk (slow) [Sigma_ladder()]
ratio : no projection; rescale each Sigma(2) term by L/Q [Sigma_ladder(), empty projection basis]
direct : no projection; open the external line exactly (ladder vertex) [Sigma_ladder_direct()]

◆ Denominators

enum class MBPT::Denominators

strong

Type of energy denominators: DFK, BW, RS, Fermi, Fermi0.

DFK : Dzuba-Flambaum-Kozlov convention (Brillouin-Wigner-like, with the target-state energy approximated by the lowest configuration). The external leg belonging to the target state is evaluated at the Fermi level (lowest state for its kappa in excited spectrum); the external leg appearing in the intermediate state keeps its actual orbital energy. Retains state dependence, with no danger of accidental enhancement.
BW : Brillouin-Wigner: the denominator is E0 - E_intermediate, where E0 is the total valence energy of the target CI level, and E_intermediate is the total zeroth-order energy of the many-body state between the two Coulomb vertices: the sum of orbital energies of the particles present, minus the holes (e.g., diagram 'a': E_int = e_v + e_y + e_n - e_a). In practice, the target-state leg of DFK is replaced by (E0 - e_s), where e_s is the other valence orbital in that diagram's intermediate state. DFK is this with E0 approximated by the leading configuration, E0 -> e_bar_target + e_s. Requires E0.
RS : Use actual orbital energies for both external legs. May be danger of accidental enhancement.
Fermi : Both external legs evaluated at the Fermi level for their kappa.
Fermi0 : As above, but assume Fermi level for all kappas the same. These often cancel, so there is no (excited-excited) term in denominator (except diagram d). Fine, since the remaining core-excited always dominates.

In each case, each diagram is averaged with its bra-ket partner, 0.5*(1/de + 1/de'), so that S^k (and hence the CI matrix) is symmetric. This is the Hermitian effective Hamiltonian, correct to this order in PT. For Fermi0 and BW the two partners coincide identically (BW uses the target energy for both bra and ket), so the average does nothing.

Energy for internal legs (hole-particle) always actual orbitals.

Function Documentation

◆ fit_tail_coefficient()

double MBPT::fit_tail_coefficient	(	const GMatrix &	Sigma,
		const DiracSpinor &	Fv,
		std::size_t	n_fit = `10`,
		double	f_cut = `1.0e-2`
	)

Fits the coefficient a_eff of the large-r polarisation tail of Sigma.

Large-r tail: \( \Sigma(r) \to a_{\rm eff}/r^4 \).

Averages \( r^4 (\Sigma F_v)(r)/F_v(r) \) (large component) over n_fit sub-grid points: the outermost ones where \( |F_v| \geq f_{\rm cut} \max|F_v| \) (stopping at a node). Further out, Fv is exponentially small and the ratio is dominated by numerical noise in Sigma. In the fit region (outside the core, typically a few to ~15 aB) \( r^4 (\Sigma F_v)/F_v \) is already close to constant.

Returns: a_eff; zero if there are no points in the fit region.

Note: Can be numerically unstable, should be treated as a test

◆ add_tail()

void MBPT::add_tail	(	DiracSpinor *	SigmaFn,
		const GMatrix &	Sigma,
		const DiracSpinor &	Fn,
		double	a_eff
	)

Adds the large-r polarisation tail of Sigma to SigmaFn.

Adds \( (a_{\rm eff}/r^4) F_n \) for r beyond the Sigma sub-grid (where \( \Sigma F_n \) is zero). Treated as a local potential. Acts on the small component only if Sigma includes g parts. Does nothing if a_eff is zero.

◆ best_omre()

double MBPT::best_omre	(	const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	valence,
		bool	print = `false`
	)

Best real part of the frequency, omre, for the direct diagram.

The value furthest from any pole of the u=0 integrand. With \( \Delta = \varepsilon_{\rm exc} - \varepsilon_{\rm core} \) (lowest excited minus highest core energy), the window \( (-\Delta, 0) \) is free of true poles; inside it lie only the weak fictitious poles at core-core energy differences (imperfect core subtraction in the polarisation-loop Gex) and, for non-lowest valence states, the small valence-valence differences (valence-line G). Only the lowest few physical excited states enter, so the valence list suffices: no basis/spectrum needed (a typical valence energy is assumed if the list is empty). One omre serves all valence states.

Parameters

core	Core states
valence	Valence states (lowest sets the window)
print	Print Delta, the in-window poles, and the chosen omre

Returns: Midpoint of the largest pole-free gap.

◆ parseSigmaLMethod() [1/2]

SigmaLMethod MBPT::parseSigmaLMethod ( const std::string & method )

Converts string (name) to SigmaLMethod enum (case-insensitive).

Also accepts "full" (basis) and "Dzuba" (ratio). Warns and defaults to ratio if unknown.

◆ parseSigmaLMethod() [2/2]

std::string MBPT::parseSigmaLMethod ( SigmaLMethod method )

Converts SigmaLMethod enum to string (name)

◆ Lkmnij()

double MBPT::Lkmnij	(	int	k,
		const DiracSpinor &	m,
		const DiracSpinor &	n,
		const DiracSpinor &	i,
		const DiracSpinor &	j,
		const Coulomb::QkTable &	qk,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		bool	include_L4,
		const Angular::SixJTable &	SJ,
		const Coulomb::LkTable *const	Lk = `nullptr`,
		std::optional< double >	e_i = `{}`,
		std::optional< double >	e_m = `{}`
	)

Full ladder integral summed over all diagrams.

Computes

\[ L^k_{mnij} = L1^k_{mnij} + L2^k_{mnij} + L3^k_{mnij} [+ L4^k_{mnij}] \]

where \( L3^k_{mnij} = L2^k_{nmji} \) and \( L4 \) involves core–core intermediate states. Lk points to the ladder table from the previous iteration; pass nullptr on the first iteration.

Parameters

k	Multipole rank
m,n	Excited (particle) orbitals
i,j	Core (hole) or valence orbitals
qk	Coulomb \( Q^k \) integral table
core	Core orbitals
excited	Excited orbitals
include_L4	Include the core–core diagram L4
SJ	6j symbol table
Lk	Ladder table from previous iteration (nullptr on first)
e_i	Optional: used in place of i.en() in energy denominators
e_m	Optional: used in place of m.en() in energy denominators

Returns: \( L^k_{mnij} \)

Note: To evaluate at a fixed external energy (for a correlation potential), use e_i or e_m (for external line in the i or m slot): the energy only enters the denominators, never the integral lookups.

◆ L1()

double MBPT::L1	(	int	k,
		const DiracSpinor &	m,
		const DiracSpinor &	n,
		const DiracSpinor &	i,
		const DiracSpinor &	j,
		const Coulomb::QkTable &	qk,
		const std::vector< DiracSpinor > &	excited,
		const Angular::SixJTable &	SJ,
		const Coulomb::LkTable *const	Lk = `nullptr`,
		std::optional< double >	e_i = `{}`
	)

Particle–particle ladder diagram L1.

\[ L1^k_{mnij} = \sum_{rs,ul} A^{kul}_{mnrsij} \frac{Q^u_{mnrs}\,(Q+L)^l_{rsij}}{\epsilon_{ij} - \epsilon_{rs}} \]

with the angular coefficient

\[ A^{kul}_{mnrsij} = (-1)^{m+n+r+s+i+j+1}\,[k] \sixj{m}{i}{k}{l}{u}{r}\sixj{n}{j}{k}{l}{u}{s} \]

Intermediate states \( r,s \) run over excited orbitals.

Parameters

k	Multipole rank
m,n	Excited (particle) orbitals
i,j	Core (hole) or valence orbitals
qk	Coulomb \( Q^k \) integral table
excited	Excited orbitals
SJ	6j symbol table
Lk	Ladder table from previous iteration (nullptr on first)
e_i	Optional: used in place of i.en() in energy denominator

Returns: \( L1^k_{mnij} \)

◆ L4()

double MBPT::L4	(	int	k,
		const DiracSpinor &	m,
		const DiracSpinor &	n,
		const DiracSpinor &	i,
		const DiracSpinor &	j,
		const Coulomb::QkTable &	qk,
		const std::vector< DiracSpinor > &	core,
		const Angular::SixJTable &	SJ,
		const Coulomb::LkTable *const	Lk = `nullptr`,
		std::optional< double >	e_m = `{}`
	)

Core–core (hole–hole) ladder diagram L4.

\[ L4^k_{mnij} = \sum_{cd,ul} A^{klu}_{mncdij} \frac{Q^u_{cdij}\,(Q+L)^l_{mncd}}{\epsilon_{cd} - \epsilon_{mn}} \]

with the angular coefficient \( A \) as in L1(). Intermediate states run over core orbitals only, making this the hole–hole counterpart of the particle–particle diagram L1. Enable via include_L4 in Lkmnij().

Parameters

k	Multipole rank
m,n	Excited (particle) orbitals
i,j	Core (hole) or valence orbitals
qk	Coulomb \( Q^k \) integral table
core	Core orbitals
SJ	6j symbol table
Lk	Ladder table from previous iteration (nullptr on first)
e_m	Optional: used in place of m.en() in energy denominator

Returns: \( L4^k_{mnij} \)

◆ L2()

double MBPT::L2	(	int	k,
		const DiracSpinor &	m,
		const DiracSpinor &	n,
		const DiracSpinor &	i,
		const DiracSpinor &	j,
		const Coulomb::QkTable &	qk,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		const Angular::SixJTable &	SJ,
		const Coulomb::LkTable *const	Lk = `nullptr`,
		std::optional< double >	e_j = `{}`,
		std::optional< double >	e_m = `{}`
	)

Particle–hole ladder diagram L2.

\[ L2^k_{mnij} = \sum_{rc,ul} (-1)^{k+u+l+1} A^{klu}_{mjcrin} \frac{Q^u_{cnir}\,(Q+L)^l_{mrcj}}{\epsilon_{cj} - \epsilon_{mr}} \]

Intermediate states: \( r \) runs over excited, \( c \) over core. The diagram \( L3 \) is the exchange partner \( L3^k_{mnij} = L2^k_{nmji} \).

Parameters

k	Multipole rank
m,n	Excited (particle) orbitals
i,j	Core (hole) or valence orbitals
qk	Coulomb \( Q^k \) integral table
core	Core orbitals
excited	Excited orbitals
SJ	6j symbol table
Lk	Ladder table from previous iteration (nullptr on first)
e_j	Optional: used in place of j.en() in energy denominator
e_m	Optional: used in place of m.en() in energy denominator

Returns: \( L2^k_{mnij} \)

◆ fill_Lk_mnib()

void MBPT::fill_Lk_mnib	(	Coulomb::LkTable *	lk,
		const Coulomb::QkTable &	qk,
		const std::vector< DiracSpinor > &	excited,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	i_orbs,
		bool	include_L4,
		const Angular::SixJTable &	sjt,
		int	max_k = `-1`,
		bool	print = `true`
	)

Fills the ladder integral table for all new index combinations.

Iterates over all combinations of excited pairs \( (m,n) \) and orbitals in i_orbs, computing \( L^k_{mnib} \) (with \( b \) in core) and storing results in lk. Only calculates new integrals. Only lowest-order.

Parameters

lk	Output ladder table (written in place)
qk	Coulomb \( Q^k \) integral table
excited	Excited orbitals
core	Core orbitals
i_orbs	Orbitals for the \( i \) index
include_L4	Include core-core diagram L4
sjt	6j symbol table
max_k	Maximum multipolarity; -1 uses qk.max_k()
print	Print progress and timing to screen

◆ Sigma_ladder()

GMatrix MBPT::Sigma_ladder	(	const DiracSpinor &	v,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		const std::vector< DiracSpinor > &	projection,
		const Coulomb::QkTable &	qk,
		const Coulomb::LkTable *	lk,
		const Angular::SixJTable &	sjt,
		bool	include_L4 = `false`,
		double	r0 = `1.0e-4`,
		double	rmax = `30.0`,
		std::size_t	stride = `4`,
		bool	include_G = `false`
	)

Ladder-diagram correction to the correlation potential, Sigma_L(e_v).

By projection (or, if projection is empty, by the L/Q ratio method).

With a non-empty projection basis, forms the ladder correlation potential by projecting the discrete ladder integrals onto the projection states of kappa_v. The exchange is folded into the Coulomb vertex via \( W = Q + P \) (as in de_valence_w):

\[ \Sigma_L = \sum_{i,amn,k} |W^k_{\cdot amn}\rangle\, \frac{L^k_{mn,i,a}} {[k][j_v]\,(\epsilon_v+\epsilon_a-\epsilon_m-\epsilon_n)}\, \langle i| + \sum_{i,nab,k} |W^k_{\cdot nab}\rangle\, \frac{L^k_{i,n,a,b}} {[k][j_v]\,(\epsilon_v+\epsilon_n-\epsilon_a-\epsilon_b)}\, \langle i| , \]

(particle-particle (a+b) and particle-hole (c+d) diagrams). The bra index \( i \) runs over the projection states of kappa_v (approximating completeness). The ladder integrals are computed on-the-fly via Lkmnij() evaluated at the fixed external energy \( \epsilon_v \) (via the e_i/e_m energy overrides, since the energy enters only the denominators). Exception: for the valence state itself, the stored table entries are used directly - they are already at the correct energy, so the i = v term is essentially free.

If projection is empty, instead uses the ratio method (following V. A. Dzuba, Phys. Rev. A 78, 042502 (2008)): no projection; each term of the regular second-order correlation potential (cf. Goldstone::Sigma_both) is rescaled by the scalar ratio of ladder to Coulomb integrals:

\[ \Sigma_L = \sum_{amn,k} |Q^k_{\cdot amn}\rangle\, \frac{L^k_{mnva}/Q^k_{mnva}} {[k][j_v]\,(\epsilon_v+\epsilon_a-\epsilon_m-\epsilon_n)}\, \langle W^k_{\cdot amn}| + \sum_{nab,k} |Q^k_{\cdot nab}\rangle\, \frac{L^k_{vnab}/Q^k_{vnab}} {[k][j_v]\,(\epsilon_v+\epsilon_n-\epsilon_a-\epsilon_b)}\, \langle W^k_{\cdot nab}| . \]

By construction the diagonal reproduces the ladder energy exactly, <v|Sigma_L|v> = de_valence_w(v), except for the dropped terms (see note). The off-diagonal (radial) structure is approximate: each term keeps the shape of the corresponding second-order term, rescaled by a scalar. All integrals come straight from the stored tables (the L^k entries with i = v are already at the valence energy), so nothing is computed on-the-fly: much faster than projection.

The sub-grid (r0, rmax, stride) defaults match Wavefunction::formSigma.

Parameters

v	Valence state (basis version; must be in the qk/lk tables)
core	Core (hole) orbitals
excited	Excited orbitals
projection	Projection basis {\|i>}; states of kappa_v are used. If empty, the ratio method is used instead (no projection)
qk	Converged Coulomb \( Q^k \) table
lk	Converged ladder \( L^k \) table. For projection, it is the internal-rung table forwarded to Lkmnij (nullptr for L(Q,Q)=L^(1)); for ratio, it supplies the L^k integrals (nullptr gives Sigma_L = 0)
sjt	6j symbol table
include_L4	Include core–core diagram in on-the-fly Lkmnij (projection only; the ratio method never re-computes L)
r0,rmax,stride	Sub-grid parameters
include_G	Include the lower (g) component of Sigma_L

Returns: Sigma_L as a coordinate-space GMatrix

Note: Ratio method: terms with \( Q^k_{mnva} = 0 \) (but \( L^k \ne 0 \)) cannot be rescaled and are dropped; inherent to the method. Note that L^k exists at both parities of k (see k_minmax_L) while Q^k exists only at Coulomb parity, so the exchange-only (wrong-parity) k channels are always dropped here; the projection and direct methods include them.

◆ Lkv_mnia()

DiracSpinor MBPT::Lkv_mnia	(	int	k,
		const DiracSpinor &	v,
		const DiracSpinor &	m,
		const DiracSpinor &	n,
		const DiracSpinor &	a,
		const Coulomb::QkTable &	qk,
		const Coulomb::YkTable &	yk,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		bool	include_L4,
		const Angular::SixJTable &	SJ,
		const Coulomb::LkTable *const	Lk = `nullptr`
	)

Vertex (ket) form of ladder integral L^k_mnia over external index i.

Returns the radial spinor \( |L^k_{mn \cdot a}\rangle \) (kappa of v) satisfying

\[ \langle x|L^k_{mn\cdot a}\rangle = L^k_{mnxa}(\epsilon_v) \]

for any \( x \) with \( \kappa_x = \kappa_v \), with the external energy fixed at \( \epsilon_v \) (as the e_i override in Lkmnij).

In each diagram the external line attaches to a single bare Coulomb line; that line is opened exactly as a radial function (Qkv_bcd). The exception is the L-part of the internal (Q+L) rung in L1 and L3, where the external line attaches to a dressed rung: for that piece we set i = v (scalar coefficient times |v>), which is exact at x = v and is the same level of treatment the scalar table gives it (entries exist only for stored orbitals).

Parameters

k	Multipole rank
v	Valence state: fixes kappa and the external energy, and supplies \|v> for the dressed-rung (i = v) piece
m,n	Excited (particle) orbitals
a	Core (hole) orbital
qk	Coulomb \( Q^k \) integral table (internal lines)
yk	Yk table spanning core+excited (opens the external line)
core	Core orbitals
excited	Excited orbitals
include_L4	Include core–core diagram L4
SJ	6j symbol table
Lk	Ladder table for the internal (Q+L) rung (nullptr: lowest order)

Returns: \( |L^k_{mn \cdot a}\rangle \)

Note: All internal lines are Hartree-Fock basis states (from the tables); only the external line is left open. This is the correct structure for acting on Brueckner orbitals.

◆ Lkv_inab()

DiracSpinor MBPT::Lkv_inab	(	int	k,
		const DiracSpinor &	v,
		const DiracSpinor &	n,
		const DiracSpinor &	a,
		const DiracSpinor &	b,
		const Coulomb::QkTable &	qk,
		const Coulomb::YkTable &	yk,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		bool	include_L4,
		const Angular::SixJTable &	SJ,
		const Coulomb::LkTable *const	Lk = `nullptr`
	)

Vertex (ket) form of ladder integral L^k_inab over external index i.

The external index i is in the m-slot. Returns the radial spinor \( |L^k_{\cdot nab}\rangle \) (kappa of v) satisfying

\[ \langle x|L^k_{\cdot nab}\rangle = L^k_{xnab}(\epsilon_v) \]

for any \( x \) with \( \kappa_x = \kappa_v \), with the external energy fixed at \( \epsilon_v \) (as the e_m override in Lkmnij).

Mirror of Lkv_mnia() for the particle-hole (c+d) diagrams: here the external line sits in the m-slot, so it is L2 and L4 whose internal (Q+L) rung contains it (i = v used for the L-part), while L1 and L3 open exactly.

Parameters

k	Multipole rank
v	Valence state: fixes kappa and the external energy, and supplies \|v> for the dressed-rung (i = v) piece
n	Excited (particle) orbital
a,b	Core (hole) orbitals
qk	Coulomb \( Q^k \) integral table (internal lines)
yk	Yk table spanning core+excited (opens the external line)
core	Core orbitals
excited	Excited orbitals
include_L4	Include core–core diagram L4
SJ	6j symbol table
Lk	Ladder table for the internal (Q+L) rung (nullptr: lowest order)

Returns: \( |L^k_{\cdot nab}\rangle \)

◆ Sigma_ladder_direct()

GMatrix MBPT::Sigma_ladder_direct	(	const DiracSpinor &	v,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		const Coulomb::QkTable &	qk,
		const Coulomb::YkTable &	yk,
		const Coulomb::LkTable *	lk,
		const Angular::SixJTable &	sjt,
		bool	include_L4 = `false`,
		double	r0 = `1.0e-4`,
		double	rmax = `30.0`,
		std::size_t	stride = `4`,
		bool	include_G = `false`
	)

Ladder correlation potential Sigma_L via the direct method.

Forms the ladder correlation potential with the external line opened exactly, rather than projected onto a basis:

\[ \Sigma_L = \sum_{amn,k} |W^k_{\cdot amn}\rangle\, \frac{1}{[k][j_v]\,(\epsilon_v+\epsilon_a-\epsilon_m-\epsilon_n)}\, \langle L^k_{mn \cdot a}| + \sum_{nab,k} |W^k_{\cdot nab}\rangle\, \frac{1}{[k][j_v]\,(\epsilon_v+\epsilon_n-\epsilon_a-\epsilon_b)}\, \langle L^k_{\cdot nab}| , \]

with the bra-side ladder vertices from Lkv_mnia() / Lkv_inab(). All internal lines are Hartree-Fock basis states; the external line is exact wherever it attaches to a bare Coulomb line (everything at lowest order in L), and is taken as |v><v| only for the dressed-rung (internal-L) attachment, which enters the energy at 4th order. Consequently Sigma_L acts correctly on Brueckner orbitals: (H + Sigma + Sigma_L)|psi_B> = e|psi_B>, rather than merely shifting by <v|Sigma_L|v>.

<v|Sigma_L|v> reproduces the ladder energy de_valence_w (up to iteration convergence of the L table). No projection basis, no Qk extension.

Parameters

v	Valence state (basis version; supplies kappa_v, en_v, and the i=v dressed-rung piece)
core	Core (hole) orbitals
excited	Excited orbitals
qk	Converged Coulomb \( Q^k \) table
yk	Yk table spanning core+excited (radial vertex functions)
lk	Converged ladder \( L^k \) table (nullptr: lowest-order L)
sjt	6j symbol table
include_L4	Include core–core diagram L4
r0,rmax,stride	Sub-grid parameters
include_G	Include the lower (g) component of Sigma_L

Returns: Sigma_L as a coordinate-space GMatrix

◆ update_Lk_mnib()

void MBPT::update_Lk_mnib	(	Coulomb::LkTable *	lk,
		const Coulomb::QkTable &	qk,
		const std::vector< DiracSpinor > &	excited,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	update_i,
		bool	include_L4,
		const Angular::SixJTable &	sjt,
		const Coulomb::LkTable *const	lk_prev,
		double	a_damp,
		bool	print
	)

Updates the ladder integral table with L(Q,Q) -> L(Q,Q+L)

Re-calculates the existing entries \( L^k_{mnib} \) of lk (those with \( i \) in update_i), using lk_prev for the internal (Q+L) rung, and stores results in lk. Designed for iterative refinement: pass the previous iteration's table as lk_prev.

Parameters

lk	Output ladder table (written in place)
qk	Coulomb \( Q^k \) integral table
excited	Excited orbitals
core	Core orbitals
update_i	Restrict re-iteration to entries whose i index is in this set (b is always core). Empty => update all. Used to converge core (update_i=core) before valence (update_i=valence).
include_L4	Include core–core diagram L4
sjt	6j symbol table
lk_prev	Ladder table from previous iteration
a_damp	Damping factor [0,1) : 0 means no damping. L -> a_dampL + (1-a_damp)L_new, where L is the existing entry of `lk` (not `lk_prev`)
print	Print progress and timing to screen

Note: Does not calculate any new integrals - assumes all already present. Just updates them (based on iterative rule: L(Q,Q) -> L(Q,Q+L))

◆ write_SigmaL()

bool MBPT::write_SigmaL	(	const std::string &	fname,
		const std::vector< SigmaLData > &	SLs,
		const Grid &	grid
	)

Writes Sigma_L (ladder) matrices to binary file.

File contains the full-grid parameters (for checking on read), followed by each Sigma_L matrix with its own sub-grid parameters and include_G flag.

Returns: True if written; false (and writes nothing) if fname is empty or 'false', or SLs is empty.

◆ read_SigmaL()

std::vector< SigmaLData > MBPT::read_SigmaL	(	const std::string &	fname,
		const std::shared_ptr< const Grid > &	grid
	)

Reads Sigma_L (ladder) matrices from binary file.

grid must match the full radial grid the matrices were calculated on. Each matrix carries its own sub-grid parameters and include_G flag (they need not match those of the base Sigma).

Returns: The Sigma_L matrices. Empty vector if fname is empty or 'false', the file doesn't exist, or on grid mismatch.

◆ k_minmax_L() [1/2]

std::pair< int, int > MBPT::k_minmax_L	(	const DiracSpinor &	a,
		const DiracSpinor &	b,
		const DiracSpinor &	c,
		const DiracSpinor &	d
	)

inline

Min and max k (multipolarity) allowed for ladder integral L^k_abcd.

Triangle rules {a,c,k} and {b,d,k} apply, plus the combined parity rule (l_a+l_b+l_c+l_d even). Unlike Q^k, there is no individual-pair parity rule linking k to the orbital parities, so k takes both parities: NOT safe to call k+=2.

Returns: {k_min, k_max}; an empty range (k_min > k_max) if no k is allowed.

◆ k_minmax_L() [2/2]

std::pair< int, int > MBPT::k_minmax_L	(	int	kap_a,
		int	kap_b,
		int	kap_c,
		int	kap_d
	)

inline

Min and max k allowed for ladder integral L^k_abcd (kappa version).

Same rules as the DiracSpinor version (see above). NOT safe to call k+=2.

◆ L3()

double MBPT::L3	(	int	k,
		const DiracSpinor &	m,
		const DiracSpinor &	n,
		const DiracSpinor &	i,
		const DiracSpinor &	j,
		const Coulomb::QkTable &	qk,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		const Angular::SixJTable &	SJ,
		const Coulomb::LkTable *const	Lk = `nullptr`,
		std::optional< double >	e_i = `{}`
	)

inline

Exchange partner of L2; equals L2 with m,n and i,j swapped.

\[ L3^k_{mnij} = L2^k_{nmji} \]

(e_i optional: used in place of i.en() in energy denominator)

◆ de_valence()

template<typename Qintegrals , typename QorLintegrals >

double MBPT::de_valence	(	const DiracSpinor &	v,
		const Qintegrals &	qk,
		const QorLintegrals &	lk,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited
	)

Second-order (or ladder) correction to the valence energy.

Computes the correlation energy shift

\[ \delta\epsilon_v = \frac{1}{[j_v]} \sum_{amn,k} \frac{Q^k_{vamn}\,(L^k_{mnva} + P^k_{nmav})} {[k]\,(\epsilon_v + \epsilon_a - \epsilon_m - \epsilon_n)} + \frac{1}{[j_v]} \sum_{nab,k} \frac{Q^k_{vnab}\,(L^k_{vnab} + P^k_{vnab})} {[k]\,(\epsilon_v + \epsilon_n - \epsilon_a - \epsilon_b)}, \]

where \( L^k \) and its exchange \( P^k \) (lk.P) both come from lk, \( a,b \) run over core, \( m,n \) over excited, and \( k \) over the Coulomb-parity range only. When lk holds plain Coulomb integrals the result is the MBPT(2) correction; when lk holds ladder integrals it is the full ladder correction.

Parameters

v	Valence orbital
qk	Coulomb \( Q^k \) integral table
lk	\( Q^k \) or ladder \( L^k \) integral table
core	Core orbitals
excited	Excited orbitals

Returns: \( \delta\epsilon_v \)

◆ de_valence_w()

template<typename Qintegrals , typename Lintegrals >

double MBPT::de_valence_w	(	const DiracSpinor &	v,
		const Qintegrals &	qk,
		const Lintegrals &	lk,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		const Angular::SixJTable *	sj = `nullptr`
	)

Ladder (or MBPT2) valence energy, antisymmetrising the FIRST integral.

Computes

\[ \delta\epsilon_v = \sum_{amn,k} \frac{W^k_{vamn}\,L^k_{mnva}}{[k]\,[j_v]\,\Delta\epsilon} + \text{(c+d)}, \]

with \( W = Q + P \) the antisymmetrised Coulomb integral (qk.W) and the ladder lk entering only through the direct integral lk.Q. Equivalent to de_valence() (which instead antisymmetrises the ladder via lk.P), since the exchange symmetry holds under the full sum. Unlike de_valence(), this works only with the Coulomb integrals in the first slot (the ladder lacks the required symmetry). No screening/eta is applied.

Parameters

v	Valence orbital
qk	Coulomb integrals supplying \( W = Q+P \) (first integral)
lk	Ladder \( L^k \) integrals (second, direct integral)
core	Core orbitals
excited	Excited orbitals
sj	Optional 6j table (speeds up the W exchange sum)

Returns: \( \delta\epsilon_v \)

◆ de_core()

template<typename Qintegrals , typename QorLintegrals >

double MBPT::de_core	(	const Qintegrals &	qk,
		const QorLintegrals &	lk,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited
	)

Second-order (or ladder) correction to the core energy.

Sums the correlation energy shift over all core orbitals. When lk holds plain Coulomb integrals the result is the MBPT(2) core correction; when lk holds ladder integrals it is the ladder correction.

Parameters

qk	Coulomb \( Q^k \) integral table
lk	\( Q^k \) or ladder \( L^k \) integral table
core	Core orbitals
excited	Excited orbitals

Returns: Total core correlation energy shift

◆ ladder()

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

Driver for the ladder-diagram calculation: the Ladder{} input block.

Calculates (and iterates to convergence) the ladder integrals Lk, writing them to the .lk file. Then constructs the ladder correlation potential, Sigma_L, for the lowest valence state of each kappa (or for each valence state, if each_valence is set), and writes these to the .sl file. The .sl file may then be read in by the Correlations block (ladder option) to include Sigma_L into the correlation potential. Runs before Correlations. Options are parsed from the Ladder{} input block.

Parameters

input	The Ladder{} input block
wf	Wavefunction: supplies the core, valence states, and basis

Warning: : Experimental, not fully tested, be careful

Note: Valence states missing from the basis are skipped (with a warning).

◆ equal() [1/2]

template<typename T >

bool MBPT::equal	(	const RadialMatrix< T > &	lhs,
		const RadialMatrix< T > &	rhs
	)

Checks if two matrices are equal (to within parts in 10^12)

◆ max_element() [1/2]

template<typename T >

double MBPT::max_element ( const RadialMatrix< T > & a )

returns maximum element (by abs)

◆ max_delta() [1/2]

template<typename T >

double MBPT::max_delta	(	const RadialMatrix< T > &	a,
		const RadialMatrix< T > &	b
	)

returns maximum difference (abs) between two matrices

◆ max_epsilon() [1/2]

template<typename T >

double MBPT::max_epsilon	(	const RadialMatrix< T > &	a,
		const RadialMatrix< T > &	b
	)

returns maximum relative difference, |(aij-bij)/(aij+bij)|, of two matrices

◆ parse_Denominators() [1/2]

std::string MBPT::parse_Denominators ( Denominators d )

Returns string representation of Denominators enum.

◆ parse_Denominators() [2/2]

Denominators MBPT::parse_Denominators ( std::string_view s )

Parses string to Denominators enum (case-insensitive); DFK if unrecognised.

◆ leg_de()

double MBPT::leg_de	(	Denominators	denominators,
		double	et_bar,
		double	et,
		double	ei_bar,
		double	ei,
		double	es,
		double	E0
	)

External-leg part of a Sigma_2 energy denominator.

For diagrams a, b, c1, and c2. Each of these diagrams has denominator (e_a - e_n) + leg_de, where a/n are the internal hole/excited states (always actual energies), and leg_de is the contribution of the two external legs, (e_target - e_intermediate). The "target" leg is the external leg whose energy slot represents the target-state energy; the "intermediate" leg is the external leg that is part of the intermediate state (the state between the two Coulomb vertices). Which energy fills each slot depends on the Denominators mode:

RS : et - ei (actual orbital energies)
Fermi : et_bar - ei_bar (Fermi-level energies, see e_bar)
Fermi0 : 0 (the legs cancel)
DFK : et_bar - ei
BW : (E0 - es) - ei. The target slot is not a single orbital energy: the whole denominator is E0 - E_intermediate, with E_int = es + ei + e_n - e_a, so the target slot becomes (E0 - es), where es is the other valence orbital in the intermediate state.

Diagram d has all four valence orbitals in its intermediate state and is handled separately (see Sigma2::S_Sigma2_d).

Parameters

denominators	Denominator mode; see MBPT::Denominators.
et_bar	Fermi-level energy (e_bar of its kappa) of the target leg.
et	Actual orbital energy of the target leg.
ei_bar	Fermi-level energy of the intermediate-state leg.
ei	Actual orbital energy of the intermediate-state leg.
es	Energy of the other valence orbital in the intermediate state (the remaining external leg); used only by BW.
E0	Total valence energy of the target CI level; used only by BW.

Returns: External-leg part of the energy denominator.

◆ e_bar()

double MBPT::e_bar	(	int	kappa_v,
		const std::vector< DiracSpinor > &	excited
	)

Returns energy of first excited state matching a given kappa.

Searches excited for the first state with the given kappa_v and returns its energy. Used to set a representative energy for a partial wave.

Parameters

kappa_v	Relativistic angular momentum quantum number.
excited	Excited (particle) states.

Returns: Energy of the matching state, if kappa is present, otherwise 0

Note: Assumes excited is sorted by energy (for each kappa); returns first (not lowest) energy; If no state with given kappa is present, returns 0. (Matrix element will be zero anyway)

◆ Sk_vwxy_SR()

bool MBPT::Sk_vwxy_SR	(	int	k,
		const DiracSpinor &	v,
		const DiracSpinor &	w,
		const DiracSpinor &	x,
		const DiracSpinor &	y
	)

Selection rule for Sk_vwxy.

Differs from the \( Q^k_{vwxy} \) selection rule due to parity.

Returns: True if \( S^k_{vwxy} \) is non-zero by selection rules.

◆ k_minmax_S() [1/2]

std::pair< int, int > MBPT::k_minmax_S	(	const DiracSpinor &	v,
		const DiracSpinor &	w,
		const DiracSpinor &	x,
		const DiracSpinor &	y
	)

Minimum and maximum k allowed by selection rules for Sk_vwxy.

Unlike \( Q^k \), \( k \) does not step by 2 for Sigma_2 matrix elements.

Returns: Pair {k_min, k_max}.

◆ k_minmax_S() [2/2]

std::pair< int, int > MBPT::k_minmax_S	(	int	twojv,
		int	twojw,
		int	twojx,
		int	twojy
	)

Overload taking 2j values directly.

◆ Sk_vwxy()

double MBPT::Sk_vwxy	(	int	k,
		const DiracSpinor &	v,
		const DiracSpinor &	w,
		const DiracSpinor &	x,
		const DiracSpinor &	y,
		const Coulomb::QkTable &	qk,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		const Angular::SixJTable &	SixJ,
		Denominators	denominators = `Denominators::DFK`,
		const std::vector< double > &	fk = `{}`,
		double	E0 = `0.0`
	)

Reduced two-body Sigma (2nd-order correlation) operator matrix element.

Computes \( S^k_{vwxy} \), the reduced matrix element of the two-body second-order correlation (Sigma_2) operator, summed over all 9 Goldstone diagrams.

\[ \begin{equation*} \begin{split} \Sigma^2_{vwxy} =& ~ \frac{g_{vnxa}\widetilde g_{awny}-g_{vnax}g_{awny}} {e_{xa}-\varepsilon_{vn}}\quad\text{(diagram 'a')}\\ &+ \frac{g_{vaxn}\widetilde g_{nway}-g_{vanx} g_{nway}} {e_{ya}-\varepsilon_{wn}} \quad\text{(diagram 'b')} \\ &-\frac{g_{vnay}g_{awxn}} {e_{ya}-\varepsilon_{vn}} -\frac{g_{vany}g_{nwxa}} {e_{xa}-\varepsilon_{wn}} \quad\text{(diagram 'c1+c2')}\\ &+ \frac{g_{vwab}g_{abxy}}{\varepsilon_{ab}-\varepsilon_{vw}} \quad\text{(diagram 'd')}. \end{split} \end{equation*} \]

The 'reduced' Sk is defined similarly to Coulomb case (Coulomb). The correlation diagrams have the same angular decomposition as the Coulomb integrals:

\[ \Sigma^2_{vwxy} = \sum_k A^k_{vwxy} S^k_{vwxy}, \]

Note: these have fewer symmetries than \( Q^k \); specifically \( S^k_{vwxy} = S^k_{wvyx} \). We call with the "Lk" symmetry (though, we should have called it "Sk"). Since each diagram is averaged with its bra-ket partner (Hermitised, see MBPT::Denominators), the bra-ket symmetry \( S^k_{vwxy} = S^k_{xyvw} \) also holds, for all denominator options.

Parameters

k	Multipolarity.
v	External spinor.
w	External spinor.
x	External spinor.
y	External spinor.
qk	Coulomb integral table (QkTable).
core	Core (hole) states for internal lines.
excited	Excited (particle) states (internal lines).
SixJ	Precomputed 6-j symbol table.
denominators	Energy denominator convention: see MBPT::Denominators.
fk	Screening factors; fk[k] scales the k-th Coulomb line. Missing (or empty) implies 1.0 (no screening).
E0	Total valence energy of the target CI level; used only by Denominators::BW.

Returns: \( S^k_{vwxy} \).

◆ calculate_Sk()

Coulomb::LkTable MBPT::calculate_Sk	(	const std::string &	filename,
		const std::vector< DiracSpinor > &	external,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		const Coulomb::QkTable &	qk,
		int	max_k,
		bool	exclude_wrong_parity_box,
		Denominators	denominators,
		bool	no_new_integrals = `false`,
		const std::vector< double > &	fk = `{}`,
		double	E0 = `0.0`
	)

Calculates (or reads in) a table of two-body Sigma_2 matrix elements.

Computes \( S^k_{vwxy} \) for all relevant combinations of states in external, using the provided core and excited bases and Coulomb table. Results are written to / read from filename (empty string or "false" disables I/O).

Parameters

filename	File to read/write the table (blank or "false" to not read/write).
external	Basis states for external legs (all ME between these are computed).
core	Core (hole) states for internal summations.
excited	Excited (particle) states for internal summations.
qk	Precomputed Coulomb integral table (QkTable).
max_k	Maximum multipolarity to include (<= 0: no limit).
exclude_wrong_parity_box	If true, excludes box diagrams with "wrong" parity.
denominators	DFK, BW, RS, Fermi, Fermi0: see MBPT::Denominators
no_new_integrals	If true, only reads existing integrals; no new computation.
fk	Screening factors; fk[k] scales the k-th Coulomb line.
E0	Target-level valence energy; used only by Denominators::BW.

Returns: LkTable containing all computed \( S^k_{vwxy} \) matrix elements.

Note: no_new_integrals - if we know all required integrals are already in the file to be read in, saves time. Otherwise, ampsci will check if any new integrals are required. This checking can take a while, particularly for large basis.

◆ average_hk()

std::vector< double > MBPT::average_hk	(	const Coulomb::LkTable &	Sk,
		const Coulomb::QkTable &	qk,
		const std::vector< DiracSpinor > &	external,
		int	max_k = `-1`
	)

Average Sigma_2 correction ratios, h_k, for each multipole k.

h_k = <S^k/Q^k>, averaged over all stored S^k integrals with |Q^k| above a small cut-off. Each distinct integral is counted once (the 4-fold symmetry of S^k is accounted for), so integrals are weighted equally. Used to extrapolate Sigma_2 to diagrams outside the tabulated set: S^k ~ h_k Q^k. Entries with no data are 0.0 (no correction).

Parameters

Sk	Table of Sigma_2 integrals (see calculate_Sk).
qk	Coulomb integral table.
external	States for external legs (typically the cis2 basis).
max_k	Maximum multipole; if negative, determined from basis.

Returns: Vector of average correction ratios, indexed by k.

◆ Sigma_vw()

template<class CoulombIntegral >

double MBPT::Sigma_vw	(	const DiracSpinor &	v,
		const DiracSpinor &	w,
		const CoulombIntegral &	qk,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		int	max_l_internal = `99`,
		std::optional< double >	ev = `std::nullopt`
	)

Matrix element of the 1-body Sigma (2nd-order correlation) operator.

Matrix element of 1-body Sigma (2nd-order correlation) operator; de_v = <v|Sigma|v>. qk (CoulombIntegral) may be YkTable or QkTable.

Computes \( \langle v | \Sigma(E) | w \rangle \) by summing over internal core and excited states using the provided Coulomb integral table.

The energy at which Sigma is evaluated:

If ev is given, it is used directly.
Otherwise \( E = \frac{1}{2}(\varepsilon_v + \varepsilon_w) \) is used.

max_l_internal truncates the angular momentum of internal lines; intended for convergence tests only.

Parameters

v	External bra spinor.
w	External ket spinor.
qk	Coulomb integral table (YkTable or QkTable).
core	Core (hole) states.
excited	Excited (particle) states.
max_l_internal	Maximum \( l \) for internal lines (default 99).
ev	Optional energy at which Sigma is evaluated.

Returns: \( \langle v | \Sigma(E) | w \rangle \).

◆ Sigma_vw_direct_exchange()

template<class CoulombIntegral >

std::pair< double, double > MBPT::Sigma_vw_direct_exchange	(	const DiracSpinor &	v,
		const DiracSpinor &	w,
		const CoulombIntegral &	qk,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		int	max_l_internal = `99`,
		std::optional< double >	ev = `std::nullopt`
	)

Direct and exchange parts of Sigma_vw(), returned separately.

Direct and exchange parts of the 1-body Sigma (2nd-order correlation) matrix element, returned separately as {direct, exchange}; <v|Sigma|w> = direct + exchange. qk may be YkTable or QkTable.

Direct and exchange parts of \( \langle v | \Sigma(E) | w \rangle \), returned separately as {direct, exchange}. As Sigma_vw(), but with the direct (QQ) and exchange (QP) contributions accumulated separately; Sigma_vw() returns their sum.

Returns: Pair {direct, exchange}.

◆ dSigma_dE_vw()

template<class CoulombIntegral >

double MBPT::dSigma_dE_vw	(	const DiracSpinor &	v,
		const DiracSpinor &	w,
		const CoulombIntegral &	qk,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		double	ev,
		int	max_l_internal = `99`,
		double	delta = `0.01`
	)

Energy derivative of the one-body correlation correction.

\( d\langle v|\Sigma(E)|w\rangle/dE \), evaluated at E = ev. Central finite difference of Sigma_vw(), with step delta.

Parameters

v,w	External spinors.
qk	Coulomb integral table (YkTable or QkTable).
core	Core (hole) states.
excited	Excited (particle) states.
ev	Energy at which the derivative is evaluated.
max_l_internal	Maximum \( l \) for internal lines (default 99).
delta	Finite-difference step, in au.

Returns: \( d\langle v|\Sigma(E)|w\rangle/dE \).

◆ equal() [2/2]

template<typename T >

bool MBPT::equal	(	const SpinorMatrix< T > &	lhs,
		const SpinorMatrix< T > &	rhs
	)

Checks if two matrices are equal (to within parts in 10^12)

◆ max_element() [2/2]

template<typename T >

double MBPT::max_element ( const SpinorMatrix< T > & a )

returns maximum element (by abs)

◆ max_delta() [2/2]

template<typename T >

double MBPT::max_delta	(	const SpinorMatrix< T > &	a,
		const SpinorMatrix< T > &	b
	)

Returns maximum difference (abs) between two matrices.

A matrix without g parts is treated as having zero g parts.

◆ max_epsilon() [2/2]

template<typename T >

double MBPT::max_epsilon	(	const SpinorMatrix< T > &	a,
		const SpinorMatrix< T > &	b
	)

Returns maximum relative difference (abs) between two matrices.

Maximum over all elements of \( |(a_{ij}-b_{ij})/(a_{ij}+b_{ij})| \). A matrix without g parts is treated as having zero g parts.

Detailed Description

Namespaces

Classes

Typedefs

Enumerations

Functions

Variables

Class Documentation

◆ MBPT::FeynmanOptions

◆ MBPT::rgrid_params

◆ MBPT::SigmaData

◆ MBPT::SigmaLData

Enumeration Type Documentation

◆ SigmaMethod

◆ HoleParticle

◆ Screening

◆ GreenStates

◆ SigmaLMethod

◆ Denominators

Function Documentation

◆ fit_tail_coefficient()

◆ add_tail()

◆ best_omre()

◆ parseSigmaLMethod() [1/2]

◆ parseSigmaLMethod() [2/2]

◆ Lkmnij()

◆ L1()

◆ L4()

◆ L2()

◆ fill_Lk_mnib()

◆ Sigma_ladder()

◆ Lkv_mnia()

◆ Lkv_inab()

◆ Sigma_ladder_direct()

◆ update_Lk_mnib()

◆ write_SigmaL()

◆ read_SigmaL()

◆ k_minmax_L() [1/2]

◆ k_minmax_L() [2/2]

◆ L3()

◆ de_valence()

◆ de_valence_w()

◆ de_core()

◆ ladder()

◆ equal() [1/2]

◆ max_element() [1/2]

◆ max_delta() [1/2]

◆ max_epsilon() [1/2]

◆ parse_Denominators() [1/2]

◆ parse_Denominators() [2/2]

◆ leg_de()

◆ e_bar()

◆ Sk_vwxy_SR()

◆ k_minmax_S() [1/2]

◆ k_minmax_S() [2/2]

◆ Sk_vwxy()

◆ calculate_Sk()

◆ average_hk()

◆ Sigma_vw()

◆ Sigma_vw_direct_exchange()

◆ dSigma_dE_vw()

◆ equal() [2/2]

◆ max_element() [2/2]

◆ max_delta() [2/2]

◆ max_epsilon() [2/2]