ok

Structure Radiation and Normalisation of states, using diagram method.

Two main functions: SR(), norm(), which return the reduced matrix elements for structure radiation and normalisation, respectively.

Structure radiation is sum of srTB()+srC(); norm() gives normalisation of states.

Typically, one should use splines for the 'legs' (outer states) of diagrams. However, code is written such that {w,v} (the valence states) always appear as the first state in the Q integrals; thus those states are always directly integrated (other states may use the existing y^k integrals, which were calculated using the spline states), unless using QkTable (see below).

For using spline states as legs: This means you must typically ensure basis is large enough to make relevant valence states "physical"

The user should check if the zeroth order <w|t|v> matches closely between valence/spline states. If not, basis/cavity is too small and SR+N probably meaningless

Option to use QkTable to calculate SR+N. Otherwise, calculates Coulomb integrals on-the-fly. If QkTable is used, leads to ~10x speed-up, but uses large amount of memory. Also, using QkTable means spline states are used for "legs" of diagrams.

norm() works differently: it directly takes in the operator to calculate norm. This is because norm depends only on Coulomb operators. In some sense, would make sense to keep norm in a separate class. However, they are typically calculated together, and use the same Coulomb integrals.

Formulas are presented in (for example): W. R. Johnson, Z. W. Liu, and J. Sapirstein, At. Data Nucl. Data Tables 64, 279 (1996), https://www.sciencedirect.com/science/article/pii/S0092640X96900248

Specific form of formulas used here are presented in: Phys. Rev. A 107, 052812 (2023)

Note: For Structure Radiation, must call StructureRad::solve_core() first.

#include <StructureRad.hpp>

Public Member Functions
	StructureRad (const std::vector< DiracSpinor > &basis, double en_core, std::pair< int, int > nminmax={0, 999}, const std::string &Qk_fname="", int k_cut=99, const std::vector< double > &fk={}, const std::vector< double > &etak={}, bool verbose=true)
	Construct a Structure Radiation (MBPT) object.

const std::vector< DiracSpinor > &	core () const
	const reference to employed subset of core orbitals

const std::vector< DiracSpinor > &	excited () const
	const reference to employed subset of excited orbitals

const std::vector< DiracSpinor > &	basis () const
	const reference to employed subset of core+excited orbitals

const Coulomb::meTable< double > &	me_Table () const
	const reference to underlying matrix element table. May be empty

const Coulomb::YkTable &	Yk () const
	const reference to Yk table (always formed by the constructor)

const std::optional< Coulomb::QkTable > &	Qk () const
	const reference to underlying Qk table. Optional; may not be set

void	solve_core (const DiracOperator::TensorOperator const h, const ExternalField::CorePolarisation const dV=nullptr, double omega=0.0)
	Prepares the operator matrix elements required for SR and BO.

double	SR (const DiracSpinor &w, const DiracSpinor &v, double omega=0.0) const
	Reduced matrix element of structure radiation, <w\|\|h_sr\|\|v>.

double	norm (const DiracSpinor &w, const DiracSpinor &v, const DiracOperator::TensorOperator const h, const ExternalField::CorePolarisation const dV=nullptr, bool conj=false) const
	Returns Normalisation of states, reduced ME: <w\|\|h\|\|v>_norm.

double	f_norm (const DiracSpinor &v) const
	Normalisation factor; defined: <w\|\|h\|\|v>_norm = <w\|\|h\|\|v>(f_v + f_w)

double	norm_factor (const DiracSpinor &w, const DiracSpinor &v) const
	Normalisation factor for pair: <w\|\|h\|\|v>_norm = norm_factor(w,v) <w\|\|h\|\|v>

double	srn (const DiracSpinor &w, const DiracSpinor &v, const DiracOperator::TensorOperator const h, const ExternalField::CorePolarisation const dV=nullptr, double omega=0.0) const
	Returns sum of SR+Norm diagrams, reduced ME: <w\|\|T+B+C+N\|\|v>.

double	f_root_scr (int k) const
	Effective screening factor for each Coulomb line, sqrt(fk); 1 if not set.

double	eta_hp (int k) const
	Effective hole-particle factor for polarisation loops; 1 if not set.

double	BO (const DiracSpinor &w, const DiracSpinor &v, double fw=1.0, double fv=1.0) const
	Returns Brueckner orbital contribution to the reduced ME.

double	Sigma_vw (const DiracSpinor &v, const DiracSpinor &w) const
	Second-order <v\|Sigma\|w>, at E = (e_v + e_w)/2 (fk and etak not applied)

Coulomb::meTable< double >	srn_table (const DiracOperator::TensorOperator const h, const std::vector< DiracSpinor > &as, const std::vector< DiracSpinor > &tbs={}, double omega=0.0, const ExternalField::CorePolarisation const dV=nullptr) const
	Constructs an me table of srn() for each pair {a,b}.

