Amplitudes Namespace Reference

Detailed Description

Physical amplitudes and observables (matrix elements, second-order amplitudes); testable functions, callable from any module.

Classes
struct	MEdata
	Result of a single matrix element calculation. More...
struct	MEoptions
	Options for the matrix_elements() list driver. More...
struct	SRNdata
	Result of a matrix element calculation with second-order MBPT corrections: structure radiation, normalisation, Brueckner orbital. More...
struct	SRNoptions
	Options for sr_matrix_elements() More...

Enumerations
enum class	Frequency { transition , fixed }
	Which frequency the operator, or the RPA, is evaluated at. More...

Functions
MEdata	matrix_element (const DiracSpinor &a, const DiracSpinor &b, const DiracOperator::TensorOperator *h, const DiracOperator::TensorOperator *h_minus=nullptr, const ExternalField::CorePolarisation *dV=nullptr, DiracOperator::MatrixElementType type=DiracOperator::MatrixElementType::Reduced, double omega=0.0)
	Single matrix element of h between states a and b, with optional RPA.
void	set_operator_frequency (DiracOperator::TensorOperator *h, DiracOperator::TensorOperator *h_minus, double omega)
	Sets the \( t_\pm \) operator pair to the frequency w.
std::vector< MEdata >	matrix_elements (const std::vector< DiracSpinor > &a_orbs, const std::vector< DiracSpinor > &b_orbs, DiracOperator::TensorOperator *h, DiracOperator::TensorOperator *h_minus, ExternalField::CorePolarisation *dV, const MEoptions &options, std::ostream &outstream=std::cout)
	Matrix elements of h for all allowed pairs from two lists of orbitals, with optional RPA; owns all frequency updates and RPA solves.
Coulomb::meTable< double >	me_table (const std::vector< DiracSpinor > &a_orbs, const std::vector< DiracSpinor > &b_orbs, const DiracOperator::TensorOperator *h, const ExternalField::CorePolarisation *dV=nullptr)
	Builds a lookup table of reduced matrix elements <a||h||b>.
Coulomb::meTable< double >	me_table (const std::vector< DiracSpinor > &a_orbs, const std::vector< DiracSpinor > &b_orbs, const DiracOperator::TensorOperator *h, const ExternalField::CorePolarisation *dV, const MBPT::StructureRad *srn, double omega, int sr_n_max=999, bool sr_norm=true)
	Builds a table of reduced matrix elements, including structure radiation and (optionally) the normalisation of states.
std::vector< SRNdata >	sr_matrix_elements (const std::vector< DiracSpinor > &orbs, DiracOperator::TensorOperator *h, MBPT::StructureRad *sr, ExternalField::CorePolarisation *dV, const SRNoptions &options, std::ostream &outstream=std::cout)
	Matrix elements with structure radiation, normalisation, and Brueckner orbital corrections, for all pairs from a list of orbitals.
std::vector< MEdata >	matrix_elements (const std::vector< DiracSpinor > &orbs, DiracOperator::TensorOperator *h, DiracOperator::TensorOperator *h_minus, ExternalField::CorePolarisation *dV, const MEoptions &options, std::ostream &outstream=std::cout)
	Matrix elements of h for all pairs from a single list of orbitals.
Coulomb::meTable< double >	me_table (const std::vector< DiracSpinor > &a_orbs, const DiracOperator::TensorOperator *h, const ExternalField::CorePolarisation *dV=nullptr)
	Builds a matrix element table for a single set of orbitals.
Coulomb::meTable< double >	me_table (const std::vector< DiracSpinor > &a_orbs, const DiracOperator::TensorOperator *h, const ExternalField::CorePolarisation *dV, const MBPT::StructureRad *srn, double omega, int sr_n_max=999, bool sr_norm=true)
	Builds a SR+N matrix element table for a single set of orbitals.
double	dSigma_dE (const DiracSpinor &v, const MBPT::CorrelationPotential &Sigma0, const MBPT::CorrelationPotential &Sigma_plus, const MBPT::CorrelationPotential &Sigma_minus, double delta)
	Normalisation correction factor for a valence state, from the energy derivative of the correlation potential.
bool	allowed_K (int K, int kt, int ks, int twoJb, int twoJa)
	Is rank K allowed: triangle rules for the operators (kt, ks) and states.
