CI Namespace Reference

Detailed Description

Functions and classes for Configuration Interaction calculations.

Main functions are:

• configuration_interaction
• run_CI

Main Classes are:

• CSF2
• PsiJPi

Classes
struct	ConfigInfo
	Configuration metadata for a single CI level. More...
class	CSF2
	Two-electron configuration state function (CSF). More...
struct	Integrals
	The integral tables required to construct the CI Hamiltonian matrix. More...
struct	Level
	Identifies one CI level: its (J, parity), and which solution. More...
class	PsiJPi
	Container for CI solutions in a single (J, parity) block. More...
struct	Sigma1Correction
	Derivative (dSigma/dE) correction data for the one-body Sigma_1 matrix elements. More...
struct	Solutions
	The result of a CI calculation: the solutions, and the integrals used to construct the CI Hamiltonian. More...

Functions
double	CSF2_Coulomb (const Coulomb::QkTable &qk, DiracSpinor::Index v, DiracSpinor::Index w, DiracSpinor::Index x, DiracSpinor::Index y, int twoJ)
	Antisymmetrised two-body Coulomb matrix element in the coupled CSF basis.
double	CSF2_Sigma2 (const Coulomb::LkTable &Sk, DiracSpinor::Index v, DiracSpinor::Index w, DiracSpinor::Index x, DiracSpinor::Index y, int twoJ, const Coulomb::QkTable *qk=nullptr, const std::vector< double > &hk={}, const Coulomb::LkTable *dSk=nullptr, double dE0=0.0)
	Two-body \( \Sigma_2 \) (MBPT) correction to CSF2_Coulomb().
double	CSF2_Breit (const Coulomb::WkTable &Bk, DiracSpinor::Index v, DiracSpinor::Index w, DiracSpinor::Index x, DiracSpinor::Index y, int twoJ)
	Antisymmetrised two-body Breit matrix element in the coupled CSF basis.
double	Sigma2_AB (const CI::CSF2 &A, const CI::CSF2 &B, int twoJ, const Coulomb::LkTable &Sk, const Coulomb::QkTable *qk=nullptr, const std::vector< double > &hk={}, const Coulomb::LkTable *dSk=nullptr, double dE0=0.0)
	Two-body \( \Sigma_2 \) correction to Hab().
double	Breit_AB (const CI::CSF2 &A, const CI::CSF2 &B, int twoJ, const Coulomb::WkTable &Bk)
	Breit correction to Hab().
double	corrected_Sigma (double Sigma, double dSigma, double dE)
	Resummed derivative (dSigma/dE) correction to a Sigma_1 matrix element (Kozlov formula).
double	corrected_Sk (double Sk, double dSk, double dE0)
	Resummed shift of a \( \Sigma_2 \) integral to a new reference energy E0 (Brillouin-Wigner denominators).
Sigma1Correction	calculate_dSdE_correction (const std::vector< DiracSpinor > &ci_basis, const std::vector< DiracSpinor > &s1_basis_core, const std::vector< DiracSpinor > &s1_basis_excited, const Coulomb::QkTable &qk)
	Builds the Sigma_1 derivative-correction tables; see Sigma1Correction.
Sigma1Correction	calculate_dSdE_correction (const std::vector< DiracSpinor > &ci_basis, const MBPT::CorrelationPotential &Sigma)
	Builds the Sigma_1 derivative-correction tables directly from a correlation potential that holds dSigma/dE matrices.
LinAlg::Matrix< double >	iterate_E0 (PsiJPi *psi, const Sigma1Correction &s1c, const Coulomb::meTable< double > &h1, const Coulomb::QkTable &qk, const Coulomb::WkTable *Bk=nullptr, const Coulomb::LkTable *Sk=nullptr, const std::vector< double > &hk={}, const Coulomb::LkTable *dSk=nullptr, double E0_sigma2=0.0, std::ostream &outstream=std::cout)
	Finds the reference energy E0 for the dSigma/dE correction, for a single (J, parity); returns the CI Hamiltonian built with it.
double	Hab (const CI::CSF2 &A, const CI::CSF2 &B, int twoJ, const Coulomb::meTable< double > &h1, const Coulomb::QkTable &qk, const Sigma1Correction *s1c=nullptr)
	CI Hamiltonian matrix element between two two-electron CSFs.
Coulomb::meTable< double >	calculate_h1_table (const std::vector< DiracSpinor > &ci_basis, const std::vector< DiracSpinor > &s1_basis_core, const std::vector< DiracSpinor > &s1_basis_excited, const Coulomb::QkTable &qk, bool include_Sigma1)
	Builds the one-body Hamiltonian matrix element table for the CI basis.
Coulomb::meTable< double >	calculate_h1_table (const std::vector< DiracSpinor > &ci_basis, const MBPT::CorrelationPotential &Sigma, bool include_Sigma1)
	Builds the one-body Hamiltonian table using a precomputed CorrelationPotential.
Coulomb::WkTable	calculate_Bk (const std::string &bk_filename, const HF::Breit *const pBr, const std::vector< DiracSpinor > &ci_basis, int max_k, bool no_new_integralsQ=false)
	Builds or loads the two-body Breit integral table.
std::vector< DiracSpinor >	basis_subset (const std::vector< DiracSpinor > &basis, const std::string &include_str, const std::string &exclude_str="")
	Returns the subset of basis matching include_str, excluding states in exclude_str.
double	ReducedME (const LinAlg::View< const double > &cA, const std::vector< CI::CSF2 > &CSFAs, int twoJA, const LinAlg::View< const double > &cB, const std::vector< CI::CSF2 > &CSFBs, int twoJB, const Coulomb::meTable< double > &h, int K_rank, int Parity)
	Reduced matrix element between two CI states (low-level overload).
double	RME_CSF2 (const CI::CSF2 &X, int twoJX, const CI::CSF2 &V, int twoJV, const Coulomb::meTable< double > &h, int K_rank)
	Reduced matrix element between two two-electron CSFs.
std::pair< std::string, double >	leading_config (const LinAlg::View< const double > &coefs, const std::vector< CSF2 > &csfs)
	Leading non-relativistic configuration of a CI state.
std::pair< int, int >	Term_S_L (int l1, int l2, int twoJ, double gJ_target)
	Determines the best-fit (S, L) term for a two-electron state by matching the g-factor.
std::pair< int, int >	Term_S_L_from_expectation (double L2, double S2, int twoJ)
	Determines the (S, L) term for a two-electron state from the expectation values of L^2 and S^2.
std::string	Term_Symbol (int two_J, int L, int two_S, int parity)
	Returns spectroscopic term symbol string, e.g. "3P_1".
std::string	Term_Symbol (int L, int two_S, int parity)
	Returns term symbol without the J subscript, e.g. "3P".
LinAlg::Matrix< double >	construct_Hci (const PsiJPi &psi, const Coulomb::meTable< double > &h1, const Coulomb::QkTable &qk, const Coulomb::WkTable *Bk=nullptr, const Coulomb::LkTable *Sk=nullptr, const Sigma1Correction *s1c=nullptr, const std::vector< double > &hk={}, const Coulomb::LkTable *dSk=nullptr, double dE0=0.0)
	Constructs the full CI Hamiltonian matrix in the CSF basis.
LinAlg::Matrix< double >	construct_Hci (const PsiJPi &psi, const Integrals &ints)
	Constructs the CI Hamiltonian matrix from a set of integral tables.
double	ReducedME (const PsiJPi &As, std::size_t iA, const PsiJPi &Bs, std::size_t iB, const Coulomb::meTable< double > &h, int K_rank, int Parity)
	Reduced matrix element between two CI states (PsiJPi overload).
Solutions	configuration_interaction (const IO::InputBlock &input, const Wavefunction &wf)
	Runs Configuration Interation: returns CI solutions for all requested J and parity values.
PsiJPi	run_CI (const std::vector< DiracSpinor > &ci_sp_basis, int twoJ, int parity, int num_solutions, std::optional< double > all_below_cm, const Coulomb::meTable< double > &h1, const Coulomb::QkTable &qk, const Coulomb::WkTable &Bk, const Coulomb::LkTable &Sk, bool include_Sigma2, bool print_details, bool read_only=false, std::ostream &outstream=std::cout, const std::string &ci_fname="", const Sigma1Correction *s1c=nullptr, const std::vector< double > &hk={}, const Coulomb::LkTable *dSk=nullptr, double E0_sigma2=0.0)
	Constructs and solves the CI eigenvalue problem for a single J,pi.