Constructor & Destructor Documentation

◆ StructureRad()

MBPT::StructureRad::StructureRad	(	const std::vector< DiracSpinor > &	basis,
		double	en_core,
		std::pair< int, int >	nminmax = `{0, 999}`,
		const std::string &	Qk_fname = `""`,
		int	k_cut = `99`,
		const std::vector< double > &	fk = `{}`,
		const std::vector< double > &	etak = `{}`,
		bool	verbose = `true`
	)

Construct a Structure Radiation (MBPT) object.

en_core defines the core boundary: states with e < en_core are treated as core states, while states with e > en_core are treated as excited states. Typically en_core is half-way between max(e_core) and min(e_valence), i.e., wf.FermiLevel().

nminmax is the pair {nmin, nmax}. Only core states with n >= nmin are used, and only excited states with n <= nmax are included in the summations. Negative-energy states (if in the basis) are always included as excited states.

Qk_fname is the filename for the QkTable. If provided, the QkTable will be read from / written to this file and then used to calculate SR. This leads to ~10x speedup in the calculation at the cost of significantly higher memory usage.

fk and etak are effective screening factors for screening and hole-particle. These should only be used as a test.

Parameters

basis	Basis orbitals (typically splines); split into core and excited states by `en_core`
en_core	Core/excited energy boundary
nminmax	{nmin, nmax}: n range of core/excited states included
Qk_fname	QkTable file; if empty, no QkTable (Yk table used)
k_cut	Maximum multipolarity k of the Coulomb integrals (see note)
fk	Effective screening factors, for each k (test only)
etak	Effective hole-particle factors, for each k (test only)
verbose	If true, prints progress to screen

Note: k_cut applies to every Coulomb line, including those inside the exchange integrals (all Q^k, P^k, W^k with k > k_cut are zero), whether or not a QkTable is used.

Member Function Documentation

◆ core()

const std::vector< DiracSpinor > & MBPT::StructureRad::core ( ) const

inline

const reference to employed subset of core orbitals

◆ excited()

const std::vector< DiracSpinor > & MBPT::StructureRad::excited ( ) const

inline

const reference to employed subset of excited orbitals

◆ basis()

const std::vector< DiracSpinor > & MBPT::StructureRad::basis ( ) const

inline

const reference to employed subset of core+excited orbitals

◆ me_Table()

const Coulomb::meTable< double > & MBPT::StructureRad::me_Table ( ) const

inline

const reference to underlying matrix element table. May be empty

◆ Yk()

const Coulomb::YkTable & MBPT::StructureRad::Yk ( ) const

inline

const reference to Yk table (always formed by the constructor)

◆ Qk()

const std::optional< Coulomb::QkTable > & MBPT::StructureRad::Qk ( ) const

inline

const reference to underlying Qk table. Optional; may not be set

◆ solve_core()

void MBPT::StructureRad::solve_core	(	const DiracOperator::TensorOperator *const	h,
		const ExternalField::CorePolarisation *const	dV = `nullptr`,
		double	omega = `0.0`
	)

Prepares the operator matrix elements required for SR and BO.

Pre-computes the full set of <i||T||j> matrix elements
Does for all core-core, core-excited, excited-excited orbitals
Optionally includes core polarisation; dV may be null
h must not be null
For frequency-dependent operators, this assumes operator has already been updated
Assumes Core Polarisation has already been solved

Parameters

h	Operator (must not be null)
dV	Core polarisation (RPA); included in the matrix elements if not null
omega	Transition frequency: selects the dV component on the internal lines (emission if negative); see CorePolarisation::dV

Note: This must be called before structure radiation (or BO) can be calculated

◆ SR()

double MBPT::StructureRad::SR	(	const DiracSpinor &	w,
		const DiracSpinor &	v,
		double	omega = `0.0`
	)		const

inline

Reduced matrix element of structure radiation, <w||h_sr||v>.

Sum of the top, bottom, and centre diagrams, for the operator (and dV) given to solve_core().

Parameters

w	Bra state
v	Ket state
omega	Operator frequency; enters the energy denominators of the top and bottom diagrams

Note: Must call solve_core() first; aborts otherwise.

Warning: Returns zero if <w||h||v> is not in the solve_core() table: i.e., if zero by selection rules, or if w or v (matched by n and kappa) is not in the basis.

◆ norm()

double MBPT::StructureRad::norm	(	const DiracSpinor &	w,
		const DiracSpinor &	v,
		const DiracOperator::TensorOperator *const	h,
		const ExternalField::CorePolarisation *const	dV = `nullptr`,
		bool	conj = `false`
	)		const

Returns Normalisation of states, reduced ME: <w||h||v>_norm.

\[ \redmatel{w}{h}{v}_{\rm norm} = (f_w + f_v) \redmatel{w}{h + \delta V}{v}, \]

