ok

Matrix element of tensor operator: J_L(qr)*C^L.

Common matrix element for scattering problems. Defined as matrix element of J_L(qr)*C^L. Note: Does NOT include the Sqrt(2L+1) term (or 4pi). JL is spherical bessel function. Uses a look-up table for jL, faster than calculating on-the-fly. Note: different to other operators: rank (rank=L) and q are variable. Must use set_L_q (which sets L and q value) prior to use. Note: this is not thread-safe, so cannot parallelise over L or q.

#include <jL.hpp>

Inheritance diagram for DiracOperator::jL:

Public Member Functions
	jL (const Grid &r_grid, const Grid &q_grid, std::size_t max_l, bool subtract_one=false)
	Contruction takes radial grid, a q grid, and a maximum L. Fills lookup table.

	jL (const jL &other)
	Constructing from existing operator: copies JL table - faster. Can copy from a jL of different type (g0,g5 etc)

jL &	operator= (const jL &)=delete

std::size_t	L () const
	Current value of L (should = rank)

std::size_t	max_L () const
	Maximum L value in table.

const auto &	q_grid () const

const auto &	r_grid () const

virtual void	set_L_q (std::size_t L, double q)
	Sets the current L and q values for use. Note: NOT thread safe!

DiracSpinor	radial_rhs (const int kappa_a, const DiracSpinor &Fb) const override
	Computes the right-hand spinor dF_b for the radial integral.

virtual double	radialIntegral (const DiracSpinor &Fa, const DiracSpinor &Fb) const override final
	Radial integral R_ab, defined by RME = angularF(a,b) * radialIntegral(a,b).

double	rme (const DiracSpinor &a, const DiracSpinor &b, std::size_t L, double q) const
	Directly calculate reduced matrix element without needing to call set_L_q - this is thread safe.

bool	is_zero (const DiracSpinor &a, const DiracSpinor &b, std::size_t L) const
	Checks if specific ME is zero (when not useing set_L_q)

virtual double	angularF (const int ka, const int kb) const override
	Angular factor A_ab linking the radial integral to the RME.

virtual double	angularCff (int, int) const override
	Angular coefficient C_ff for the f_a*f_b term of the radial integral.

virtual double	angularCgg (int, int) const override
	Angular coefficient C_gg for the g_a*g_b term of the radial integral.

virtual double	angularCfg (int, int) const override
	Angular coefficient C_fg for the f_a*g_b term of the radial integral.

virtual double	angularCgf (int, int) const override
	Angular coefficient C_gf for the g_a*f_b term of the radial integral.

virtual std::string	name () const override
	Returns "name" of operator (e.g., 'E1')

std::string	units () const override final
	Returns units of operator as a string (usually au, may be MHz, etc.)

Public Member Functions inherited from DiracOperator::TensorOperator
	TensorOperator (const TensorOperator &)=default

TensorOperator &	operator= (const TensorOperator &)=default

	TensorOperator (TensorOperator &&)=default

TensorOperator &	operator= (TensorOperator &&)=default

bool	freqDependantQ () const
	Returns true if the operator is frequency-dependent (requires updateFrequency() calls).

bool	isZero (int ka, int kb) const
	Returns true if <a\|h\|b> = 0 by rank/parity selection rules.

bool	isZero (const DiracSpinor &Fa, const DiracSpinor &Fb) const
	Overload taking DiracSpinors; forwards to isZero(ka, kb).

bool	selectrion_rule (int twoJA, int piA, int twoJB, int piB) const
	Returns true if the matrix element is non-zero by angular momentum and parity selection rules (arguments are 2j and pi as integers).

virtual void	updateFrequency (const double)
	Updates the operator for a new frequency omega.

virtual void	updateRank (int)

const std::vector< double > &	getv () const
	Returns a const ref to the stored vector v.

double	getc () const
	Returns the "overall" constant c.

bool	imaginaryQ () const
	returns true if operator is imaginary (has imag MEs)

int	rank () const
	Rank k of operator.

int	parity () const
	returns parity, as integer (+1 or -1)

int	symm_sign (const DiracSpinor &Fa, const DiracSpinor &Fb) const
	returns relative sign between <a\|\|x\|\|b> and <b\|\|x\|\|a>