bool	operator== (const CSF2 &A, const CSF2 &B)
bool	operator!= (const CSF2 &A, const CSF2 &B)
std::optional< Level >	parse_level (std::string_view str)
	Parses the text form of a CI level reference; see Level.
std::string	to_string (const Level &level)
	Text form of a CI level reference, e.g., "2+:3"; see Level.
std::vector< CSF2 >	form_CSFs (int twoJ, int parity, const std::vector< DiracSpinor > &cisp_basis)
	Forms all two-electron CSFs with given total J and parity.
double	LS_amplitude (int n1, int l1, int twoj1, int n2, int l2, int twoj2, int L, int S, int twoJ)
	jj -> LS recoupling amplitude for an antisymmetrised two-electron CSF.
std::pair< double, double >	expectation_L2S2 (const LinAlg::View< const double > &coefs, const std::vector< CSF2 > &csfs, int twoJ)
	Expectation values of L^2 and S^2 for a two-electron CI state.
LinAlg::Vector< double >	TPsi_reduced (const std::vector< CSF2 > &CSFs, int twoJ, const PsiJPi &Psi0, std::size_t i0, const Coulomb::meTable< double > &h, int K_rank)
	Action of a one-body operator on a CI state, in the CSF basis.
PsiJPi	solve_mixed_state (const PsiJPi &Psi0, std::size_t i0, const PsiJPi &target, const LinAlg::Matrix< double > &Hci, const Coulomb::meTable< double > &h, int K_rank, double omega=0.0)
	Solves the CI mixed-states (Sternheimer) equation for a one-body operator.
PsiJPi	project_out (PsiJPi dPsi, const PsiJPi &levels, const std::vector< std::size_t > &indices)
	Removes CI levels from a mixed state, so that it is orthogonal to them.
PsiJPi	solve_mixed_state (const PsiJPi &Psi0, std::size_t i0, int twoJ, int parity, const std::vector< DiracSpinor > &ci_sp_basis, const Coulomb::meTable< double > &h, int K_rank, const Coulomb::meTable< double > &h1, const Coulomb::QkTable &qk, const Coulomb::WkTable *Bk=nullptr, const Coulomb::LkTable *Sk=nullptr, double omega=0.0)
	Solves the CI mixed-states equation; constructs the CI matrix internally.
std::pair< double, double >	A_K_coefs (int K, int kt, int ks, int twoJb, int twoJn, int twoJa)
	Angular coefficients of the two terms of the second-order amplitude \( A^K \).
double	z_component (int K, int kt, int ks, int twoJb, int twoJa, int two_m)
	Converts the reduced amplitude \( A^K \) to its contribution to the z-component of the amplitude.
int	symm_sign (const DiracOperator::TensorOperator *h, int twoJA, int twoJB)
	Relative sign between <A||h||B> and <B||h||A>, for CI states with total angular momenta 2J_A and 2J_B (cf DiracOperator::TensorOperator::symm_sign)
double	sigma_rme (const PsiJPi &Psi_b, std::size_t ib, const PsiJPi &Psi_a, std::size_t ia, const std::vector< DiracSpinor > &ci_basis)
	Reduced matrix element of the Pauli spin operator between two CI states, \( \redmatel{b}{\sigma}{a} \).
std::pair< double, double >	A_K (int K, const PsiJPi &Psi_b, std::size_t ib, const 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 Integrals &ints, const std::vector< Level > &levels_to_remove={}, std::ostream &outstream=std::cout)
	Second-order amplitude \( A^K \) between two CI states, evaluated with CI mixed states.
double	A_K_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.
double	A_K_cv (int K, const PsiJPi &Psi_b, std::size_t ib, const PsiJPi &Psi_a, std::size_t ia, const DiracOperator::TensorOperator *t, const DiracOperator::TensorOperator *s, double omega, double omega_s, const std::vector< DiracSpinor > &core, const std::vector< DiracSpinor > &ci_basis, const ExternalField::CorePolarisation *dVt=nullptr, const ExternalField::CorePolarisation *dVs=nullptr)
	Core-valence contribution to \( A^K \): the Pauli blocking of the core excitations by the valence electrons.

---

Class Documentation

CI::ConfigInfo

struct CI::ConfigInfo

Class Members
string	config {}	Dominant configuration label (typically non-relativistic notation)
double	ci2 {0.0}	Squared CI coefficient of the dominant configuration (or sum over non-rel degenerates)
double	gJ {0.0}
double	L {-1.0}	Approximate orbital angular momentum L (-1 if not assigned)
double	twoS {-1.0}	Twice the approximate spin S (-1 if not assigned)
double	L2 {-1.0}	Expectation value of L^2 for the CI state (-1 if not computed)
double	S2 {-1.0}	Expectation value of S^2 for the CI state (-1 if not computed)

CI::Level

struct CI::Level

Class Members
int	twoJ {0}	Twice the total angular momentum, 2J.
int	parity {1}	Parity: +1 or -1.
size_t	index {0}	Which solution, counting from zero, in order of energy.

CI::Solutions

struct CI::Solutions

Class Members
vector< PsiJPi >	levels {}	One entry per {J, parity} requested.
Integrals	integrals {}	Integral tables used to construct the CI Hamiltonians.

Function Documentation

CSF2_Coulomb()

double CI::CSF2_Coulomb	(	const Coulomb::QkTable &	qk,
		DiracSpinor::Index	v,
		DiracSpinor::Index	w,
		DiracSpinor::Index	x,
		DiracSpinor::Index	y,
		int	twoJ
	)

Antisymmetrised two-body Coulomb matrix element in the coupled CSF basis.

Evaluates the angular-reduced, antisymmetrised Coulomb interaction between two two-electron CSFs \( |vw; J\rangle \) and \( |xy; J\rangle \):

\[ \langle vw; J \| g \| xy; J \rangle = \eta_{vw}\eta_{xy} \sum_k (-1)^{j_v+j_x+k+J} \begin{Bmatrix} j_v & j_w & J \\ j_y & j_x & k \end{Bmatrix} Q^k_{vwxy} + \text{exchange}, \]

where \( \eta_{ab} = 1/\sqrt{2} \) if \( a = b \) (identical-particle normalisation) and 1 otherwise, and \( Q^k \) are the Coulomb integrals stored in qk.

Parameters

qk	Table of Coulomb \( Q^k \) integrals.
v,w	Indices of the bra single-particle states.
x,y	Indices of the ket single-particle states.
twoJ	Twice the total angular momentum 2J of the coupled pair.

Returns

Antisymmetrised, angular-reduced two-body Coulomb matrix element.

CSF2_Sigma2()

double CI::CSF2_Sigma2	(	const Coulomb::LkTable &	Sk,
		DiracSpinor::Index	v,
		DiracSpinor::Index	w,
		DiracSpinor::Index	x,
		DiracSpinor::Index	y,
		int	twoJ,
		const Coulomb::QkTable *	qk = nullptr,
		const std::vector< double > &	hk = {},
		const Coulomb::LkTable *	dSk = nullptr,
		double	dE0 = 0.0
	)

Two-body \( \Sigma_2 \) (MBPT) correction to CSF2_Coulomb().

Evaluates the same angular reduction as CSF2_Coulomb(), but using the two-body \( \Sigma_2 \) integrals \( S^k \) stored in Sk in place of the Coulomb \( Q^k \) integrals. Adds the second-order MBPT correction to the two-electron interaction.

Parameters

Sk	Table of two-body \( \Sigma_2 \) ( \( L^k \)) integrals.
v,w	Indices of the bra single-particle states.
x,y	Indices of the ket single-particle states.
twoJ	Twice the total angular momentum 2J of the coupled pair.

Returns

Antisymmetrised two-body \( \Sigma_2 \) matrix element.

CSF2_Breit()

double CI::CSF2_Breit	(	const Coulomb::WkTable &	Bk,
		DiracSpinor::Index	v,
		DiracSpinor::Index	w,
		DiracSpinor::Index	x,
		DiracSpinor::Index	y,
		int	twoJ
	)

Antisymmetrised two-body Breit matrix element in the coupled CSF basis.

Evaluates the same angular reduction as CSF2_Coulomb(), but using the Breit \( B^k \) integrals stored in Bk.

Parameters

Bk	Table of Breit \( W^k \) integrals.
v,w	Indices of the bra single-particle states.
x,y	Indices of the ket single-particle states.
twoJ	Twice the total angular momentum 2J of the coupled pair.