int	smallest_allowed_K (int kt, int ks, int twoJb, int twoJa)
	The smallest rank K allowed for the amplitude; negative if there is none.
double	sos_valence (int K, const DiracSpinor &Fb, const DiracSpinor &Fa, const DiracOperator::TensorOperator *t, const DiracOperator::TensorOperator *s, double omega, double omega_s, const std::vector< DiracSpinor > &spectrum, const ExternalField::CorePolarisation *dVt=nullptr, const ExternalField::CorePolarisation *dVs=nullptr, const Coulomb::meTable< double > &t_me={}, const Coulomb::meTable< double > &s_me={})
	Valence part of the second-order amplitude A^K, by sum-over-states.
double	sos_core (int K, int twoJ, const DiracOperator::TensorOperator *t, const DiracOperator::TensorOperator *s, double omega, double omega_s, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &excited, const ExternalField::CorePolarisation *dVt=nullptr, const ExternalField::CorePolarisation *dVs=nullptr)
	Contribution to A^K from the polarisation of the closed core, by sum-over-states.
std::pair< double, double >	ms_valence (int K, const DiracSpinor &Fb, const DiracSpinor &Fa, const DiracOperator::TensorOperator *t, const DiracOperator::TensorOperator *s, double omega, double omega_s, const ExternalField::TDHF *dVt, const ExternalField::TDHF *dVs, const MBPT::CorrelationPotential *Sigma=nullptr, const std::vector< DiracSpinor > &project_onto={}, std::ostream &outstream=std::cout)
	Valence part of the second-order amplitude A^K, evaluated with mixed states (TDHF method).
double	ms_core (int K, int twoJ, const DiracOperator::TensorOperator *t, const DiracOperator::TensorOperator *s, double omega, double omega_s, const std::vector< DiracSpinor > &core, const ExternalField::TDHF *dVt, const ExternalField::TDHF *dVs, const MBPT::CorrelationPotential *Sigma=nullptr)
	Contribution to A^K from the polarisation of the closed core, evaluated with mixed states.
double	sos_ci (int K, const CI::PsiJPi &Psi_b, std::size_t ib, const CI::PsiJPi &Psi_a, std::size_t ia, const DiracOperator::TensorOperator *t, const Coulomb::meTable< double > &t_me, const DiracOperator::TensorOperator *s, const Coulomb::meTable< double > &s_me, double omega, double omega_s, const std::vector< CI::PsiJPi > &ciwfs, const std::vector< CI::Level > &levels_to_remove={}, std::ostream &outstream=std::cout)
	Second-order amplitude \( A^K \) between two CI states, by sum-over-states over the solved CI levels.

Enumeration Type Documentation

Frequency

enum class Amplitudes::Frequency

strong

Which frequency the operator, or the RPA, is evaluated at.

There is no obvious default: the right choice depends on the calculation, so it must be given explicitly.

• transition: the driver sets the frequency itself, to each pair's own transition frequency \( \omega_{ab} = \en_a - \en_b \). The operator is updated (h at \( +|\omega_{ab}| \), h_minus at \( -|\omega_{ab}| \)), and the RPA is re-solved, at every element. This is the physically correct frequency for a transition.
• fixed: the driver does not touch the frequency. The operator is assumed to already be at the intended frequency, and the RPA to already have been solved there, by the caller. Use for a fixed external field, or to reproduce a fixed-frequency calculation.

Function Documentation

matrix_element()

MEdata Amplitudes::matrix_element	(	const DiracSpinor &	a,
		const DiracSpinor &	b,
		const DiracOperator::TensorOperator *	h,
		const DiracOperator::TensorOperator *	h_minus = nullptr,
		const ExternalField::CorePolarisation *	dV = nullptr,
		DiracOperator::MatrixElementType	type = DiracOperator::MatrixElementType::Reduced,
		double	omega = 0.0
	)

Single matrix element of h between states a and b, with optional RPA.

Pure evaluation: assumes h, h_minus, and dV are already at the correct frequency. If omega is negative and h_minus is given, h_minus is used for the matrix element (sign-sensitive frequency-dependent operators, e.g. E1v: h holds \( t_+ \) at \( +|\omega| \), h_minus holds \( t_- \) at \( -|\omega| \)).

