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: the lowest-order value and each correction, separately. More...
struct	MEoptions
	Options for the matrix element drivers. More...

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

Functions
void	solve_at_frequency (DiracOperator::TensorOperator *h, ExternalField::CorePolarisation *dV, MBPT::StructureRad *sr, double omega, const MEoptions &options)
	Sets the operator frequency, solves the RPA, and re-builds the structure radiation tables, for a transition at frequency omega, as the options ask. The one place the frequencies are handled.
MEdata	matrix_element (const DiracSpinor &a, const DiracSpinor &b, const DiracOperator::TensorOperator *h, const ExternalField::CorePolarisation *dV=nullptr, const MBPT::StructureRad *sr=nullptr, DiracOperator::MatrixElementType type=DiracOperator::MatrixElementType::Reduced, double omega=0.0, double omega_sr=0.0)
	Single matrix element of h between states a and b, with its RPA and structure radiation corrections.
double	sr_frequency (const ExternalField::CorePolarisation *dV, double omega, const MEoptions &options)
	The frequency for the structure radiation denominators of a transition at omega: the transition frequency if the RPA follows each transition, otherwise the frequency the RPA was solved at (zero if there is no RPA).
double	frequency_derivative_term (const DiracSpinor &a, const DiracSpinor &b, DiracOperator::TensorOperator *h, const MBPT::StructureRad *sr, double omega)
	Frequency-derivative term of <a||h||b> for a frequency-dependent operator: (dt/dw) dw2, the shift from the second-order energy shifts dw2 moving the transition frequency.
std::vector< std::pair< const DiracSpinor *, const DiracSpinor * > >	orbital_pairs (const std::vector< DiracSpinor > &a_orbs, const std::vector< DiracSpinor > &b_orbs, const DiracOperator::TensorOperator *h, const MEoptions &options)
	The pairs of orbitals (a, b) of <a||h||b> that the drivers calculate, in order.
std::vector< MEdata >	matrix_elements (const std::vector< DiracSpinor > &a_orbs, const std::vector< DiracSpinor > &b_orbs, DiracOperator::TensorOperator *h, ExternalField::CorePolarisation *dV, MBPT::StructureRad *sr, const MEoptions &options, std::ostream &outstream=std::cout)
	Matrix elements of h for all allowed pairs from two lists of orbitals, with optional RPA and structure radiation; 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)
	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< MEdata >	matrix_elements (const std::vector< DiracSpinor > &orbs, DiracOperator::TensorOperator *h, ExternalField::CorePolarisation *dV, MBPT::StructureRad *sr, 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)
	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 \); see solve_at_frequency. 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 (the structure radiation tables are then built once, by the driver). Use for a fixed external field, or to reproduce a fixed-frequency calculation.

Function Documentation

solve_at_frequency()

void Amplitudes::solve_at_frequency	(	DiracOperator::TensorOperator *	h,
		ExternalField::CorePolarisation *	dV,
		MBPT::StructureRad *	sr,
		double	omega,
		const MEoptions &	options
	)

Sets the operator frequency, solves the RPA, and re-builds the structure radiation tables, for a transition at frequency omega, as the options ask. The one place the frequencies are handled.

In this order:

• with operator_omega = transition: the operator is set to \( +|\omega| \), since the RPA must be solved with the operator at positive frequency;
• with rpa_omega = transition: the RPA is re-solved at omega (cleared first if the previous solution is poorly converged, or if first-order RPA was asked for, which must not iterate from a previous solution);
• with operator_omega = transition: the operator is then set to the signed omega, its physical frequency for this transition;
• with structure radiation, and either of the above: the SR tables are re-built, so that their internal lines carry the operator and the RPA as just set.

With both frequencies fixed, nothing is touched: the caller has set the operator, solved the RPA, and built the SR tables.

Parameters

h	The operator (frequency updated if frequency-dependent).
dV	RPA; nullptr for none.
sr	Structure radiation; nullptr for none.
omega	The transition frequency, signed: e_a - e_b for <a||h||b>.
options	See MEoptions (the two Frequency choices, rpa_iterations, print).

matrix_element()

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

Single matrix element of h between states a and b, with its RPA and structure radiation corrections.

Pure evaluation: assumes h is at the signed transition frequency, and dV and sr solved there (see solve_at_frequency). The RPA correction is the component at the transition frequency: the emission one ( \( \delta V^\dagger \)) if omega is negative; see CorePolarisation::dV. The structure radiation denominators are taken at omega_sr: its frequency dependence is very weak, so when the RPA is at a fixed frequency this is usually the frequency the RPA was solved at (see sr_frequency). The Brueckner orbital correction (MBPT::StructureRad::BO) is not included here.

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.
dV	RPA correction, already solved; nullptr for none.
sr	Structure radiation, tables already built; nullptr for none.
type	Form of matrix element (Reduced, Stretched, HFConstant).
omega	Transition frequency, e_a - e_b: selects the RPA component, and is recorded in the output.
omega_sr	Frequency for the structure radiation denominators.