Returns

Antisymmetrised two-body Breit matrix element.

Sigma2_AB()

double CI::Sigma2_AB	(	const CI::CSF2 &	A,
		const CI::CSF2 &	B,
		int	twoJ,
		const Coulomb::LkTable &	Sk,
		const Coulomb::QkTable *	qk = nullptr,
		const std::vector< double > &	hk = {},
		const Coulomb::LkTable *	dSk = nullptr,
		double	dE0 = 0.0
	)

Two-body \( \Sigma_2 \) correction to Hab().

Evaluates the MBPT \( \Sigma_2 \) contribution to the CI matrix element using CSF2_Sigma2(). Add to Hab() to form the full CI+MBPT Hamiltonian matrix element.

Parameters

A,B	The two CSFs.
twoJ	Twice the total angular momentum 2J.
Sk	Table of \( \Sigma_2 \) ( \( L^k \)) integrals.

Returns

\( \Sigma_2 \) correction to \( H_{AB} \).

Breit_AB()

double CI::Breit_AB	(	const CI::CSF2 &	A,
		const CI::CSF2 &	B,
		int	twoJ,
		const Coulomb::WkTable &	Bk
	)

Breit correction to Hab().

Evaluates the two-body Breit contribution to the CI matrix element using CSF2_Breit(). Add to Hab() to include the Breit interaction.

Parameters

A,B	The two CSFs.
twoJ	Twice the total angular momentum 2J.
Bk	Table of Breit \( W^k \) integrals.

Returns

Breit correction to \( H_{AB} \).

corrected_Sigma()

double CI::corrected_Sigma	(	double	Sigma,
		double	dSigma,
		double	dE
	)

Resummed derivative (dSigma/dE) correction to a Sigma_1 matrix element (Kozlov formula).

Returns the corrected matrix element,

\[ \Sigma \to \Sigma \left[ 1 - \delta E \, (d\Sigma/dE)/\Sigma \right]^{-1}, \]

which resums the linear expansion \( \Sigma(\epsilon_0 + \delta E) \approx \Sigma + \delta E \, d\Sigma/dE \).

Guard: if the corrected value exceeds \( |\Sigma| \), then \( \delta E \, (d\Sigma/dE) \) is approaching \( \Sigma \) - the pole of the resummed (Pade) form, at \( \delta E = \Sigma / (d\Sigma/dE) \) - where the expression diverges; the correction is distrusted, and the uncorrected \( \Sigma \) is returned.

Parameters

Sigma	Matrix element \( \langle a|\Sigma_1|b\rangle \).
dSigma	Energy derivative, \( \langle a|d\Sigma_1/dE|b\rangle \).
dE	Energy shift \( \delta E \) from the energy Sigma_1 was evaluated at.

Returns

Corrected matrix element.

corrected_Sk()

double CI::corrected_Sk	(	double	Sk,
		double	dSk,
		double	dE0
	)

Resummed shift of a \( \Sigma_2 \) integral to a new reference energy E0 (Brillouin-Wigner denominators).

\( S^k \to (S^k)^2 / (S^k - \delta E_0\, dS^k/dE_0) \), the same resummation as corrected_Sigma. Exact for a single energy denominator (each term is \( N/(D + \delta E_0) \)), so it holds up much better than the linear expansion when \( \delta E_0 \) is a sizeable fraction of the denominator.

Unlike corrected_Sigma there is no guard against the correction enhancing \( |S^k| \): for \( \Sigma_2 \) the shift legitimately goes either way. The only guard is on crossing the pole of the resummed form, at \( \delta E_0 = S^k / (dS^k/dE_0) \) (detected by the denominator changing sign), beyond which the expansion is meaningless: the unshifted \( S^k \) is returned there.

Parameters

Sk	Integral \( S^k \), evaluated at the reference E0.
dSk	Derivative \( dS^k/dE_0 \), at the same reference.
dE0	Shift from the reference, \( E_0 - E_0^{\rm ref} \).

Returns

Shifted \( S^k \).

calculate_dSdE_correction() [1/2]

Sigma1Correction CI::calculate_dSdE_correction	(	const std::vector< DiracSpinor > &	ci_basis,
		const std::vector< DiracSpinor > &	s1_basis_core,
		const std::vector< DiracSpinor > &	s1_basis_excited,
		const Coulomb::QkTable &	qk
	)

Builds the Sigma_1 derivative-correction tables; see Sigma1Correction.

For each same-kappa pair in ci_basis, computes the (uncorrected) Sigma_1 matrix element and its energy derivative (central finite difference of MBPT::Sigma_vw). Sigma_1 is evaluated at the energy of the first state of each kappa in ci_basis - the same convention as calculate_h1_table(), so the stored S1 matches the Sigma_1 included in the h1 table.

Parameters

ci_basis	Basis states for which table entries are needed.
s1_basis_core	Core states used as internal lines for Sigma_1.
s1_basis_excited	Excited states used as internal lines for Sigma_1.
qk	Table of Coulomb \( Q^k \) integrals.

Returns

Filled Sigma1Correction tables (E0 left unset; see iterate_E0).

calculate_dSdE_correction() [2/2]

Sigma1Correction CI::calculate_dSdE_correction	(	const std::vector< DiracSpinor > &	ci_basis,
		const MBPT::CorrelationPotential &	Sigma
	)

Builds the Sigma_1 derivative-correction tables directly from a correlation potential that holds dSigma/dE matrices.

Fills S1 from the actual Sigma of the correlation potential (via CorrelationPotential::SigmaFv - so it matches the Sigma_1 in the h1 table exactly, whatever the method: Goldstone, Feynman, all-orders), and dS1 from its stored dSigma/dE matrices (see the Correlations option derivative). Much faster than the qk-table overload (matrix applications only), and consistent with the all-orders Sigma.

The reference energies e_sigma are taken from the stored Sigma data (the energy each Sigma matrix was formed at, first entry of each kappa).

Parameters

ci_basis	Basis states for which table entries are needed.
Sigma	Correlation potential; must hold dSigma/dE matrices (see CorrelationPotential::has_derivative()).

Returns

Filled Sigma1Correction tables (E0 left unset; see iterate_E0).

Note

The ladder part (Sigma_L) is included in S1 (via SigmaFv) but has no energy derivative, so it is absent from dS1.

iterate_E0()

LinAlg::Matrix< double > CI::iterate_E0	(	PsiJPi *	psi,
		const Sigma1Correction &	s1c,
		const Coulomb::meTable< double > &	h1,
		const Coulomb::QkTable &	qk,
		const Coulomb::WkTable *	Bk = nullptr,
		const Coulomb::LkTable *	Sk = nullptr,
		const std::vector< double > &	hk = {},
		const Coulomb::LkTable *	dSk = nullptr,
		double	E0_sigma2 = 0.0,
		std::ostream &	outstream = std::cout
	)

Finds the reference energy E0 for the dSigma/dE correction, for a single (J, parity); returns the CI Hamiltonian built with it.

E0 is the lowest energy of this (J, parity), which is only known once we have solved - so it is found self-consistently. If s1c has no E0 set (i.e., 0.0), the iteration starts from the lowest zeroth-order configuration energy. The first pass builds the CI Hamiltonian at the current E0 and diagonalises it for the lowest level. E0 enters only through the (small) dSigma/dE correction, so the state itself hardly changes from pass to pass: later passes rebuild the Hamiltonian with the updated E0 and take the new E0 as the expectation value of the new Hamiltonian in the (unchanged) state - first-order perturbation theory in the change to the Hamiltonian. Only the first pass is diagonalised.

Parameters

psi	Solved for the lowest level (updated in place).
s1c	Correction tables. E0 belongs to a single (J, parity), while these tables are shared between them, so a local copy is made: s1c is left as it was.
h1	One-body matrix element table (includes Sigma_1).
qk	Coulomb \( Q^k \) table.
Bk	Pointer to Breit table; ignored if nullptr.
Sk	Pointer to \( \Sigma_2 \) table; ignored if nullptr.
hk	Average S^k/Q^k ratios; see MBPT::average_hk.
dSk	Pointer to the dS^k/dE0 table (Brillouin-Wigner Sigma_2); ignored if nullptr. Sigma_2 is shifted to the current E0 each pass, so both corrections move together as E0 converges.
E0_sigma2	Reference E0 that Sk and dSk were tabulated at.
outstream	Stream for the per-pass output.

Returns

CI Hamiltonian matrix, built with the converged E0.