The MatrixElementType factor is calculated with h (it is purely angular), and stored in the returned MEdata rather than applied.

Parameters

a,b	States: <a||h||b>.
h	The tensor operator.
h_minus	Operator at negative frequency; nullptr if not required.
dV	RPA correction, already solved; nullptr for none.
type	Form of matrix element (Reduced, Stretched, HFConstant).
omega	Transition frequency (only selects h vs h_minus, and is recorded in the output).

Returns

MEdata holding t0, dv, factor, omega, and labels.

set_operator_frequency()

void Amplitudes::set_operator_frequency	(	DiracOperator::TensorOperator *	h,
		DiracOperator::TensorOperator *	h_minus,
		double	omega
	)

Sets the \( t_\pm \) operator pair to the frequency w.

Sets h to \( +|\omega| \) and, if given, h_minus to \( -|\omega| \); does nothing for frequency-independent operators, for which updateFrequency() must not be called. Call this before matrix_elements when using Frequency::fixed.

Note

Only for the \( t_\pm \) pair. Where there is no \( t_- \) (e.g. sr_matrix_elements), the operator takes the signed frequency directly, not \( |\omega| \).

matrix_elements() [1/2]

std::vector< MEdata > Amplitudes::matrix_elements	(	const std::vector< DiracSpinor > &	a_orbs,
		const std::vector< DiracSpinor > &	b_orbs,
		DiracOperator::TensorOperator *	h,
		DiracOperator::TensorOperator *	h_minus,
		ExternalField::CorePolarisation *	dV,
		const MEoptions &	options,
		std::ostream &	outstream = std::cout
	)

Matrix elements of h for all allowed pairs from two lists of orbitals, with optional RPA; owns all frequency updates and RPA solves.

Calculates \( \redmatel{a}{h}{b} \) for each pair allowed by the selection rules, with the bra states taken from a_orbs and the ket states from b_orbs, diagonal first, then off-diagonal.

Selection rules: pairs with isZero() are skipped; diagonal elements only for even-parity operators.

When a_orbs and b_orbs are the same list (as in the single-list overload below), each pair is calculated once: for odd-parity operators, only elements with the even-parity state on the right are included (unless calculate_both); for even-parity operators, only the upper triangle. Two distinct lists have no such pairing, so every pair is calculated and calculate_both has no effect.

Frequency handling (see Frequency): with transition, the operator is updated at each pair's transition frequency, h at \( +|\omega_{ab}| \) and h_minus at \( -|\omega_{ab}| \), and the RPA is re-solved there (cleared first when poorly converged, or when rpa_iterations is 1 so that first-order RPA is not iterated from a previous solution). With fixed, neither is touched: the caller must set the operator frequency (see set_operator_frequency) and solve the RPA before calling.

Parameters

a_orbs	Bra states (index a).
b_orbs	Ket states (index b).
h	The tensor operator.
h_minus	Operator at negative frequency (e.g. a clone of h for E1v); nullptr if not required.
dV	RPA. nullptr for no RPA.
options	See MEoptions.
outstream	Stream for progress output.

Returns

Vector of MEdata, one per calculated matrix element.

me_table() [1/4]

Coulomb::meTable< double > Amplitudes::me_table	(	const std::vector< DiracSpinor > &	a_orbs,
		const std::vector< DiracSpinor > &	b_orbs,
		const DiracOperator::TensorOperator *	h,
		const ExternalField::CorePolarisation *	dV = nullptr
	)

Builds a lookup table of reduced matrix elements <a||h||b>.

Fills and returns a Coulomb::meTable<double> with reduced matrix elements

\[ t_{ab} = \redmatel{a}{h}{b} + \delta V_{ab} \]

for all non-zero pairs from a_orbs and b_orbs.

The symmetry-conjugate \( \redmatel{b}{h}{a} \) is also stored, via symm_sign(). Filled with OpenMP parallelisation.

This is a pure table builder: h must already be at the intended frequency, and dV already solved there.

Parameters

a_orbs	Bra states.
b_orbs	Ket states.
h	Pointer to the (const) tensor operator.
dV	Optional RPA correction. If nullptr, not applied.

Returns

meTable containing t_ab for all non-zero pairs (and conjugates).

me_table() [2/4]

