Detailed Description

Many-body perturbation theory.

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

Classes
class	Feynman
	Class to construct Feynman diagrams, Green's functions and polarisation op. More...

class	Goldstone
	Class to construct Feynman diagrams, Green's functions and polarisation op. More...

class	RadialMatrix

struct	SigmaLData
	Ladder correlation potential matrix, Sigma_L, for a single (kappa, n) state, evaluated at energy en. More...

class	SpinorMatrix

class	StructureRad
	Calculates Structure Radiation + 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 }

enum class	HoleParticle { exclude , include , include_k0 }
	Options for including hole-particle interaction. include mean all k; include_k0 means k=0 term only. More...

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

enum class	GreenStates { both , core , excited }
	Which states to include in 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	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.

SigmaLMethod	parseSigmaLMethod (const std::string &method)
	Converts string (name) to SigmaLMethod enum (case-insensitive); warns and defaults to ladder if unknown.

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), by projection (or, if `projection` is empty, by the L/Q ratio method).

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 the ladder integral L^k_mnia over the 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 the ladder integral L^k_inab over the external index i (the m-slot).

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 (open external line) 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)
	Returns 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.

std::pair< int, int >	k_minmax_L (int kap_a, int kap_b, int kap_c, int kap_d)
	Returns min and max k (multipolarity) allowed for ladder integral L^k_abcd (kappa version) - see above. NOT safe to call k+=2.

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 matrix's 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 matrixs

template<typename T >
double	max_epsilon (const RadialMatrix< T > &a, const RadialMatrix< T > &b)
	returns maximum relative diference [aij-bij/(aij+bij)] (abs) between 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); returns 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 (diagrams a, b, c1, c2).

std::pair< std::vector< DiracSpinor >, std::vector< DiracSpinor > >	split_basis (const std::vector< DiracSpinor > &basis, double E_Fermi, int min_n_core=1, int max_n_excited=999)
	Splits the basis into the core (holes) and excited states.

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 \( S^k_{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 \( S^k_{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 \( \langle v \| \Sigma(E) \| w \rangle \), returned separately as {direct, exchange}.

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, \( d\langle v\|\Sigma(E)\|w\rangle/dE \), evaluated at E = `ev`.

template<typename T >
bool	equal (const SpinorMatrix< T > &lhs, const SpinorMatrix< T > &rhs)
	Checks if two matrix's 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 matrixs

template<typename T >
double	max_epsilon (const SpinorMatrix< T > &a, const SpinorMatrix< T > &b)
	returns maximum relative diference [aij-bij/(aij+bij)] (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::SigmaLData

struct MBPT::SigmaLData

Class Members
int	kappa
int	n
double	en
GMatrix	SL

Enumeration Type Documentation

◆ HoleParticle

enum class MBPT::HoleParticle

strong

Options for including hole-particle interaction. include mean all k; include_k0 means 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 Green's function:

◆ 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 orbtials.

Function Documentation

◆ 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 = e_lowest_excited - e_core_max, 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). Returns the midpoint of the largest pole-free gap; one omre serves all valence states. Set print=true to show Delta, the in-window poles, and the chosen omre.

◆ parseSigmaLMethod() [1/2]

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

Converts string (name) to SigmaLMethod enum (case-insensitive); warns and defaults to ladder 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.

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} \) 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 Qk info 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). 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 the ladder integral L^k_mnia over the 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).

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 the ladder integral L^k_inab over the external index i (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.

◆ 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 (open external line) 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)

Iterates over all combinations of excited pairs \( (m,n) \) and orbitals in i_orbs, computing \( L^k_{mnib} \) and storing results in lk. Designed for iterative refinement: pass the previous iteration's table as lk_prev.

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))

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
print	Print Qk info to screen

◆ 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 false (and writes nothing) if fname is empty or 'false'.

◆ 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.

Returns empty vector if file doesn't exist or on grid mismatch: 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).

◆ 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

Returns 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.

◆ 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

Returns min and max k (multipolarity) allowed for ladder integral L^k_abcd (kappa 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 = \sum_{mnc} \frac{Q^k_{vmcn}\,L^k_{mncv}}{\epsilon_v + \epsilon_c - \epsilon_m - \epsilon_n} \]

(schematic). 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 each valence state, and writes these to the .sl file. The .sl file may then be read in by the Correlations block (ladder_file option) to include Sigma_L into the correlation potential. Runs before Correlations. Options are parsed from the Ladder{} input block.

◆ equal() [1/2]

template<typename T >

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

Checks if two matrix's 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 matrixs

◆ max_epsilon() [1/2]

template<typename T >

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

returns maximum relative diference [aij-bij/(aij+bij)] (abs) between 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); returns 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 (diagrams a, b, c1, 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.

◆ split_basis()

std::pair< std::vector< DiracSpinor >, std::vector< DiracSpinor > > MBPT::split_basis	(	const std::vector< DiracSpinor > &	basis,
		double	E_Fermi,
		int	min_n_core = `1`,
		int	max_n_excited = `999`
	)

Splits the basis into the core (holes) and excited states.

States with energy below E_Fermi are considered core/holes. Only core states with \( n \geq \) min_n_core, and excited states with \( n \leq \) max_n_excited are kept.

Note: Negative energy states not dealt with! Assumed not to be present in basis. Fix?; Replace all instances with DiracSpinor::split_by_energy

Parameters

basis	Full set of single-particle basis states.
E_Fermi	Energy threshold separating core from excited states.
min_n_core	Minimum principal quantum number for core states.
max_n_excited	Maximum principal quantum number for excited states.

Returns: Pair {core, excited} of DiracSpinor vectors.

◆ 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.

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)

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

◆ 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 \( S^k_{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 \( S^k_{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 disables I/O).

Parameters

filename	File to read/write the table. (blank for "false" to not 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.
exclude_wrong_parity_box	If true, excludes box diagrams with "wrong" parity.
denominators	DFK, RS, Fermi, Fermi0: see MBPT::Denominators
no_new_integrals	If true, only reads existing intergals; no new computation.
fk	Screening factors; fk[k] scales the k-th Coulomb line.
E0	Target-level valence energy; used only by Denominators::BW.

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 a requred. This checking can take a while, particularly for large basis.

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

◆ 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 \( \langle v | \Sigma(E) | w \rangle \), returned separately as {direct, exchange}.

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.

As Sigma_vw(), but with the direct (QQ) and exchange (QP) contributions accumulated separately; Sigma_vw() returns their sum.

◆ 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 matrix's 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 matrixs

◆ max_epsilon() [2/2]

template<typename T >

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

returns maximum relative diference [aij-bij/(aij+bij)] (abs) between two matrices

Detailed Description

Namespaces

Classes

Typedefs

Enumerations

Functions

Variables

Class Documentation

◆ MBPT::SigmaLData

Enumeration Type Documentation

◆ HoleParticle

◆ Screening

◆ GreenStates

◆ SigmaLMethod

◆ Denominators

Function Documentation

◆ 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()

◆ split_basis()

◆ 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]