double	angularCxy (uint8_t x, uint8_t y, int kappa_a, int kappa_b) const
	Dispatches to angularCff/fg/gf/gg based on component indices x, y.

virtual std::unique_ptr< TensorOperator >	clone () const
	Returns a polymorphic copy of the operator at its current state, or nullptr if cloning is not supported by the derived class.

double	rme3js (int twoja, int twojb, int two_mb=1, int two_q=0) const
	3j-symbol factor linking the full ME to the RME.

double	rme3js (const DiracSpinor &Fa, const DiracSpinor &Fb, int two_mb=1, int two_q=0) const
	Overload of rme3js taking DiracSpinors.

DiracSpinor	reduced_rhs (const int ka, const DiracSpinor &Fb) const
	Returns angularF(ka,kb) * radial_rhs(ka,Fb); spinor-valued RME action on Fb, used in perturbation theory/TDHF.

DiracSpinor	reduced_lhs (const int ka, const DiracSpinor &Fb) const
	As reduced_rhs but for the conjugate direction; Fb * reduced_lhs(ka, Fb) = <b\|\|h\|\|a>.

double	reducedME (const DiracSpinor &Fa, const DiracSpinor &Fb) const
	Returns the reduced matrix element <a\|\|h\|\|b> = A_ab * R_ab.

double	fullME (const DiracSpinor &Fa, const DiracSpinor &Fb, std::optional< int > two_ma=std::nullopt, std::optional< int > two_mb=std::nullopt, std::optional< int > two_q=std::nullopt) const
	Returns "full" matrix element, for optional (ma, mb, q) [taken as int 2*]. If not specified, returns z-component (q=0), with ma=mb=min(ja,jb)

double	matel_factor (MatrixElementType type, int twoJa, int twoJb) const
	Returns the factor to convert a reduced ME to a different form (Reduced, Stretched, or HFConstant); see MatrixElementType.

double	matel_factor (MatrixElementType type, const DiracSpinor &Fa, const DiracSpinor &Fb) const
	Overload of matel_factor taking DiracSpinors.

Protected Attributes
std::size_t	m_max_l

LinAlg::Matrix< std::vector< double > >	m_j_lq_r

const Grid *	m_r_grid

const Grid *	m_q_grid

bool	m_subtract_one

Protected Attributes inherited from DiracOperator::TensorOperator
int	m_rank

Parity	m_parity

Realness	m_Realness

bool	m_freqDependantQ {false}

double	m_constant

std::vector< double >	m_vec

Additional Inherited Members
Protected Member Functions inherited from DiracOperator::TensorOperator
	TensorOperator (int rank_k, Parity pi, double constant=1.0, const std::vector< double > &vec={}, Realness RorI=Realness::real, bool freq_dep=false)
	Constructs a specific tensor operator. Called by derived classes.

Constructor & Destructor Documentation

◆ jL() [1/2]

DiracOperator::jL::jL	(	const Grid &	r_grid,
		const Grid &	q_grid,
		std::size_t	max_l,
		bool	subtract_one = `false`
	)

inline

Contruction takes radial grid, a q grid, and a maximum L. Fills lookup table.

◆ jL() [2/2]

DiracOperator::jL::jL ( const jL & other )

inline

Constructing from existing operator: copies JL table - faster. Can copy from a jL of different type (g0,g5 etc)

Member Function Documentation

◆ L()

std::size_t DiracOperator::jL::L ( ) const

inline

Current value of L (should = rank)

◆ max_L()

std::size_t DiracOperator::jL::max_L ( ) const

inline

Maximum L value in table.

◆ set_L_q()

virtual void DiracOperator::jL::set_L_q	(	std::size_t	L,
		double	q
	)

inlinevirtual

Sets the current L and q values for use. Note: NOT thread safe!

Reimplemented in DiracOperator::ig5jL, and DiracOperator::ig0g5jL.

◆ radial_rhs()

DiracSpinor DiracOperator::jL::radial_rhs	(	const int	kappa_a,
		const DiracSpinor &	Fb
	)		const

inlineoverridevirtual

Computes the right-hand spinor dF_b for the radial integral.

Returns \( \delta F_b \) such that the radial integral satisfies:

\[ R_{ab} = F_a \cdot \delta F_b = \int_0^\infty \left(f_a\,\delta f_b + g_a\,\delta g_b\right)\,{\rm d}r \]