Coulomb::meTable< double > Amplitudes::me_table	(	const std::vector< DiracSpinor > &	a_orbs,
		const std::vector< DiracSpinor > &	b_orbs,
		const DiracOperator::TensorOperator *	h,
		const ExternalField::CorePolarisation *	dV,
		const MBPT::StructureRad *	srn,
		double	omega,
		int	sr_n_max = 999,
		bool	sr_norm = true
	)

Builds a table of reduced matrix elements, including structure radiation and (optionally) the normalisation of states.

As above, but each element also carries the second-order corrections,

\[ t_{ab} = \redmatel{a}{h}{b} + \delta V_{ab} + \delta_{\rm SR}^{ab}. \]

Parameters

a_orbs	Bra states.
b_orbs	Ket states.
h	Pointer to the (const) tensor operator.
dV	Optional RPA correction. If nullptr, not applied.
srn	Structure radiation/normalisation. If nullptr, not applied (and the table is as the plain overload above).
omega	Frequency for the structure radiation denominators. Their frequency dependence is very weak, so in practice this is usually taken as the frequency the RPA was solved at.
sr_n_max	SR+N is applied only to pairs with both n <= sr_n_max. SR+N is meaningful only between physical states, so this limits it to the low-n part of a large basis, where the states are not cavity states. Does not affect the internal lines of the diagrams (see MBPT::StructureRad) [999].
sr_norm	If false, only the structure radiation is added, not the normalisation of states [true].

Returns

meTable containing t_ab for all non-zero pairs (and conjugates).

sr_matrix_elements()

std::vector< SRNdata > Amplitudes::sr_matrix_elements	(	const std::vector< DiracSpinor > &	orbs,
		DiracOperator::TensorOperator *	h,
		MBPT::StructureRad *	sr,
		ExternalField::CorePolarisation *	dV,
		const SRNoptions &	options,
		std::ostream &	outstream = std::cout
	)

Matrix elements with structure radiation, normalisation, and Brueckner orbital corrections, for all pairs from a list of orbitals.

For each pair allowed by the selection rules (as matrix_elements), evaluates the lowest-order matrix element, RPA, and the second-order MBPT corrections via MBPT::StructureRad: SR (top+bottom+centre diagrams), normalisation of states, and (optionally) the Brueckner orbital correction.

Frequency handling (see Frequency): with operator_omega = transition, the operator is evaluated at each pair's transition frequency, where its frequency dependence is important, and its BO term then includes the frequency-derivative correction \( ({\rm d}t/{\rm d}\omega)\,\delta\omega^{(2)} \) for off-diagonal elements. With rpa_omega = transition, the RPA and the SR tables are re-solved at each transition frequency; with fixed, the caller must have solved the RPA already, and the SR denominators are taken at the frequency it was solved at (zero if there is no RPA).

The caller constructs (and owns) the MBPT::StructureRad object, which holds the basis, Qk integrals, and screening options; solve_core is called here, and re-called whenever the frequency changes.

Parameters

orbs	Orbitals for the external legs; all pairs considered.
h	The tensor operator.
sr	StructureRad object (mutated: solve_core is called).
dV	RPA. nullptr for no RPA.
options	See SRNoptions.
outstream	Stream for per-element progress output.

Returns

Vector of SRNdata, one per calculated matrix element.

matrix_elements() [2/2]

std::vector< MEdata > Amplitudes::matrix_elements ( const std::vector< DiracSpinor > &  orbs, DiracOperator::TensorOperator *  h, DiracOperator::TensorOperator *  h_minus, ExternalField::CorePolarisation *  dV, const MEoptions &  options, std::ostream &  outstream = std::cout  )

inline

Matrix elements of h for all pairs from a single list of orbitals.

Convenience overload; calls matrix_elements(orbs, orbs, ...) with both bra and ket taken from orbs, so each pair is calculated once.

me_table() [3/4]

Coulomb::meTable< double > Amplitudes::me_table ( const std::vector< DiracSpinor > &  a_orbs, const DiracOperator::TensorOperator *  h, const ExternalField::CorePolarisation *  dV = nullptr  )

inline

Builds a matrix element table for a single set of orbitals.

Convenience overload; calls me_table(a_orbs, a_orbs, ...) with both bra and ket taken from a_orbs.

me_table() [4/4]