Hab()

double CI::Hab	(	const CI::CSF2 &	A,
		const CI::CSF2 &	B,
		int	twoJ,
		const Coulomb::meTable< double > &	h1,
		const Coulomb::QkTable &	qk,
		const Sigma1Correction *	s1c = nullptr
	)

CI Hamiltonian matrix element between two two-electron CSFs.

Computes \( H_{AB} = \langle A | \hat{H} | B \rangle \) using the Slater-Condon rules, including one-body terms from h1 (which may already incorporate \( \Sigma_1 \) corrections) and the two-body Coulomb interaction via CSF2_Coulomb().

Does NOT include \( \Sigma_2 \) or Breit corrections; add those via Sigma2_AB() and Breit_AB() respectively.

Parameters

A,B	The two CSFs.
twoJ	Twice the total angular momentum 2J.
h1	Table of one-body matrix elements \( \langle a | h_1 | b \rangle \).
qk	Table of Coulomb \( Q^k \) integrals.
s1c	Optional derivative (dSigma/dE) correction to Sigma_1; applied to each one-body matrix element (with the spectator orbital energy) if given. See Sigma1Correction.

Returns

CI Hamiltonian matrix element \( H_{AB} \).

calculate_h1_table() [1/2]

Coulomb::meTable< double > CI::calculate_h1_table	(	const std::vector< DiracSpinor > &	ci_basis,
		const std::vector< DiracSpinor > &	s1_basis_core,
		const std::vector< DiracSpinor > &	s1_basis_excited,
		const Coulomb::QkTable &	qk,
		bool	include_Sigma1
	)

Builds the one-body Hamiltonian matrix element table for the CI basis.

Constructs a lookup table of single-particle matrix elements \( \langle a | h_1 | b \rangle \) for all pairs \( a, b \) in ci_basis. The diagonal elements are the HF single-particle energies.

If include_Sigma1 is true, the one-body MBPT \( \Sigma_1 \) correction is computed from the Coulomb integrals in qk using s1_basis_core and s1_basis_excited as the internal lines of the MBPT diagrams and added to the diagonal.

Parameters

ci_basis	Basis states for which table entries are needed.
s1_basis_core	Core states used as internal lines for \( \Sigma_1 \).
s1_basis_excited	Excited states used as internal lines for \( \Sigma_1 \).
qk	Table of Coulomb \( Q^k \) integrals.
include_Sigma1	If true, add one-body MBPT \( \Sigma_1 \) corrections.

Returns

Table of \( \langle a | h_1 | b \rangle \) matrix elements.

Warning

Assumes ci_basis states are Hartree-Fock eigenstates, so off-diagonal HF terms vanish.

calculate_h1_table() [2/2]

Coulomb::meTable< double > CI::calculate_h1_table	(	const std::vector< DiracSpinor > &	ci_basis,
		const MBPT::CorrelationPotential &	Sigma,
		bool	include_Sigma1
	)

Builds the one-body Hamiltonian table using a precomputed CorrelationPotential.

Overload of calculate_h1_table() that uses a CorrelationPotential object (i.e., a precomputed \( \Sigma_1 \) operator) instead of computing MBPT diagrams on the fly. Preferred when a CorrelationPotential is available, as it is generally faster and more complete.

Parameters

ci_basis	Basis states for which table entries are needed.
Sigma	Precomputed one-body correlation potential \( \Sigma_1 \).
include_Sigma1	If true, include \( \Sigma_1 \) corrections from Sigma.

Returns

Table of \( \langle a | h_1 | b \rangle \) matrix elements.

calculate_Bk()

Coulomb::WkTable CI::calculate_Bk	(	const std::string &	bk_filename,
		const HF::Breit *const	pBr,
		const std::vector< DiracSpinor > &	ci_basis,
		int	max_k,
		bool	no_new_integralsQ = false
	)

Builds or loads the two-body Breit integral table.

Computes Breit \( W^k \) integrals for all pairs in ci_basis using the Breit operator pBr. Results are cached to/from bk_filename.

If pBr is nullptr or no_new_integralsQ is true, no new integrals are computed; only cached values are loaded.

Parameters

bk_filename	Filename for caching the \( W^k \) table.
pBr	Pointer to Breit operator; if nullptr, returns empty table.
ci_basis	Basis for which Breit integrals are needed.
max_k	Maximum multipolarity k to include.
no_new_integralsQ	If true, skip computing any new integrals.

Returns

Table of Breit \( W^k \) integrals.

basis_subset()

std::vector< DiracSpinor > CI::basis_subset	(	const std::vector< DiracSpinor > &	basis,
		const std::string &	include_str,
		const std::string &	exclude_str = ""
	)

Returns the subset of basis matching include_str, excluding states in exclude_str.

Filters basis to retain only states described by the ampsci basis-string notation (e.g., "20spdf") that are not part of the frozen core.

Parameters

basis	Full single-particle basis to filter.
include_str	Basis-string specifying which states to keep; if empty, all states in basis are kept (subject to the frozen-core exclusion).
exclude_str	Basis-string specifying core states to exclude.

Returns

Filtered basis vector.

ReducedME() [1/2]

double CI::ReducedME	(	const LinAlg::View< const double > &	cA,
		const std::vector< CI::CSF2 > &	CSFAs,
		int	twoJA,
		const LinAlg::View< const double > &	cB,
		const std::vector< CI::CSF2 > &	CSFBs,
		int	twoJB,
		const Coulomb::meTable< double > &	h,
		int	K_rank,
		int	Parity
	)

Reduced matrix element between two CI states (low-level overload).

Evaluates the reduced matrix element of a rank-K_rank tensor operator between two CI states:

\[ \redmatel{A}{T^K}{B} = \sum_{ij} c_i^A \, c_j^B \, \redmatel{\text{CSF}_i}{T^K}{\text{CSF}_j}, \]

where the single-particle reduced matrix elements are looked up from h.

Parameters

cA,cB	CI expansion coefficient vectors for states A and B.
CSFAs,CSFBs	CSF bases for states A and B respectively.
twoJA,twoJB	Twice the total angular momentum of states A and B.
h	Lookup table of single-particle reduced matrix elements.
K_rank	Rank of the tensor operator.
Parity	Parity of the operator (+1 or -1).

Returns

Reduced matrix element \( \redmatel{A}{T^K}{B} \).

RME_CSF2()

double CI::RME_CSF2	(	const CI::CSF2 &	X,
		int	twoJX,
		const CI::CSF2 &	V,
		int	twoJV,
		const Coulomb::meTable< double > &	h,
		int	K_rank
	)

Reduced matrix element between two two-electron CSFs.

Evaluates \( \redmatel{X; J_X}{T^K}{V; J_V} \) for a rank-K_rank one-body tensor operator using the standard 6j angular reduction, accounting for identical-particle normalisation factors.

Warning

This function may not handle all cases correctly; results should be verified for non-trivial configurations.

leading_config()

std::pair< std::string, double > CI::leading_config	(	const LinAlg::View< const double > &	coefs,
		const std::vector< CSF2 > &	csfs
	)

Leading non-relativistic configuration of a CI state.

Sums |c|^2 over the relativistic CSFs belonging to each non-relativistic configuration, and returns the configuration with the largest total weight.

Parameters

coefs	CI expansion coefficients (one per CSF).
csfs	The CSF basis (matching coefs).

Returns

Pair {configuration label (non-rel notation), total |c|^2 weight}.

Term_S_L()

std::pair< int, int > CI::Term_S_L	(	int	l1,
		int	l2,
		int	twoJ,
		double	gJ_target
	)

Determines the best-fit (S, L) term for a two-electron state by matching the g-factor.

Iterates over all allowed (S, L) combinations for given orbital angular momenta l1, l2 and total twoJ /2, and returns the pair whose Lande g-factor is closest to gJ_target.

Parameters

l1,l2	Orbital angular momenta of the two electrons.
twoJ	Twice the total angular momentum 2J.
gJ_target	Target g-factor to match.

Returns

Best-fit {2S, L} pair.

Term_S_L_from_expectation()

std::pair< int, int > CI::Term_S_L_from_expectation	(	double	L2,
		double	S2,
		int	twoJ
	)

Determines the (S, L) term for a two-electron state from the expectation values of L^2 and S^2.

Returns the (S, L) pair, subject to the triangle condition with J, that minimises the combined distance |L(L+1) - L2| + |S(S+1) - S2|. For a mixed state, this is the nearest (dominant) term; the purity is indicated by L2, S2 themselves. See expectation_L2S2.