Returns

MEdata holding t0, dv, sr, norm, factor, omega, labels.

sr_frequency()

double Amplitudes::sr_frequency	(	const ExternalField::CorePolarisation *	dV,
		double	omega,
		const MEoptions &	options
	)

The frequency for the structure radiation denominators of a transition at omega: the transition frequency if the RPA follows each transition, otherwise the frequency the RPA was solved at (zero if there is no RPA).

The structure radiation follows the RPA: its frequency dependence is very weak, so it is taken at the same frequency.

frequency_derivative_term()

double Amplitudes::frequency_derivative_term	(	const DiracSpinor &	a,
		const DiracSpinor &	b,
		DiracOperator::TensorOperator *	h,
		const MBPT::StructureRad *	sr,
		double	omega
	)

Frequency-derivative term of <a||h||b> for a frequency-dependent operator: (dt/dw) dw2, the shift from the second-order energy shifts dw2 moving the transition frequency.

\( \delta t = ({\rm d}t/{\rm d}\omega)\,\delta\omega^{(2)} \). Off-diagonal elements only; zero for a frequency-independent operator. The derivative is taken numerically, which is why the operator is not const. It is left as it was, at omega.

Parameters

a,b	States: <a||h||b>.
h	The operator, at the signed transition frequency omega.
sr	Structure radiation (holds Sigma for the energy shifts).
omega	Transition frequency, e_a - e_b.

orbital_pairs()

std::vector< std::pair< const DiracSpinor *, const DiracSpinor * > > Amplitudes::orbital_pairs	(	const std::vector< DiracSpinor > &	a_orbs,
		const std::vector< DiracSpinor > &	b_orbs,
		const DiracOperator::TensorOperator *	h,
		const MEoptions &	options
	)

The pairs of orbitals (a, b) of <a||h||b> that the drivers calculate, in order.

Diagonal elements first (only for even-parity operators, and only for states in both lists), then off-diagonal; pairs forbidden by the selection rules are skipped. When a_orbs and b_orbs are the same list, each pair is calculated once: for odd-parity operators only elements with the even-parity state on the right, for even-parity operators only the upper triangle (unless calculate_both). Two distinct lists have no such pairing, so every pair is included.

Parameters

a_orbs	Bra states (a).
b_orbs	Ket states (b).
h	The tensor operator (selection rules).
options	diagonal, off_diagonal, calculate_both.

Returns

Pointers into a_orbs and b_orbs: {a, b} for each element.

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,
		ExternalField::CorePolarisation *	dV,
		MBPT::StructureRad *	sr,
		const MEoptions &	options,
		std::ostream &	outstream = std::cout
	)

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

For each pair of orbital_pairs: solve_at_frequency, then matrix_element. Nothing is printed beyond the RPA progress (if options.print, which also prints each transition as the RPA is re-solved).

Parameters

a_orbs	Bra states (index a).
b_orbs	Ket states (index b).
h	The tensor operator.
dV	RPA. nullptr for no RPA.
sr	Structure radiation (mutated: its tables are built). nullptr for none.
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
	)

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

Fills and returns a Coulomb::meTable<double> with the bare reduced matrix elements \( \redmatel{a}{h}{b} \) 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. h must already be at the intended frequency. For a table with RPA (and structure radiation), use the overload taking omega.

Parameters

a_orbs	Bra states.
b_orbs	Ket states.
h	Pointer to the (const) tensor operator.

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}, \]

and every element is at the one frequency omega: the RPA correction is \( \delta V(\omega) \) throughout (the emission component, \( \delta V^\dagger \), for negative omega; see CorePolarisation::dV), and the structure radiation is at omega, absorbed by the ket. The reversed pair \( (b, a) \) is then not the symmetry conjugate of \( (a, b) \), and is calculated as well (except at \( \omega = 0 \), where it is).

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 of the table: that of the operator component the table stands for. Sets the sign of the RPA correction, and the structure radiation denominators.
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, both orderings.

matrix_elements() [2/2]

std::vector< MEdata > Amplitudes::matrix_elements ( const std::vector< DiracSpinor > &  orbs, DiracOperator::TensorOperator *  h, ExternalField::CorePolarisation *  dV, MBPT::StructureRad *  sr, 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  )

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, as \( \delta V(\omega) \) of the operator's own frequency at every vertex (the emission component for a negative frequency; see CorePolarisation::dV), whichever side of the leg the intermediate state lies.

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; each is \( \delta V(\omega) \) at the frequency of its own operator (cf. sos_valence). 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/Pi 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 J/Pi, 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.