Coulomb::meTable< double > Amplitudes::me_table ( const std::vector< DiracSpinor > &  a_orbs, const DiracOperator::TensorOperator *  h, const ExternalField::CorePolarisation *  dV, const MBPT::StructureRad *  srn, double  omega, int  sr_n_max = 999, bool  sr_norm = true  )

inline

Builds a SR+N matrix element table for a single set of orbitals.

Convenience overload; calls me_table(a_orbs, a_orbs, ...) with both bra and ket taken from a_orbs.

dSigma_dE()

double Amplitudes::dSigma_dE	(	const DiracSpinor &	v,
		const MBPT::CorrelationPotential &	Sigma0,
		const MBPT::CorrelationPotential &	Sigma_plus,
		const MBPT::CorrelationPotential &	Sigma_minus,
		double	delta
	)

Normalisation correction factor for a valence state, from the energy derivative of the correlation potential.

The normalisation of a Brueckner orbital differs from unity at second order; the correction factor for state \( v \) is

\[ \frac{{\rm d}\Sigma_v}{{\rm d}\en} = \lambda_v \frac{\matel{v}{\Sigma(\en_v+\delta) - \Sigma(\en_v-\delta)}{v}} {2\delta}, \]

evaluated by central difference, with \( \lambda_v \) the scaling factor of the correlation potential. The normalisation correction to a matrix element is then \( \delta t^{\rm Norm}_{ab} = \frac{1}{2}(t_{ab} + \delta V_{ab}) ({\rm d}\Sigma_a/{\rm d}\en + {\rm d}\Sigma_b/{\rm d}\en) \). This is the non-perturbative alternative to the sum-over-states normalisation of MBPT::StructureRad::norm.

Parameters

v	Valence state.
Sigma0	The correlation potential of the wavefunction (supplies the lambda scaling factor).
Sigma_plus	Correlation potential formed at e_v + delta.
Sigma_minus	Correlation potential formed at e_v - delta.
delta	The energy step the potentials were formed at.

Returns

d Sigma_v / d e (dimensionless; zero if the potentials hold no Sigma of this kappa).

Note

getSigma falls back to the nearest-n Sigma of matching kappa, so a state the potentials were not formed for still returns a value (evaluated with that Sigma).

allowed_K()

bool Amplitudes::allowed_K	(	int	K,
		int	kt,
		int	ks,
		int	twoJb,
		int	twoJa
	)

Is rank K allowed: triangle rules for the operators (kt, ks) and states.

Second-order (in the external field) amplitudes for a single-valence atom.

For a transition \( a \to b \) due to two one-body operators, \( t \) (rank \( k_t \), frequency \( \omega \)) and \( s \) (rank \( k_s \), frequency \( \omega_s \)), the amplitude with definite projections \( q_1, q_2 \) of the two operators is

\[ A^{k_tk_s}_{q_1q_2} = \sum_n \left[ \frac{\matel{b}{t_{q_1}}{n}\matel{n}{s_{q_2}}{a}}{\en_a + \omega_s - \en_n} + \frac{\matel{b}{s_{q_2}}{n}\matel{n}{t_{q_1}}{a}}{\en_a + \omega - \en_n} \right], \]

the sum over \( n \) running over the magnetic quantum numbers too, and \( \en_b = \en_a + \omega + \omega_s \). The operators are coupled to rank \( K \), with \( Q = q_1+q_2 = m_b-m_a \) and \( [K] \equiv 2K+1 \),

\[ A^K_Q = \sum_{q_1q_2}\braket{k_tq_1\,k_sq_2}{KQ}\,A^{k_tk_s}_{q_1q_2} = (-1)^{k_t-k_s+Q}\sqrt{[K]}\sum_{q_1q_2} \threej{k_t}{k_s}{K}{q_1}{q_2}{-Q}\,A^{k_tk_s}_{q_1q_2}, \]

and \( A^K \) follows from the Wigner-Eckart theorem,

\[ A^K_Q = (-1)^{j_b-m_b}\threej{j_b}{K}{j_a}{-m_b}{Q}{m_a}\,A^K \qquad {\rm with} \qquad A^K = \sum_n \left[ c_1(j_n)\, \frac{\redmatel{b}{t}{n}\redmatel{n}{s}{a}}{\en_a + \omega_s - \en_n} + c_2(j_n)\, \frac{\redmatel{b}{s}{n}\redmatel{n}{t}{a}}{\en_a + \omega - \en_n} \right], \]