Parameters

L2	Expectation value of L^2.
S2	Expectation value of S^2.
twoJ	Twice the total angular momentum 2J.

Returns

Best-fit {S, L} pair.

Term_Symbol() [1/2]

std::string CI::Term_Symbol	(	int	two_J,
		int	L,
		int	two_S,
		int	parity
	)

Returns spectroscopic term symbol string, e.g. "3P_1".

Term_Symbol() [2/2]

std::string CI::Term_Symbol	(	int	L,
		int	two_S,
		int	parity
	)

Returns term symbol without the J subscript, e.g. "3P".

construct_Hci() [1/2]

LinAlg::Matrix< double > CI::construct_Hci	(	const PsiJPi &	psi,
		const Coulomb::meTable< double > &	h1,
		const Coulomb::QkTable &	qk,
		const Coulomb::WkTable *	Bk = nullptr,
		const Coulomb::LkTable *	Sk = nullptr,
		const Sigma1Correction *	s1c = nullptr,
		const std::vector< double > &	hk = {},
		const Coulomb::LkTable *	dSk = nullptr,
		double	dE0 = 0.0
	)

Constructs the full CI Hamiltonian matrix in the CSF basis.

Builds the symmetric matrix \( H_{AB} \) for all CSF pairs in psi, calling Hab() for each element and optionally adding Breit and \( \Sigma_2 \) corrections.

Parameters

psi	CI solution container holding the CSF basis and J/parity.
h1	One-body matrix element table (may include \( \Sigma_1 \)).
qk	Coulomb \( Q^k \) integral table.
Bk	Pointer to Breit \( W^k \) table; ignored if nullptr.
Sk	Pointer to \( \Sigma_2 \) \( L^k \) table; ignored if nullptr.
s1c	Pointer to derivative (dSigma/dE) correction for Sigma_1; ignored if nullptr. See Sigma1Correction.
hk	Average S^k/Q^k ratios; if non-empty, diagrams with no stored S^k use S^k = hk[k]*Q^k. See MBPT::average_hk.
dSk	Pointer to the dS^k/dE0 table; ignored if nullptr. If given, each stored S^k is shifted to the target E0. See CI::corrected_Sk.
dE0	Shift of the target E0 from the one Sk was tabulated at.

Returns

Full CI Hamiltonian matrix in the CSF basis.

construct_Hci() [2/2]

LinAlg::Matrix< double > CI::construct_Hci	(	const PsiJPi &	psi,
		const Integrals &	ints
	)

Constructs the CI Hamiltonian matrix from a set of integral tables.

Overload of construct_Hci() taking the tables as Integrals; the Breit and \( \Sigma_2 \) corrections are included if those tables are non-empty.

If the iterative (dSigma/dE) correction is included, the reference energy E0 (which also sets the Sigma_2 BW shift, if present) belongs to this (J, parity) block: the lowest solved energy of psi is used if it has solutions; otherwise the stored fallback is used, which is the ground-state energy of the CI run. Blocks never solved directly - e.g., the target block of a mixed state - thus have their corrections evaluated at the run's ground-state energy; the difference is higher order.

Parameters

psi	CI solution container holding the CSF basis and J/parity.
ints	Integral tables, e.g., from Wavefunction::CI_integrals().

Returns

Full CI Hamiltonian matrix in the CSF basis.

ReducedME() [2/2]

double CI::ReducedME ( const PsiJPi &  As, std::size_t  iA, const PsiJPi &  Bs, std::size_t  iB, const Coulomb::meTable< double > &  h, int  K_rank, int  Parity  )

inline

Reduced matrix element between two CI states (PsiJPi overload).

Convenience wrapper around the low-level ReducedME() overload. Extracts expansion coefficients and CSF lists from As and Bs for the requested solution indices iA and iB.

Parameters

As,Bs	CI solution containers for the two states.
iA,iB	Solution indices within As and Bs.
h	Lookup table of single-particle reduced matrix elements.
K_rank	Rank of the tensor operator.
Parity	Parity of the operator (+1 or -1).

Returns

Reduced matrix element \( \redmatel{A}{T^K}{B} \).

configuration_interaction()

Solutions CI::configuration_interaction	(	const IO::InputBlock &	input,
		const Wavefunction &	wf
	)

Runs Configuration Interation: returns CI solutions for all requested J and parity values.

Reads options from input, the CI Input Block, and constructs the CI basis from wf, computes the required Coulomb (and optionally Breit and two-body MBPT) integrals, then calls run_CI() for each requested (J, parity) pair.

The returned vector contains one PsiJPi per {J, parity} combination, each holding the eigenvalues and CI expansion coefficients for the requested number of solutions.

} Always check for up-to-date options from command line: $ ampsci -i CI See also run_CI, which this function calls

Parameters

input	Input block containing CI options.
wf	Fully initialised Wavefunction object supplying the orbital basis and radial grid.

Returns

Solutions: one PsiJPi per (J, parity) combination requested, together with the integral tables used to build the CI Hamiltonians.

Note

If run with read_only, no integrals are calculated, and the returned integral tables are empty (see Integrals::availableQ()).

run_CI()

PsiJPi CI::run_CI	(	const std::vector< DiracSpinor > &	ci_sp_basis,
		int	twoJ,
		int	parity,
		int	num_solutions,
		std::optional< double >	all_below_cm,
		const Coulomb::meTable< double > &	h1,
		const Coulomb::QkTable &	qk,
		const Coulomb::WkTable &	Bk,
		const Coulomb::LkTable &	Sk,
		bool	include_Sigma2,
		bool	print_details,
		bool	read_only = false,
		std::ostream &	outstream = std::cout,
		const std::string &	ci_fname = "",
		const Sigma1Correction *	s1c = nullptr,
		const std::vector< double > &	hk = {},
		const Coulomb::LkTable *	dSk = nullptr,
		double	E0_sigma2 = 0.0
	)

Constructs and solves the CI eigenvalue problem for a single J,pi.

Builds the CI+MBPT Hamiltonian matrix in the basis of two-electron configuration state functions (CSFs) with total angular momentum twoJ /2 and parity parity, then solves the eigenvalue problem to obtain CI energies and expansion coefficients.

The Hamiltonian includes:

• one-body terms from h1
• two-body Coulomb interaction via the \( Q^k \) integrals in qk
• optionally, two-body Breit corrections via \( B^k \) integrals in Bk (used if Bk is non-empty)
• optionally, two-body MBPT \( \Sigma_2 \) corrections via \( S^k \) integrals in Sk (used if include_Sigma2 is true, and Sk is non-empty)

The number of solutions returned is controlled by num_solutions and all_below: if all_below is set it takes precedence and all eigenstates with total energy below the threshold are found.

Parameters

ci_sp_basis	Single-particle basis states spanning the CI space.
twoJ	Twice the total angular momentum, 2J (must be a non-negative even integer for two-electron systems).
parity	Parity of the block: +1 (even) or -1 (odd).
num_solutions	Number of lowest eigenstates to find. Ignored if all_below_cm is set. Pass 0 to find all solutions.
all_below_cm	If set, find all eigenstates with total energy below this value (in cm^-1). Overrides num_solutions.
h1	Table of one-body Hamiltonian matrix elements between single-particle basis states.
qk	Table of two-body Coulomb \( Q^k \) integrals.
Bk	Table of two-body Breit \( B^k \) integrals. Ignored (treated as absent) if the table is empty.
Sk	Table of two-body MBPT \( \Sigma_2 \) ( \( S^k \)) integrals. Only used when include_Sigma2 is true.
include_Sigma2	If true, add two-body MBPT corrections from Sk to the CI Hamiltonian.
print_details	If true, print a breakdown of the leading configurations for each solution. Leads to very large output if num_solutions is large
read_only	If true, only the solutions already in the file are read; the eigenvalue problem is not solved, and nothing is written [default: false].
outstream	Output stream for progress and results [default: stdout].
ci_fname	Filename for reading/writing CI solutions ("" disables).
s1c	Pointer to derivative (dSigma/dE) correction for Sigma_1; ignored if nullptr. See Sigma1Correction. Its reference energy E0 is the lowest energy of this (J, parity), so is found self-consistently: see iterate_E0.
hk	Average S^k/Q^k ratios, for extrapolating Sigma_2 beyond the cis2 basis; empty for no extrapolation. See MBPT::average_hk.
dSk	Pointer to the dS^k/dE0 table (Brillouin-Wigner Sigma_2); ignored if nullptr. See CI::corrected_Sk.
E0_sigma2	Reference E0 that Sk and dSk were tabulated at.

