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

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	Denominators { RS , Fermi , Fermi0 }
	Type of energy demoninators: RS, Fermi, Fermi0. More...

Functions
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, const std::vector< double > &fk={})
	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, const std::vector< double > &fk={})
	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, const std::vector< double > &fk={})
	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, const std::vector< double > &fk={})
	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, bool print=true, const std::vector< double > &fk={})
	Fills the ladder integral table for all required index combinations.

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 > &i_orbs, bool include_L4, const Angular::SixJTable &sjt, const Coulomb::LkTable const lk_prev, double a_damp, bool print, const std::vector< double > &fk)
	Updates the ladder integral table with L(Q,Q) -> L(Q,Q+L)

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, const std::vector< double > &fk={})
	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, const std::vector< double > &fk={}, const std::vector< double > &etak={})
	Second-order (or ladder) correction to the valence energy.

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.

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 Fermi0 if unrecognised.

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::Fermi0)
	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)
	Calculates (or reads in) a table of two-body Sigma_2 matrix elements.

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<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
const auto	vroot = [](auto x) { return std::sqrt(x); }

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:

◆ Denominators

enum class MBPT::Denominators

strong

Type of energy demoninators: RS, Fermi, Fermi0.

RS : Use energy denominator ascociated with actual orbital for external legs. Symmeterised so that diagram is symmetric. May be danger of accidental enhancement.
Fermi : Use energy from the "Fermi" level (i.e., lowest state for given kappa in excited spectrum).
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.

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

Function Documentation

◆ 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`,
		const std::vector< double > &	fk = `{}`
	)

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)
fk	Optional screening factors \( f_k \)

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

◆ 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`,
		const std::vector< double > &	fk = `{}`
	)

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)
fk	Optional screening factors \( f_k \)

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`,
		const std::vector< double > &	fk = `{}`
	)

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. Typically small and omitted by default; 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)
fk	Optional screening factors \( f_k \)

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`,
		const std::vector< double > &	fk = `{}`
	)

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)
fk	Optional screening factors \( f_k \)

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,
		bool	print = `true`,
		const std::vector< double > &	fk = `{}`
	)

Fills the ladder integral table for all required 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.=

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
print	Print Qk info to screen
fk	Optional screening factors \( f_k \)

◆ 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 > &	i_orbs,
		bool	include_L4,
		const Angular::SixJTable &	sjt,
		const Coulomb::LkTable *const	lk_prev,
		double	a_damp,
		bool	print,
		const std::vector< double > &	fk
	)

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
i_orbs	Orbitals for the \( i \) index
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
fk	Optional screening factors \( f_k \)

◆ 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`,
		const std::vector< double > &	fk = `{}`
	)

inline

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

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

◆ 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,
		const std::vector< double > &	fk = `{}`,
		const std::vector< double > &	etak = `{}`
	)

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
fk	Optional screening factors \( f_k \)
etak	Optional enhancement factors \( \eta_k \)

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

◆ 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 Fermi0 if unrecognised.

◆ 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::Fermi0`
	)

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_{xyvw} \). We call with the "Lk" symmetry (though, we should have called it "Sk")

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 (RS or BW) not implemented correctly .

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

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	RS, Fermi, Fermi0: see MBPT::Denominators
no_new_integrals	If true, only reads existing intergals; no new computation.

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.

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

◆ 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

Enumeration Type Documentation

◆ HoleParticle

◆ Screening

◆ GreenStates

◆ Denominators

Function Documentation

◆ Lkmnij()

◆ L1()

◆ L4()

◆ L2()

◆ fill_Lk_mnib()

◆ update_Lk_mnib()

◆ L3()

◆ de_valence()

◆ de_core()

◆ equal() [1/2]

◆ max_element() [1/2]

◆ max_delta() [1/2]

◆ max_epsilon() [1/2]

◆ parse_Denominators() [1/2]

◆ parse_Denominators() [2/2]

◆ split_basis()

◆ e_bar()

◆ Sk_vwxy_SR()

◆ k_minmax_S() [1/2]

◆ k_minmax_S() [2/2]

◆ Sk_vwxy()

◆ calculate_Sk()

◆ Sigma_vw()

◆ equal() [2/2]

◆ max_element() [2/2]

◆ max_delta() [2/2]

◆ max_epsilon() [2/2]