the coefficients being those of CI::A_K_coefs (evaluated with the single-particle j in place of J), which also gives the sign convention of the coupling and the specific cases; CI::z_component converts \( A^K \) to the z-component. With \( t = s = d \) (E1) and \( [j] \equiv 2j+1 \):

\[ \alpha_0 = \frac{A^0}{\sqrt{3[j_a]}}, \qquad \alpha_2 = -\sqrt{\frac{2j(2j-1)}{3(j+1)(2j+1)(2j+3)}}\;A^2, \qquad \beta = \frac{A^1}{\sqrt{2}\,\redmatel{b}{\bm\sigma}{a}}, \]

for \( K = 0, 2, 1 \) ( \( \alpha_2 \) requires \( j_b = j_a = j \ge 1 \)), and with \( t = d \), \( s = h_W \) (PNC, \( k_s = 0 \), so \( K = 1 \)),

\[ E_{\rm PNC} = A^1_0 = (-1)^{j_b-m}\threej{j_b}{1}{j_a}{-m}{0}{m}\,A^1 . \]

This is the single-valence analogue of CI::A_K; it covers static, dynamic, and transition polarisabilities ( \( t = s = E1 \)), and PNC amplitudes ( \( s \) = PNC operator).

Two methods are provided for the valence sum: sum-over-states (SOS) over a given spectrum, and mixed states (MS, solving the inhomogeneous equation via ExternalField::TDHF), which is complete (no truncation of the sum). The contribution of core excitations (closed core: K = 0, diagonal only) is likewise available both ways. SOS and MS must agree (to basis completeness); the comparison is a strong check on the numerics.

smallest_allowed_K()

int Amplitudes::smallest_allowed_K	(	int	kt,
		int	ks,
		int	twoJb,
		int	twoJa
	)

The smallest rank K allowed for the amplitude; negative if there is none.

sos_valence()

double Amplitudes::sos_valence	(	int	K,
		const DiracSpinor &	Fb,
		const DiracSpinor &	Fa,
		const DiracOperator::TensorOperator *	t,
		const DiracOperator::TensorOperator *	s,
		double	omega,
		double	omega_s,
		const std::vector< DiracSpinor > &	spectrum,
		const ExternalField::CorePolarisation *	dVt = nullptr,
		const ExternalField::CorePolarisation *	dVs = nullptr,
		const Coulomb::meTable< double > &	t_me = {},
		const Coulomb::meTable< double > &	s_me = {}
	)

Valence part of the second-order amplitude A^K, by sum-over-states.

Evaluates the sum above directly over the states of spectrum.

If the spectrum contains states below the Fermi level (e.g. Wavefunction::spectrum() does), the core-valence (Pauli blocking) part of the amplitude is included automatically through those terms, as in the polarisability module. The polarisation of the closed core is separate: see sos_core.

Matrix elements: taken from t_me / s_me if present in the table, otherwise calculated directly as \( \redmatel{}{h}{} + \delta V \). The tables (if given) should be formed at the operator's frequency, and may contain RPA and structure radiation. RPA enters on both vertices.

Parameters

K	Rank of the amplitude.
Fb,Fa	Final and initial valence states.
t,s	The two operators.
omega	Frequency of t. For a real transition carried entirely by t this is e_b - e_a, and s is static.
omega_s	Frequency of s. Energy conservation requires omega + omega_s = e_b - e_a; it is -omega for a dynamic polarisability (b = a).
spectrum	Intermediate states summed over.
dVt,dVs	RPA for each operator, solved at the frequency of that operator. May be nullptr. Ignored for table entries.
t_me,s_me	Optional tables of single-particle reduced matrix elements; empty (default) to calculate directly.

Returns

\( A^K \) (valence part).

Note

Degenerate denominators arise for even-parity operator pairs (an intermediate state degenerate with \( \en_a + \omega \); e.g. n = a for a static diagonal amplitude). Such terms diverge, and nothing is skipped here: remove the offending states from the sum and treat them separately (cf. the project-out treatment in the pnc module).

sos_core()