Returns

PsiJPi (PsiJPi) containing the CI eigenvalues and expansion coefficients for the requested solutions.

parse_level()

std::optional< Level > CI::parse_level

(

std::string_view

str

)

Parses the text form of a CI level reference; see Level.

Accepts 2+:3 (standard) and e2:3; the index is optional. Surrounding whitespace is ignored. Only integer J is accepted, since only two-electron CI is implemented.

Parameters

str	Text form, e.g., 2+:3, e2:3, 0+.

Returns

The level; empty if str is not a valid level reference.

to_string()

std::string CI::to_string

(

const Level &

level

)

Text form of a CI level reference, e.g., "2+:3"; see Level.

form_CSFs()

std::vector< CSF2 > CI::form_CSFs	(	int	twoJ,
		int	parity,
		const std::vector< DiracSpinor > &	cisp_basis
	)

Forms all two-electron CSFs with given total J and parity.

Iterates over all pairs of single-particle states in cisp_basis and retains those whose angular momenta can be coupled to total \( J = \) twoJ /2 and whose combined parity equals parity. Duplicate pairs are excluded by construction.

Parameters

twoJ	Twice the total angular momentum 2J.
parity	Total parity: +1 (even) or -1 (odd).
cisp_basis	Single-particle basis from which CSFs are constructed.

Returns

Sorted list of all valid two-electron CSFs for the given J and parity.

LS_amplitude()

double CI::LS_amplitude	(	int	n1,
		int	l1,
		int	twoj1,
		int	n2,
		int	l2,
		int	twoj2,
		int	L,
		int	S,
		int	twoJ
	)

jj -> LS recoupling amplitude for an antisymmetrised two-electron CSF.

Returns the amplitude of the antisymmetrised jj-coupled CSF \( |\{(n_1 l_1 j_1)(n_2 l_2 j_2)\}; J\rangle \) (orbitals in stored, i.e., sorted, order) onto the antisymmetrised LS-coupled state \( |\{(n_1 l_1)(n_2 l_2)\} L S; J\rangle \) of the same non-relativistic configuration:

\[ A(L,S) = \eta \sqrt{[j_1][j_2][L][S]} \begin{Bmatrix} l_1 & l_2 & L \\ 1/2 & 1/2 & S \\ j_1 & j_2 & J \end{Bmatrix} \]

Taken in the non-relativistic limit: the radial orbitals of \( j = l \pm 1/2 \) are treated as identical (overlap = 1).

For a common non-relativistic shell ( \( n_1 = n_2 \), \( l_1 = l_2 \)) only L+S even terms exist (Pauli). When additionally \( j_1 \neq j_2 \) the L+S odd components cancel in the antisymmetrisation and the even ones carry \( \eta = \sqrt{2} \); otherwise \( \eta = 1 \). In all cases \( \sum_{LS} A^2 = 1 \).

Parameters

n1,l1,twoj1	Quantum numbers of the first stored orbital.
n2,l2,twoj2	Quantum numbers of the second stored orbital.
L,S	Total orbital and spin angular momenta of the LS term.
twoJ	Twice the total angular momentum 2J.

Returns

Recoupling amplitude A(L,S); zero if forbidden.

Note

The sign convention follows the stored (sorted) orbital order; since nk_index sorting keeps the (n, l) order identical for all CSFs of one non-relativistic configuration, relative signs between such CSFs are consistent.

expectation_L2S2()

std::pair< double, double > CI::expectation_L2S2	(	const LinAlg::View< const double > &	coefs,
		const std::vector< CSF2 > &	csfs,
		int	twoJ
	)

Expectation values of L^2 and S^2 for a two-electron CI state.

Recouples each CSF to LS coupling (see LS_amplitude) and accumulates, per non-relativistic configuration g, \( B_g(L,S) = \sum_{I \in g} c_I A_I(L,S) \), giving

\[ \langle L^2 \rangle = \sum_{g,L,S} B_g(L,S)^2 \, L(L+1), \qquad \langle S^2 \rangle = \sum_{g,L,S} B_g(L,S)^2 \, S(S+1). \]

These are expectation values of the CI state, not eigenvalues: deviation from L(L+1), S(S+1) measures the LS-purity of the state.

Parameters

coefs	CI expansion coefficients (one per CSF).
csfs	The CSF basis (matching coefs).
twoJ	Twice the total angular momentum 2J.

Returns

Pair {<L^2>, <S^2>}.

Note

Non-relativistic limit: radial overlaps between j = l +- 1/2 orbitals are set to 1, so for a normalised state the total LS weight is exactly 1 and no renormalisation is required.

TPsi_reduced()

LinAlg::Vector< double > CI::TPsi_reduced	(	const std::vector< CSF2 > &	CSFs,
		int	twoJ,
		const PsiJPi &	Psi0,
		std::size_t	i0,
		const Coulomb::meTable< double > &	h,
		int	K_rank
	)

Action of a one-body operator on a CI state, in the CSF basis.

Forms the vector

\[ T_I = \sum_K \redmatel{I}{T^{(K)}}{K} \, c^{(0)}_K, \]

where \( I \) runs over the CSFs in CSFs (which have total angular momentum twoJ /2), \( K \) runs over the CSFs of the reference state \( \Psi_0 \) (solution i0 of Psi0), and \( c^{(0)}_K \) are its CI expansion coefficients.

The CSF matrix elements are reduced (see RME_CSF2), formed from the single-particle reduced matrix elements in h; so is the result.

Parameters

CSFs	CSFs spanning the block the operator maps into.
twoJ	Twice the total angular momentum, 2J, of CSFs.
Psi0	CI solutions containing the reference state.
i0	Index of the reference solution within Psi0.
h	Table of single-particle reduced matrix elements of T.
K_rank	Rank of the tensor operator T.

Returns

Vector \( T_I \), of length CSFs .size().

solve_mixed_state() [1/2]

PsiJPi CI::solve_mixed_state	(	const PsiJPi &	Psi0,
		std::size_t	i0,
		const PsiJPi &	target,
		const LinAlg::Matrix< double > &	Hci,
		const Coulomb::meTable< double > &	h,
		int	K_rank,
		double	omega = 0.0
	)

Solves the CI mixed-states (Sternheimer) equation for a one-body operator.

Finds the first-order correction to the CI state \( \Psi_0 \) (solution i0 of Psi0, with energy \( E_0 \)) due to the one-body operator \( T^{(K)} \), expanded over the CSFs of a single (J, parity) block:

\[ \ket{\delta\Psi} = \sum_I c_I \ket{I; J^\pi}. \]

The coefficients solve the linear system

\[ \sum_J \left[ \matel{I}{H}{J} - (E_0 + \omega) \, \delta_{IJ} \right] c_J = - \sum_K \redmatel{I}{T^{(K)}}{K} \, c^{(0)}_K, \]

where \( H \) is the CI Hamiltonian in the block defined by target (given as the matrix Hci, e.g., from construct_Hci), and the right-hand side is formed by TPsi_reduced.

Since the right-hand side is reduced (in \( T \)), so is the solution: for any CI state \( A \) in the same block, the mixed state satisfies

\[ \sum_I c^A_I \, c_I = \frac{\redmatel{A}{T^{(K)}}{\Psi_0}}{E_0 + \omega - E_A}, \]

i.e., the sum over the entire spectrum of that block, without finding (or summing over) the individual CI solutions.

Returned as a PsiJPi holding a single "solution", the coefficients \( c_I \). Its stored energy is \( E_0 \), that of the reference state (not an eigenvalue).

Parameters

Psi0	CI solutions containing the reference state.
i0	Index of the reference solution within Psi0.
target	Defines the block the mixed state lives in (2J, parity, and CSF list); its solutions, if any, are not used.
Hci	CI Hamiltonian matrix in the CSF basis of target.
h	Table of single-particle reduced matrix elements of T.
K_rank	Rank of the tensor operator T.
omega	Frequency: the mixed state due to a time-dependent operator, \( T e^{-i\omega t} \), has denominators \( E_0 + \omega - E_A \) [0].

Returns

PsiJPi for the target block, holding the single mixed state.

See also

project_out, to remove individual levels from the mixed state.

Note