The default implementation constructs \( \delta F_b \) using the stored radial function \( v(r) \) and the angular coefficients:

\[ \delta F_b(r) = c\,v(r) \begin{pmatrix} C_{ff}\,f_b(r) + C_{fg}\,g_b(r) \\ C_{gf}\,f_b(r) + C_{gg}\,g_b(r) \end{pmatrix} \]

This is used by reduced_rhs() to build \( \langle a \| \hat{h} \| b \rangle \) as a spinor-valued quantity, enabling perturbation theory and TDHF. Override this for operators whose radial structure cannot be expressed in this standard form.

Parameters

kappa_a	Relativistic quantum number \( \kappa_a \) of the bra state (needed to evaluate the angular coefficients).
Fb	Ket DiracSpinor \( F_b \) .

Returns: DiracSpinor \( \delta F_b \) .

Warning: If this is overridden, radialIntegral() should also be overridden consistently (and vice versa), so that reducedME() and reduced_rhs() remain consistent.

Reimplemented from DiracOperator::TensorOperator.

◆ radialIntegral()

virtual double DiracOperator::jL::radialIntegral	(	const DiracSpinor &	Fa,
		const DiracSpinor &	Fb
	)		const

inlinefinaloverridevirtual

Radial integral R_ab, defined by RME = angularF(a,b) * radialIntegral(a,b).

Returns the radial part \( R_{ab} \) of the reduced matrix element:

\[ \langle a \| \hat{h} \| b \rangle = A_{ab} \cdot R_{ab} \]

where \( A_{ab} \) is angularF().

The default implementation evaluates \( R_{ab} = F_a \cdot \delta F_b \) , using the default radial structure:

\[ R_{ab} = c\int_0^\infty v(r)\left( C_{ff}\,f_a f_b + C_{fg}\,f_a g_b + C_{gf}\,g_a f_b + C_{gg}\,g_a g_b \right)\,{\rm d}r \]

Override this for operators that do not fit this standard form. If radial_rhs() is also overridden, both must remain mutually consistent.

Warning: If radial_rhs() is overridden but radialIntegral() is not (or vice versa), reducedME() and reduced_rhs() will give inconsistent results.

Reimplemented from DiracOperator::TensorOperator.

◆ rme()

double DiracOperator::jL::rme	(	const DiracSpinor &	a,
		const DiracSpinor &	b,
		std::size_t	L,
		double	q
	)		const

inline

Directly calculate reduced matrix element without needing to call set_L_q - this is thread safe.

◆ is_zero()

bool DiracOperator::jL::is_zero	(	const DiracSpinor &	a,
		const DiracSpinor &	b,
		std::size_t	L
	)		const

inline

Checks if specific ME is zero (when not useing set_L_q)

◆ angularF()

virtual double DiracOperator::jL::angularF	(	const int	,
		const int
	)		const

inlineoverridevirtual

Angular factor A_ab linking the radial integral to the RME.

All derived operators must implement this. It gives the purely angular part of the reduced matrix element:

\[ \langle a \| \hat{h} \| b \rangle \equiv A_{ab} \cdot R_{ab} \]

where \( R_{ab} \) is returned by radialIntegral(). For most operators, \( A_{ab} \) is a product of Clebsch-Gordan / 3j coefficients and depends only on \( \kappa_a, \kappa_b \) (and the rank \( k \) and parity \( \pi \) of the operator).

Note: This is a pure virtual function – every derived operator must provide an implementation.

Implements DiracOperator::TensorOperator.

Reimplemented in DiracOperator::ig0g5jL, and DiracOperator::ig5jL.

◆ angularCff()

virtual double DiracOperator::jL::angularCff	(	int	kappa_a,
		int	kappa_b
	)		const

inlineoverridevirtual

Angular coefficient C_ff for the f_a*f_b term of the radial integral.

The default radial integral is structured as:

\[ R_{ab} = c\int_0^\infty v(r)\left( C_{ff}\,f_a f_b + C_{fg}\,f_a g_b + C_{gf}\,g_a f_b + C_{gg}\,g_a g_b \right)\,{\rm d}r \]

These coefficients are often constants, but may depend on \( \kappa_a, \kappa_b \) for operators with angular-momentum-dependent coupling between large and small components (e.g., spin-dependent operators). Override in derived classes as needed.

Parameters