double Amplitudes::sos_core	(	int	K,
		int	twoJ,
		const DiracOperator::TensorOperator *	t,
		const DiracOperator::TensorOperator *	s,
		double	omega,
		double	omega_s,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	excited,
		const ExternalField::CorePolarisation *	dVt = nullptr,
		const ExternalField::CorePolarisation *	dVs = nullptr
	)

Contribution to A^K from the polarisation of the closed core, by sum-over-states.

Delegates to CI::A_K_core, which is the same quantity: a core electron excited by one operator and de-excited by the other,

\[ A^0_{\rm core} = \sqrt{[J]}\sum_c \sqrt{[j_c]} \sum_m \left[ c_1\,\frac{\redmatel{c}{t}{m}\redmatel{m}{s}{c}}{\en_c+\omega_s-\en_m} + c_2\,\frac{\redmatel{c}{s}{m}\redmatel{m}{t}{c}}{\en_c+\omega-\en_m} \right]. \]

Non-zero only for K = 0 (closed core), which requires kt = ks, and only for a diagonal amplitude (b = a); the caller is responsible for the diagonal condition.

Parameters

K	Rank; zero returned unless K = 0.
twoJ	2J of the valence state (the sqrt([J]) prefactor).
t,s	The two operators.
omega,omega_s	Frequency of each operator; see sos_valence.
core	Core states c.
excited	Particle states m: basis states above the Fermi level.
dVt,dVs	RPA for each operator. May be nullptr.

Returns

\( A^0_{\rm core} \).

Note

RPA enters once: of the two matrix elements, only the one acting on the core orbital is dressed (see the note on CI::A_K_core).

ms_valence()

std::pair< double, double > Amplitudes::ms_valence	(	int	K,
		const DiracSpinor &	Fb,
		const DiracSpinor &	Fa,
		const DiracOperator::TensorOperator *	t,
		const DiracOperator::TensorOperator *	s,
		double	omega,
		double	omega_s,
		const ExternalField::TDHF *	dVt,
		const ExternalField::TDHF *	dVs,
		const MBPT::CorrelationPotential *	Sigma = nullptr,
		const std::vector< DiracSpinor > &	project_onto = {},
		std::ostream &	outstream = std::cout
	)

Valence part of the second-order amplitude A^K, evaluated with mixed states (TDHF method).

The sums over intermediate states are performed with mixed states (solutions of the inhomogeneous Dirac equation, via ExternalField::TDHF::solve_dPsi), so they are complete: no truncation of the spectrum. All intermediate states of a given kappa share the same angular coefficient, so one mixed state per kappa channel per term is required.

Each sum is formed in two independent ways: with the mixed states of \( s \), and with those of \( t \), returned as the two elements of the pair. They must agree; the difference is a check on the numerics (cf. CI::A_K).

RPA enters on both vertices: the mixed states are solved with \( t + \delta V \) (the inner vertex), and the outer matrix element is dressed with the \( \delta V \) of the outer operator. With RPA not solved, \( \delta V = 0 \) and the amplitude is at the HF level.

Parameters

K	Rank of the amplitude.
Fb,Fa	Final and initial valence states.
t,s	The two operators. The core-valence (Pauli blocking) part of the sum is included (the mixed states are complete). It may be separated by projecting the mixed states onto the span of the core states: pass the HF core as project_onto to obtain just that part (cf. the orthogonality treatment in the pnc module).
omega,omega_s	Frequency of each operator; see sos_valence.
dVt,dVs	TDHF object for each operator (required, not null): it provides the mixed-state solver even when RPA is off (unsolved TDHF gives dV = 0). Solve at the operator's frequency for RPA. Must be TDHF (or TDHFbasis): the diagram method cannot drive mixed states.
Sigma	Optional correlation potential, included in the mixed-state solutions (use with Brueckner valence states).
project_onto	If non-empty, each mixed state is projected onto the span of these states before the outer matrix element; pass the HF core for the core-valence part of the amplitude. Empty (default): the full sum.
outstream	Stream for warnings.

Returns

The two evaluations of \( A^K \): with the mixed states of s, and with those of t.

Note

If an intermediate state is degenerate with a denominator, the mixed-state equation is singular and the amplitude divergent. This cannot happen for odd-parity operators (the intermediate states have opposite parity to a and b); see the note on CI::A_K.

ms_core()