If target has the same J and parity as Psi0, and \( \omega = 0 \), the matrix on the left is singular, since \( \Psi_0 \) itself has zero eigenvalue. In that case, \( \Psi_0 \) is projected out (equivalent to subtracting \( \redmatel{\Psi_0}{T}{\Psi_0} \) from the right-hand side).

Any other state degenerate with \( E_0 + \omega \) also makes the system singular. Its term in the sum over states is divergent, and must be dealt with separately, as in degenerate perturbation theory.

If the operator cannot connect the two blocks (triangle rule or parity), the right-hand side vanishes, and the mixed state is zero.

project_out()

PsiJPi CI::project_out	(	PsiJPi	dPsi,
		const PsiJPi &	levels,
		const std::vector< std::size_t > &	indices
	)

Removes CI levels from a mixed state, so that it is orthogonal to them.

A mixed state is implicitly a sum over the entire spectrum of its (J, parity):

\[ \ket{\delta\Psi} = \sum_A \ket{A} \frac{\redmatel{A}{T^{(K)}}{\Psi_0}}{E_0 + \omega - E_A}. \]

Subtracting the projection onto the listed levels removes exactly their terms from that sum, so they may be treated separately: e.g., with experimental energies or matrix elements.

Parameters

dPsi	Mixed state, from solve_mixed_state (taken by value).
levels	Solved CI levels of the same (J, parity): the eigenstates of the CI Hamiltonian used for the mixed state.
indices	Which solutions of levels to remove.

Returns

The mixed state, orthogonal to the listed levels.

Note

Removing a level degenerate with \( E_0 + \omega \) does not help: the mixed state itself does not exist in that case (see solve_mixed_state).

solve_mixed_state() [2/2]

PsiJPi CI::solve_mixed_state	(	const PsiJPi &	Psi0,
		std::size_t	i0,
		int	twoJ,
		int	parity,
		const std::vector< DiracSpinor > &	ci_sp_basis,
		const Coulomb::meTable< double > &	h,
		int	K_rank,
		const Coulomb::meTable< double > &	h1,
		const Coulomb::QkTable &	qk,
		const Coulomb::WkTable *	Bk = nullptr,
		const Coulomb::LkTable *	Sk = nullptr,
		double	omega = 0.0
	)

Solves the CI mixed-states equation; constructs the CI matrix internally.

Convenience overload of solve_mixed_state: forms the CSFs for the requested (twoJ, parity) block from ci_sp_basis, constructs the CI Hamiltonian matrix via construct_Hci, then solves the mixed-states equation.

Use the other overload if the CI matrix for the target block is already available (e.g., when several operators are considered).

Parameters

Psi0	CI solutions containing the reference state.
i0	Index of the reference solution within Psi0.
twoJ	Twice the total angular momentum, 2J, of the mixed state.
parity	Parity of the mixed state: +1 or -1. This is the parity of Psi0 times that of the operator.
ci_sp_basis	Single-particle basis used to construct the CSFs.
h	Table of single-particle reduced matrix elements of T.
K_rank	Rank of the tensor operator T.
h1	One-body matrix element table (may include Sigma_1).
qk	Coulomb Q^k integral table.
Bk	Pointer to Breit W^k table; ignored if nullptr.
Sk	Pointer to Sigma_2 L^k table; ignored if nullptr.
omega	Frequency; see the other overload [0].

Returns

PsiJPi for the (twoJ, parity) block, holding the single mixed state.

A_K_coefs()

std::pair< double, double > CI::A_K_coefs	(	int	K,
		int	kt,
		int	ks,
		int	twoJb,
		int	twoJn,
		int	twoJa
	)

Angular coefficients of the two terms of the second-order amplitude \( A^K \).

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}}{E_a + \omega_s - E_n} + \frac{\matel{b}{s_{q_2}}{n}\matel{n}{t_{q_1}}{a}}{E_a + \omega - E_n} \right], \]

the sum over \( n \) running over the magnetic quantum numbers too, and \( E_b = E_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 the reduced amplitude follows from the Wigner-Eckart theorem (the convention of DiracOperator::TensorOperator::rme3js):

\[ A^K_Q = (-1)^{J_b-m_b}\threej{J_b}{K}{J_a}{-m_b}{Q}{m_a}\,A^K . \]

\( A^K \) is the reduced matrix element of \( [t \times s]^{K} \) (both terms). In terms of reduced matrix elements of the two operators,

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

with the coefficients returned here

\[ c_1 = (-1)^{K}\sqrt{[K]}\,(-1)^{J_b+J_a} \sixj{K}{k_s}{k_t}{J_n}{J_b}{J_a}, \qquad c_2 = (-1)^{k_t+k_s}\sqrt{[K]}\,(-1)^{J_b+J_a} \sixj{K}{k_t}{k_s}{J_n}{J_b}{J_a}. \]

For a real transition the whole frequency is usually carried by \( t \), so that \( \omega = E_b - E_a \) and \( s \) is static. For the dynamic polarisability of a single state, \( b = a \) and \( \omega_s = -\omega \), giving the two denominators \( E_a \mp \omega - E_n \).

Sign convention

The coupling above is the standard Clebsch-Gordan one, \( t \) first. The alternative definition

\[ \tilde A^K_Q = (-1)^{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} = (-1)^K A^K_Q \]

differs by \( (-1)^K \). It cancels in anything rebuilt from \( A^K_Q \) (z_component), so only affects quantities taken directly from \( A^K \) at odd \( K \): the sign of \( \beta \).

Specific cases

With \( t = s = d \) (E1) and \( [J] \equiv 2J+1 \):

\[ \alpha_0 = \frac{A^0}{\sqrt{3[J_a]}} \quad (K=0,\ b=a), \qquad \alpha_2 = -\sqrt{\frac{2J(2J-1)}{3(J+1)(2J+1)(2J+3)}}\;A^2 \quad (K=2,\ J_b=J_a=J\ge1), \]

\[ \beta = \frac{A^1}{\sqrt{2}\,\redmatel{b}{\bm\sigma}{a}} \quad (K=1), \]

see sigma_rme. With \( t = d \) and \( s = h_W \) (PNC, \( k_s = 0 \), so \( K = 1 \)), at \( m_a = m_b = m \):

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

Parameters

K	Rank of the amplitude.
kt,ks	Ranks of the \( t \) and \( s \) operators.
twoJb	2J of the final state.
twoJn	2J of the intermediate states.
twoJa	2J of the initial state.

Returns

\( \{c_1, c_2\} \).

z_component()

double CI::z_component	(	int	K,
		int	kt,
		int	ks,
		int	twoJb,
		int	twoJa,
		int	two_m
	)

Converts the reduced amplitude \( A^K \) to its contribution to the z-component of the amplitude.

Undoing the coupling of A_K_coefs,

\[ A^{k_tk_s}_{q_1q_2} = \sum_{KQ}\braket{k_tq_1\,k_sq_2}{KQ}\,A^K_Q , \]

so the z-component ( \( m_a = m_b = m \), and \( q_1 = q_2 = 0 \), hence \( Q = 0 \)) is the sum over ranks

\[ A_{zz} \equiv A^{k_tk_s}_{00} = \sum_K \braket{k_t 0\,,\,k_s 0}{K 0} \, (-1)^{J_b-m}\threej{J_b}{K}{J_a}{-m}{0}{m} A^K, \]

the factor returned here being that of the \( K \) term. The Clebsch-Gordan coefficient is unity for \( k_s = 0 \) (as for a PNC amplitude), but not in general.

Parameters

K	Rank of the amplitude.
kt,ks	Ranks of the \( t \) and \( s \) operators.
twoJb,twoJa	2J of the final and initial states.
two_m	Twice the z-component of the angular momentum.

symm_sign()

int CI::symm_sign	(	const DiracOperator::TensorOperator *	h,
		int	twoJA,
		int	twoJB
	)

Relative sign between <A||h||B> and <B||h||A>, for CI states with total angular momenta 2J_A and 2J_B (cf DiracOperator::TensorOperator::symm_sign)

sigma_rme()

double CI::sigma_rme	(	const PsiJPi &	Psi_b,
		std::size_t	ib,
		const PsiJPi &	Psi_a,
		std::size_t	ia,
		const std::vector< DiracSpinor > &	ci_basis
	)

Reduced matrix element of the Pauli spin operator between two CI states, \( \redmatel{b}{\sigma}{a} \).

This is the factor that defines the vector transition polarisability, beta: the rank-1 part of the amplitude is written

\[ A^{1} = i\,\beta\,(\epsilon^L\times\epsilon^S)\cdot \matel{J_bM_b}{\bm\sigma}{J_aM_a}, \qquad \beta = \frac{A^1}{\sqrt{2}\,\redmatel{b}{\bm\sigma}{a}}, \]

the matrix element expressing the Wigner-Eckart factor of the rank-1 amplitude (any rank-1 operator would do; \( \bm\sigma \) is conventional).

Evaluated directly, as \( \bm\sigma = 2S \): the single-particle table of DiracOperator::s, contracted with the CI expansions by ReducedME. Nothing is assumed about L and S, which are not good quantum numbers for relativistic CI states.

Parameters

Psi_b,ib	Final CI state (solution ib of Psi_b).
Psi_a,ia	Initial CI state.
ci_basis	Single-particle basis of the CI expansion.

Returns

\( \redmatel{b}{\bm\sigma}{a} \).

Note

For a single valence electron the convention is to drop the radial overlap, so that \( \redmatel{7s}{\bm\sigma}{6s} = 2S_{\kappa\kappa} \) rather than zero (as in the dcp module, via Angular::S_kk). There is no consistent analogue for a multi-configuration state: forcing the radial overlaps to unity spoils even the diagonal matrix element, by adding cross-configuration terms.

So this vanishes for two states with no configuration in common (e.g., 3s3p and 3s4p). Then beta does not parameterise the rank-1 amplitude, and \( A^1 \) should be used directly.

A_K()

std::pair< double, double > CI::A_K	(	int	K,
		const PsiJPi &	Psi_b,
		std::size_t	ib,
		const 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 Integrals &	ints,
		const std::vector< Level > &	levels_to_remove = {},
		std::ostream &	outstream = std::cout
	)

Second-order amplitude \( A^K \) between two CI states, evaluated with CI mixed states.

Evaluates \( A^K \); see A_K_coefs. The sums over the intermediate spectrum are performed with the CI mixed states of solve_mixed_state, so they are complete: there is no sum over individual CI solutions, and no truncation of the spectrum. All intermediate states of a given (J, parity) share the same angular coefficient, so one mixed state per (J, parity) and per term is required.

Each sum is formed in two independent ways: with the mixed states of \( s \), and with those of \( t \). Writing \( \ket{A_s} \) for the state \( a \) plus its mixed state due to \( s \), these are the first-order parts of

\[ \redmatel{B_s}{t}{A_s} \qquad{\rm and}\qquad \redmatel{B_t}{s}{A_t}, \]

returned as the two elements of the pair. They must agree; the difference is a check on the numerics.

Covers, e.g., static, dynamic, and transition polarisabilities ( \( t = s = E1 \)) and PNC amplitudes ( \( s \) = PNC operator).

Parameters

K	Rank of the amplitude. It vanishes unless \( |k_t-k_s| \le K \le k_t+k_s \) and \( (J_b, K, J_a) \) satisfy the triangle rule.
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). For a frequency-dependent operator or RPA, the table should have been formed at omega.
s,s_me	The \( s \) operator, and its table (formed at omega_s).
omega	Frequency of \( t \). For a real transition carried entirely by \( t \) this is \( E_b - E_a \), for which the second denominator above is just \( E_b - E_n \).
omega_s	Frequency of \( s \). Energy conservation requires \( \omega + \omega_s = E_b - E_a \); it is zero for a transition carried entirely by \( t \), and \( -\omega \) for a dynamic polarisability.
ints	Integral tables, used to construct the CI Hamiltonian of each intermediate (J, parity); e.g., Wavefunction::CI_integrals().
levels_to_remove	CI levels to be removed from the intermediate states (see project_out), so that they may be treated separately - e.g., with experimental energies. See Level; the CI problem for those (J, parity) is solved here, as far as required.
outstream	Stream for progress and the intermediate sums.