with f from f_norm(). Does not require solve_core().

Parameters

w	Bra state
v	Ket state
h	Operator (must not be null)
dV	Core polarisation (RPA); included if not null
conj	The dV component: emission (true) for a negative frequency; see CorePolarisation::dV

◆ f_norm()

double MBPT::StructureRad::f_norm ( const DiracSpinor & v ) const

inline

Normalisation factor; defined: <w||h||v>_norm = <w||h||v>(f_v + f_w)

◆ norm_factor()

double MBPT::StructureRad::norm_factor	(	const DiracSpinor &	w,
		const DiracSpinor &	v
	)		const

Normalisation factor for pair: <w||h||v>_norm = norm_factor(w,v) <w||h||v>

◆ srn()

double MBPT::StructureRad::srn	(	const DiracSpinor &	w,
		const DiracSpinor &	v,
		const DiracOperator::TensorOperator *const	h,
		const ExternalField::CorePolarisation *const	dV = `nullptr`,
		double	omega = `0.0`
	)		const

inline

Returns sum of SR+Norm diagrams, reduced ME: <w||T+B+C+N||v>.

SR() + norm().

Parameters

w	Bra state
v	Ket state
h	Operator, for norm() only (must not be null)
dV	Core polarisation (RPA), for norm() only; may be null
omega	Operator frequency: for the SR() denominators, and the sign selects the dV component in norm()

Note: Must call solve_core() first (aborts otherwise). h and dV used for norm() only

◆ f_root_scr()

double MBPT::StructureRad::f_root_scr ( int k ) const

inline

Effective screening factor for each Coulomb line, sqrt(fk); 1 if not set.

◆ eta_hp()

double MBPT::StructureRad::eta_hp ( int k ) const

inline

Effective hole-particle factor for polarisation loops; 1 if not set.

◆ BO()

double MBPT::StructureRad::BO	(	const DiracSpinor &	w,
		const DiracSpinor &	v,
		double	fw = `1.0`,
		double	fv = `1.0`
	)		const

Returns Brueckner orbital contribution to the reduced ME.

From 10.1006/adnd.1996.0024:

\[ f_v \redmatel{w}{h}{\delta v} + f_w \redmatel{\delta w}{h}{v}, \qquad \ket{\delta v} = \sum_{i \neq v} \frac{\ket{i}\matel{i}{\Sigma}{v}}{\varepsilon_v - \varepsilon_i}, \]

where \( \Sigma \) is the second-order correlation potential at energy \( \varepsilon_v \) (with fk and etak, if given), and the sum runs over core and excited basis states with \( \kappa_i = \kappa_v \). The \( \redmatel{w}{h}{i} \) are taken from the solve_core() table (so include dV, if given there).

Parameters

w	Bra state
v	Ket state
fw	Factor multiplying the \( \redmatel{\delta w}{h}{v} \) term
fv	Factor multiplying the \( \redmatel{w}{h}{\delta v} \) term

Note: Must call solve_core() first; aborts otherwise.

Warning: Terms are zero if w or v (matched by n and kappa) is not in the basis.

◆ Sigma_vw()

double MBPT::StructureRad::Sigma_vw	(	const DiracSpinor &	v,
		const DiracSpinor &	w
	)		const

Second-order <v|Sigma|w>, at E = (e_v + e_w)/2 (fk and etak not applied)

◆ srn_table()

Coulomb::meTable< double > MBPT::StructureRad::srn_table	(	const DiracOperator::TensorOperator *const	h,
		const std::vector< DiracSpinor > &	as,
		const std::vector< DiracSpinor > &	tbs = `{}`,
		double	omega = `0.0`,
		const ExternalField::CorePolarisation *const	dV = `nullptr`
	)		const

Constructs an me table of srn() for each pair {a,b}.

For each a in as and b in tbs (or in as, if tbs is empty) with non-zero <a||h||b>. The {b,a} element is also stored, using the symmetry of h.

Parameters

h	Operator (must not be null): selection rules, symmetry, and norm()
as	First set of states
tbs	Second set of states; if empty, `as` is used
omega	Operator frequency, for SR() only
dV	Core polarisation (RPA), for norm() only; may be null

Returns: Table of srn(a, b, h, dV, omega).

Note: Must call solve_core() first.

The documentation for this class was generated from the following files:

MBPT/StructureRad.hpp
MBPT/StructureRad.cpp

ok

Public Member Functions

Constructor & Destructor Documentation

◆ StructureRad()

Member Function Documentation

◆ core()

◆ excited()

◆ basis()

◆ me_Table()

◆ Yk()

◆ Qk()

◆ solve_core()

◆ SR()

◆ norm()

◆ f_norm()

◆ norm_factor()

◆ srn()

◆ f_root_scr()

◆ eta_hp()

◆ BO()

◆ Sigma_vw()

◆ srn_table()