double Amplitudes::ms_core	(	int	K,
		int	twoJ,
		const DiracOperator::TensorOperator *	t,
		const DiracOperator::TensorOperator *	s,
		double	omega,
		double	omega_s,
		const std::vector< DiracSpinor > &	core,
		const ExternalField::TDHF *	dVt,
		const ExternalField::TDHF *	dVs,
		const MBPT::CorrelationPotential *	Sigma = nullptr
	)

Contribution to A^K from the polarisation of the closed core, evaluated with mixed states.

As sos_core, with the sum over excited states m performed with mixed states of the core orbitals (generalises the TDHF core polarisability of the polarisability module to two operators). Non-zero only for K = 0 and a diagonal amplitude (caller's responsibility).

RPA enters once, through the mixed-state solve (the vertex acting on the core orbital); the outer matrix element is bare - the same counting as sos_core.

The mixed states of a core orbital include (occupied) core intermediate states; for the diagonal amplitude (omega_s = -omega) their contributions cancel in the sum over the core, so this agrees with sos_core, which sums over excited states only.

Parameters

K	Rank; zero returned unless K = 0.
twoJ	2J of the valence state.
t,s	The two operators.
omega,omega_s	Frequency of each operator.
core	Core states c.
dVt,dVs	TDHF object for each operator (required); see ms_valence.
Sigma	Optional correlation potential for the mixed states. The intermediate states here are core excitations, for which the valence Sigma is not appropriate: normally nullptr (matching sos_core, which uses the HF basis).

Returns

\( A^0_{\rm core} \).

sos_ci()

double Amplitudes::sos_ci	(	int	K,
		const CI::PsiJPi &	Psi_b,
		std::size_t	ib,
		const CI::PsiJPi &	Psi_a,
		std::size_t	ia,
		const DiracOperator::TensorOperator *	t,
		const Coulomb::meTable< double > &	t_me,
		const DiracOperator::TensorOperator *	s,
		const Coulomb::meTable< double > &	s_me,
		double	omega,
		double	omega_s,
		const std::vector< CI::PsiJPi > &	ciwfs,
		const std::vector< CI::Level > &	levels_to_remove = {},
		std::ostream &	outstream = std::cout
	)

Second-order amplitude \( A^K \) between two CI states, by sum-over-states over the solved CI levels.

Evaluates \( A^K \) (see CI::A_K_coefs) directly,

\[ A^K = \sum_n \left[ c_1(J_n)\, \frac{\redmatel{b}{t}{n}\redmatel{n}{s}{a}}{E_a + \omega_s - E_n} + c_2(J_n)\, \frac{\redmatel{b}{s}{n}\redmatel{n}{t}{a}}{E_a + \omega - E_n} \right], \]

the sum running over the solutions of each allowed intermediate (J, parity). This is the sum-over-states analogue of CI::A_K (the generalisation of the CI_Pol module): only the levels actually solved in the CI{} block are available, so the sum is truncated - both to the (J, parity) blocks that were solved (a note is printed for any that are missing) and to num_solutions levels within each. CI::A_K, by contrast, is complete; the comparison shows how much of the amplitude the low levels carry.

The reduced matrix elements are the CI contractions of the single-particle tables (CI::ReducedME). Corrections to the matrix elements (RPA, structure radiation, normalisation of states) enter through those tables.

Parameters

K	Rank of the amplitude.
Psi_b,ib	Final CI state (solution ib of Psi_b).
Psi_a,ia	Initial CI state.
t,t_me	The \( t \) operator, and its table of single-particle reduced matrix elements (which may include RPA, structure radiation, normalisation), formed at omega.
s,s_me	The \( s \) operator, and its table (formed at omega_s).
omega,omega_s	Frequency of each operator; see CI::A_K.
ciwfs	The solved CI blocks, e.g., Wavefunction::CIwfs(): the intermediate states of the sum.
levels_to_remove	CI levels skipped in the sum, so that they may be treated separately - e.g., with experimental energies.
outstream	Stream for the per-(J, parity) contributions.

Returns

\( A^K \) (valence part).

Note

The polarisation of the closed core and the core-valence terms lie outside the CI space: see CI::A_K_core and CI::A_K_cv.

A level degenerate with a denominator ( \( E_a + \omega_s \) or \( E_a + \omega \)) makes its term diverge, and nothing is skipped here: remove the offending level (levels_to_remove) and treat it separately.