Returns

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

Note

The parity selection rule \( \pi_a\pi_b = \pi_t\pi_s \) must hold, else the amplitude is zero.

If a state of an intermediate (J, parity) is degenerate with the denominator ( \( E_a + \omega_s \) or \( E_a + \omega \)), the mixed-states equation is singular, and its term in \( A^K \) is divergent. This cannot happen for operators of odd parity (as for polarisabilities and PNC), since then the intermediate states have the opposite parity to \( a \) and \( b \).

Corrections to the matrix elements (RPA, structure radiation, normalisation of states) enter through the single-particle tables.

A_K_core()

double CI::A_K_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.

The intermediate states of A_K carry no core hole. This is the missing term: a core electron \( c \) excited by one operator and de-excited by the other,

\[ \redmatel{c}{[t\times s]^0}{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], \]

\[ A^0_{\rm core} = \sqrt{[J]}\sum_c \sqrt{[j_c]}\, \redmatel{c}{[t\times s]^0}{c}, \]

with \( c_1, c_2 \) from A_K_coefs at \( (j_c, j_m, j_c) \).

The core is closed, \( J=0 \), so this is non-zero only for \( K=0 \) (which requires \( k_t=k_s \)) and only for \( b=a \), the valence factor being \( \braket{B}{A} \). For \( t=s=E1 \) it is the core polarisability. Apart from \( \sqrt{[J]} \) it is the same for every CI level, so it need only be evaluated once.

Parameters

K	Rank of the amplitude; zero unless \( K=0 \).
twoJ	2J of the CI state. The diagonal condition is left to the caller: zero unless the final and initial states are the same.
t,s	The two operators.
omega,omega_s	Frequency of each operator; see A_K.
core	Hole states \( c \); e.g., Wavefunction::core().
excited	Particle states \( m \): basis states above the Fermi level. Not restricted to the CI basis, and states occupied by the valence electrons are not removed - see A_K_cv.
dVt,dVs	RPA for each operator, solved at the frequency of that 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. Dressing both counts each RPA chain twice, since the sum over \( c \) already runs over every link of the chain. Same counting as the polarisability module.

No structure radiation: both lines here are core lines

A_K_cv()

double CI::A_K_cv	(	int	K,
		const PsiJPi &	Psi_b,
		std::size_t	ib,
		const PsiJPi &	Psi_a,
		std::size_t	ia,
		const DiracOperator::TensorOperator *	t,
		const DiracOperator::TensorOperator *	s,
		double	omega,
		double	omega_s,
		const std::vector< DiracSpinor > &	core,
		const std::vector< DiracSpinor > &	ci_basis,
		const ExternalField::CorePolarisation *	dVt = nullptr,
		const ExternalField::CorePolarisation *	dVs = nullptr
	)

Core-valence contribution to \( A^K \): the Pauli blocking of the core excitations by the valence electrons.

A_K_core sums over every particle state, including those occupied by the valence electrons. That excitation is blocked; the path that replaces it is one operator exciting a core electron to \( v' \), the other dropping a valence electron from \( v \) into the hole. This is a one-body operator in the valence space,

\[ \redmatel{v'}{[t\times s]^K_{cv}}{v} = \sum_c \left[ c_2\,\frac{\redmatel{v'}{s}{c}\redmatel{c}{t}{v}} {\en_{v'} - \omega_s - \en_c} + c_1\,\frac{\redmatel{v'}{t}{c}\redmatel{c}{s}{v}} {\en_{v'} - \omega - \en_c} \right], \qquad A^K_{cv} = \redmatel{B}{[t\times s]^K_{cv}}{A}, \]

with \( c_1, c_2 \) from A_K_coefs at \( (j_{v'}, j_c, j_v) \). The two terms take the coefficient of the opposite ordering, since the hole reverses the roles of the operators; its sign is cancelled by the reversed energy denominator. The contraction with the CI states is ReducedME.

For \( v'=v \) this is the blocking counter term, of weight \( n_v/[j_v] \) for occupation \( n_v \). The \( v' \ne v \) terms contribute for any \( K \), and between different CI states.

Parameters

K	Rank of the amplitude.
Psi_b,ib	Final CI state.
Psi_a,ia	Initial CI state.
t,s	The two operators.
omega,omega_s	Frequency of each operator; see A_K.
core	Hole states \( c \); e.g., Wavefunction::core().
ci_basis	Single-particle basis of the CI expansion.
dVt,dVs	RPA for each operator, solved at the frequency of that operator. May be nullptr.

Returns

\( A^K_{cv} \).

Note

RPA enters twice, on both matrix elements - unlike A_K_core. The blocked pair is a link of the RPA chain, and the chain continues on either side of it: dressing one vertex removes the chains that end on the blocked pair, but not those that pass through it.

No structure radiation; cf A_K_core.