ok

Stores Wavefunction (set of valence orbitals, grid, HF etc.)

Also holds the Hartree-Fock potential (and core), basis, spectrum, correlation potential (Sigma), and CI solutions.

Construction:

Set of GridParameters [see Maths/Grid]
Nuclear::Nucleus [see Potentials/NuclearPotentials]
var_alpha = \(\lambda\), \(\alpha = \lambda\alpha_0\)
run_label: Optional label for output identity - for distinguishing outputs with different parameters

Copying:: A copy has its own HF object (and core), but its Sigma (if formed) still refers to the HF object of the original (see copySigma). Applying Sigma does not use the HF object, but forming or writing Sigma does: these then use the original core (e.g., without any changes to the copy's HF), and the original must outlive the copy.

#include <Wavefunction.hpp>

Public Member Functions
	Wavefunction (std::shared_ptr< const Grid > grid, const Nuclear::Nucleus &nucleus, double var_alpha=1.0, const std::string &run_label="")
	Construct with a Grid [shared resource] and a nucleus.

	Wavefunction (const GridParameters &gridparams, const Nuclear::Nucleus &nucleus, double var_alpha=1.0, const std::string &run_label="")
	As above, but Grid is constructed here using given parameters.

	Wavefunction ()
	Default grid and nucleus (see GridParameters, Nuclear::Nucleus)

const Grid &	grid () const
	Returns a const reference to the radial grid.

std::shared_ptr< const Grid >	grid_sptr () const
	Copy of shared_ptr to grid [shared resource].

double	alpha () const
	Local value of fine-structure constant.

double	dalpha2 () const
	Variation in alpha^2 : x = (alpha/alpha_0)^2 - 1.

const Nuclear::Nucleus &	nucleus () const
	Returns Nuclear::nucleus object (contains nuc. parameters)

int	Znuc () const
	Nuclear charge, Z.

int	Anuc () const
	Nuclear mass number, A.

double	get_rrms () const
	Nuclear rms charge radius, in fm (femtometres)

const std::vector< DiracSpinor > &	core () const
	Core orbitals (frozen HF core); empty if HF not set.

const std::vector< DiracSpinor > &	valence () const
	Valence orbitals (HF or Brueckner orbitals)

std::vector< DiracSpinor > &	valence ()

const std::vector< DiracSpinor > &	hf_valence () const
	HF valence orbitals: copy of valence, as last solved by solve_valence.

const std::vector< DiracSpinor > &	basis () const
	Basis, eigenstates of HF potential. Used for MBPT.

std::vector< DiracSpinor > &	basis ()

const std::vector< DiracSpinor > &	spectrum () const
	Spectrum: like basis, but includes correlations.

std::vector< DiracSpinor > &	spectrum ()

const std::vector< CI::PsiJPi > &	CIwfs () const
	CI solutions, one for each {J, parity}; empty unless CI was run.

const CI::PsiJPi *	CIwf (int J, int parity) const
	CI solutions for given J and parity; nullptr if not found.

const CI::Integrals &	CI_integrals () const
	Integral tables used to construct the CI Hamiltonians.

const std::vector< double > &	vnuc () const
	Nuclear potential.

const HF::HartreeFock *	vHF () const
	Returns ptr to Hartree Fock (class); nullptr if not set.

HF::HartreeFock *	vHF ()

std::vector< double >	vlocal (int l=0) const
	Local part of potential, e.g., Vl = Vnuc + Vdir + Vrad_el(l).

std::vector< double >	Hmag (int l=0) const
	QED Magnetic form factor.

const QED::RadPot *	vrad () const
	Pointer to QED radiative potential. May be nullptr.

QED::RadPot *	vrad ()

const MBPT::CorrelationPotential *	Sigma () const
	Returns ptr to (const) Correlation Potential, Sigma; nullptr if none.

MBPT::CorrelationPotential *	Sigma ()

int	Ncore () const
	Number of electrons in the core.

const DiracSpinor *	getState (int n, int k) const
	Finds requested state (n, kappa); returns nullptr if not found.

const DiracSpinor *	getState (std::string_view state) const
	As above, but takes 'short symbol' (e.g., 6s+, 6p-)

double	FermiLevel () const
	Returns energy location of the "Fermi Level".

double	energy_gap () const
	Energy gap between lowest valence + highest core state.

std::string	coreConfiguration () const
	Returns full core configuration.

std::string	coreConfiguration_nice () const
	Returns core configuration, in nice output notation.

std::string	atom () const
	String of atom info (e.g., "Cs, Z=55 A=133")

std::string	atomicSymbol () const
	e.g., "Cs"

std::string	identity () const
	Atomic symbol, including core ionisation degree and run_label.

const std::string &	run_label () const
	Optional label for output identity (see identity)

int	ion_degree (int num_val) const
	0 for neutral, 1 for singly-ionised etc. (num_val: valence electrons)

std::string	ion_symbol (int num_val) const
	I for neutral, II for singly-ionised etc. (num_val: valence electrons)

int	Zion () const
	Effective charge (for core) = Z-N_core.

void	printCore () const
	Prints table of core orbitals + energies etc.

void	printValence (const std::vector< DiracSpinor > &tmp_orbitals={}) const
	Prints table of valence orbitals + energies etc.

void	printBasis (const std::vector< DiracSpinor > &the_basis) const
	Prints table of Basis/Spectrum orbitals, compares to HF orbitals.

bool	isInCore (int n, int k) const
	Check if a state (n, kappa) is in the core list.

bool	isInValence (int n, int k) const
	Check if a state (n, kappa) is in the valence list.

std::vector< double >	coreDensity () const
	Calculates rho(r) = sum_c psi^2(r) for core states, c={n,k,m}.

double	coreEnergyHF () const
	Calculates HF core energy.

void	set_HF (HF::Method method=HF::Method::HartreeFock, const std::string &in_core="", std::optional< HF::Breit::Params > breit_params=std::nullopt, double eps_HF=1.0e-14, bool print=true)
	Initialises HF object and populates core orbitals.

void	solve_core (bool print=true)
	Performs Hartree-Fock procedure for core (nothing if HF not set)

void	solve_core (HF::Method method, const std::string &in_core="", std::optional< HF::Breit::Params > breit_params=std::nullopt, double eps_HF=1.0e-14, bool print=true)
	This version will first set_HF(), then solve_core()

void	solve_core (const std::string &method, const std::string &in_core="", std::optional< HF::Breit::Params > breit_params=std::nullopt, double eps_HF=1.0e-14, bool print=true)
	This is an overloaded member function, provided for convenience. It differs from the above function only in what argument(s) it accepts.

void	solve_valence (const std::string &in_valence_str="", const bool print=true)
	Performs Hartree-Fock procedure for valence; populates valence.

void	solve_exotic (const std::string &in_exotic_str, double mass=PhysConst::m_muon, bool screening=true, bool print=true)
	Solves for exotic atoms (e.g., muonic), including screening.

void	hartreeFockBrueckner (const bool print=true)
	Forms Brueckner valence orbitals; replaces existing valence states.

void	fitSigma_hfBrueckner (const std::string &valence_list, const std::vector< double > &fit_energies, bool separate_lambda=false)
	First, fits Sigma to energies, then forms fitted Brueckner orbitals.

void	radiativePotential (QED::RadPot::Scale s, double rcut, double scale_rN, const std::vector< double > &x_spd, bool do_readwrite=true, bool print=true)
	OLD: deprecated.

void	radiativePotential (const IO::InputBlock &qed_input, bool do_readwrite, bool print)
	Calculates radiative potential, adds to HF potential.

void	formBasis (const SplineBasis::Parameters &params)
	Calculates + populates basis [see BSplineBasis].

void	formSpectrum (const SplineBasis::Parameters &params)
	Calculates + populates Spectrum [see BSplineBasis].

void	formSigma (int nmin_core=1, double r0=1.0e-4, double rmax=30.0, int stride=4, bool each_valence=false, bool include_G=false, bool include_Breit=false, int n_max_breit=0, const std::vector< double > &lambdas={}, const std::vector< double > &fk={}, const std::vector< double > &etak={}, bool read_write=true, const std::string &fname="", bool FeynmanQ=false, bool ScreeningQ=false, bool hole_particleQ=false, int lmax=6, double omre=-0.2, double w0=0.01, double wratio=1.5, bool complex_green=false, const std::optional< IO::InputBlock > &ek=std::nullopt, const std::string &ladder_file="", bool derivative=false, bool print_de=false, bool fk_both_lines=false, bool tail_extrapolation=false)
	Forms + stores correlation potential Sigma.

void	copySigma (const MBPT::CorrelationPotential *const Sigma)
	Copies correlation potential Sigma (nothing if `Sigma` is null).

void	update_Vnuc (const std::vector< double > &v_new)
	Replaces nuclear potential, in both Wavefunction and HartreeFock.

std::tuple< double, double >	lminmax_core_range (int l, double eps=0.0) const
	Returns [min,max] r values for which the core density is > cutoff.

double	H0ab (const DiracSpinor &Fa, const DiracSpinor &Fb) const
	Returns <a\|H\|b> for Hamiltonian H (NOT Sigma, Breit, or exchange).

double	H0ab (const DiracSpinor &Fa, const DiracSpinor &dFa, const DiracSpinor &Fb, const DiracSpinor &dFb) const
	As above, but with the radial derivatives given.

double	Hab (const DiracSpinor &Fa, const DiracSpinor &Fb) const
	Returns <a\|H\|b>, including exchange, Sigma, and Breit.

void	ConfigurationInteraction (const IO::InputBlock &input)
	Runs the CI+MBPT routines; stores wavefunctions.

nlohmann::json	output_to_json (const std::string &out_name="ampsci_output.json")
	Writes wavefunction information to json file.

Constructor & Destructor Documentation

◆ Wavefunction() [1/3]

Wavefunction::Wavefunction	(	std::shared_ptr< const Grid >	grid,
		const Nuclear::Nucleus &	nucleus,
		double	var_alpha = `1.0`,
		const std::string &	run_label = `""`
	)

Construct with a Grid [shared resource] and a nucleus.

Nucleus: isotope data etc. Optional fractional variation in alpha: alpha = var_alpha * alpha_0, alpha_0 =~ 1/137. Forms the nuclear potential; orbitals are formed later (see set_HF, solve_core).

Parameters

grid	Radial grid (shared resource)
nucleus	Nuclear parameters (isotope, charge distribution)
var_alpha	Ratio alpha/alpha_0
run_label	Optional label for output identity (see identity)

◆ Wavefunction() [2/3]

Wavefunction::Wavefunction	(	const GridParameters &	gridparams,
		const Nuclear::Nucleus &	nucleus,
		double	var_alpha = `1.0`,
		const std::string &	run_label = `""`
	)

As above, but Grid is constructed here using given parameters.

◆ Wavefunction() [3/3]

Wavefunction::Wavefunction ( )

inline

Default grid and nucleus (see GridParameters, Nuclear::Nucleus)

Member Function Documentation

◆ grid()

const Grid & Wavefunction::grid ( ) const

inline

Returns a const reference to the radial grid.

◆ grid_sptr()

std::shared_ptr< const Grid > Wavefunction::grid_sptr ( ) const

inline

Copy of shared_ptr to grid [shared resource].

Used when we want to construct a new object that shares this grid.

◆ alpha()

double Wavefunction::alpha ( ) const

inline

Local value of fine-structure constant.

◆ dalpha2()

double Wavefunction::dalpha2 ( ) const

inline

Variation in alpha^2 : x = (alpha/alpha_0)^2 - 1.

◆ nucleus()

const Nuclear::Nucleus & Wavefunction::nucleus ( ) const

inline

Returns Nuclear::nucleus object (contains nuc. parameters)

◆ Znuc()

int Wavefunction::Znuc ( ) const

inline

Nuclear charge, Z.

◆ Anuc()

int Wavefunction::Anuc ( ) const

inline

Nuclear mass number, A.

◆ get_rrms()

double Wavefunction::get_rrms ( ) const

inline

Nuclear rms charge radius, in fm (femtometres)

◆ core()

const std::vector< DiracSpinor > & Wavefunction::core ( ) const

inline

Core orbitals (frozen HF core); empty if HF not set.

◆ valence()

const std::vector< DiracSpinor > & Wavefunction::valence ( ) const

inline

Valence orbitals (HF or Brueckner orbitals)

◆ hf_valence()

const std::vector< DiracSpinor > & Wavefunction::hf_valence ( ) const

inline

HF valence orbitals: copy of valence, as last solved by solve_valence.

◆ basis()

const std::vector< DiracSpinor > & Wavefunction::basis ( ) const

inline

Basis, eigenstates of HF potential. Used for MBPT.

Includes Breit and QED (if they are included), but not correlations.

◆ spectrum()

const std::vector< DiracSpinor > & Wavefunction::spectrum ( ) const

inline

Spectrum: like basis, but includes correlations.

◆ CIwfs()

const std::vector< CI::PsiJPi > & Wavefunction::CIwfs ( ) const

inline

CI solutions, one for each {J, parity}; empty unless CI was run.

◆ CIwf()

const CI::PsiJPi * Wavefunction::CIwf	(	int	J,
		int	parity
	)		const

inline

CI solutions for given J and parity; nullptr if not found.

J is integer (CI is for two-electron systems); parity is +1 or -1.

◆ CI_integrals()

const CI::Integrals & Wavefunction::CI_integrals ( ) const

inline

Integral tables used to construct the CI Hamiltonians.

Empty unless CI was run; see CI::Integrals.

◆ vnuc()

const std::vector< double > & Wavefunction::vnuc ( ) const

inline

Nuclear potential.

Only provide const version, since HF and WF version of vnuc must be kept in sync (see update_Vnuc).

◆ vHF()

const HF::HartreeFock * Wavefunction::vHF ( ) const

inline

Returns ptr to Hartree Fock (class); nullptr if not set.

◆ vlocal()

std::vector< double > Wavefunction::vlocal ( int l = 0 ) const

Local part of potential, e.g., Vl = Vnuc + Vdir + Vrad_el(l).

Can be l-dependent. Returns a copy. Just Vnuc if HF is not set.

◆ Hmag()

std::vector< double > Wavefunction::Hmag ( int l = 0 ) const

QED Magnetic form factor.

May return empty vector (no QED, or HF not set). Not typically l-dependent, but may be in future. Returns a copy.

◆ vrad()

const QED::RadPot * Wavefunction::vrad ( ) const

inline

Pointer to QED radiative potential. May be nullptr.

◆ Sigma()

const MBPT::CorrelationPotential * Wavefunction::Sigma ( ) const

inline

Returns ptr to (const) Correlation Potential, Sigma; nullptr if none.

◆ Ncore()

int Wavefunction::Ncore ( ) const

Number of electrons in the core.

◆ getState() [1/2]

const DiracSpinor * Wavefunction::getState	(	int	n,
		int	k
	)		const

Finds requested state (n, kappa); returns nullptr if not found.

Searches valence first, then core.

◆ getState() [2/2]

const DiracSpinor * Wavefunction::getState ( std::string_view state ) const

As above, but takes 'short symbol' (e.g., 6s+, 6p-)

◆ FermiLevel()

double Wavefunction::FermiLevel ( ) const

Returns energy location of the "Fermi Level".

Energy half way between core/valence. Defined: 0.5*( max(e_core) + min(e_valence)), plus a small offset (1e-5; or 1e-3 if there are no valence states, in which case min(e_valence) is taken as zero). If there is no core: 1.1*min(e_valence). Should be -ve.

◆ energy_gap()

double Wavefunction::energy_gap ( ) const

Energy gap between lowest valence + highest core state.

e(v) - e(c) [should be positive]. Zero if there are no core or no valence states.

◆ coreConfiguration()

std::string Wavefunction::coreConfiguration ( ) const

inline

Returns full core configuration.

◆ coreConfiguration_nice()

std::string Wavefunction::coreConfiguration_nice ( ) const

inline

Returns core configuration, in nice output notation.

◆ atom()

std::string Wavefunction::atom ( ) const

inline

String of atom info (e.g., "Cs, Z=55 A=133")

◆ atomicSymbol()

std::string Wavefunction::atomicSymbol ( ) const

inline

e.g., "Cs"

◆ identity()

std::string Wavefunction::identity ( ) const

Atomic symbol, including core ionisation degree and run_label.

Used as the base name for output and cache files (Sigma, Qk, CI, etc.). Form: symbol + Zion + Breit/QED flag ("b", "q", or "bq") + alpha variation ("_x", with x = (alpha/alpha_0)^2 to 6 significant figures, if |dalpha2()| > 1e-5) + "_run_label" (if not empty). E.g., "Cs1", "Cs1bq_label", "Cs1_1.01".

◆ run_label()

const std::string & Wavefunction::run_label ( ) const

inline

Optional label for output identity (see identity)

◆ ion_degree()

int Wavefunction::ion_degree ( int num_val ) const

inline

0 for neutral, 1 for singly-ionised etc. (num_val: valence electrons)

◆ ion_symbol()

std::string Wavefunction::ion_symbol ( int num_val ) const

inline

I for neutral, II for singly-ionised etc. (num_val: valence electrons)

◆ Zion()

int Wavefunction::Zion ( ) const

inline

Effective charge (for core) = Z-N_core.

◆ printCore()

void Wavefunction::printCore ( ) const

Prints table of core orbitals + energies etc.

◆ printValence()

void Wavefunction::printValence ( const std::vector< DiracSpinor > & tmp_orbitals = {} ) const

Prints table of valence orbitals + energies etc.

Can optionally give it any list of orbitals to print (default: valence).

◆ printBasis()

void Wavefunction::printBasis ( const std::vector< DiracSpinor > & the_basis ) const

Prints table of Basis/Spectrum orbitals, compares to HF orbitals.

◆ isInCore()

bool Wavefunction::isInCore	(	int	n,
		int	k
	)		const

Check if a state (n, kappa) is in the core list.

◆ isInValence()

bool Wavefunction::isInValence	(	int	n,
		int	k
	)		const

Check if a state (n, kappa) is in the valence list.

◆ coreDensity()

std::vector< double > Wavefunction::coreDensity ( ) const

Calculates rho(r) = sum_c psi^2(r) for core states, c={n,k,m}.

◆ coreEnergyHF()

double Wavefunction::coreEnergyHF ( ) const

Calculates HF core energy.

Sum of core orbital energies (weighted by occupation), minus the double-counted Coulomb (direct + exchange) and Breit (if included) interactions. QED (if included in the core) enters through the orbital energies. Zero if HF is not set.

◆ set_HF()

void Wavefunction::set_HF	(	HF::Method	method = `HF::Method::HartreeFock`,
		const std::string &	in_core = `""`,
		std::optional< HF::Breit::Params >	breit_params = `std::nullopt`,
		double	eps_HF = `1.0e-14`,
		bool	print = `true`
	)

Initialises HF object and populates core orbitals.

Does not solve HF equations (see solve_core). Replaces any existing HF object, including any radiative potential (QED is added later; see radiativePotential).

Parameters

method	Mean-field method (HartreeFock, Hartree, KohnSham..)
in_core	Core configuration, "[X],nLm,nLm,...", e.g., "[Xe]" or "[Xe],6s1". A ':' may replace one of the commas: shells after it are in the core, but may also be valence states (see solve_valence)
breit_params	Breit parameters (nullopt: no Breit). HartreeFock method only: ignored (with a warning) otherwise
eps_HF	HF convergence goal (if >= 1, taken as the exponent: 10^-eps_HF)
print	Print method details

Warning: Aborts if the core configuration is invalid, or has more than Z electrons.

◆ solve_core() [1/3]

void Wavefunction::solve_core ( bool print = true )

Performs Hartree-Fock procedure for core (nothing if HF not set)

◆ solve_core() [2/3]

void Wavefunction::solve_core	(	HF::Method	method,
		const std::string &	in_core = `""`,
		std::optional< HF::Breit::Params >	breit_params = `std::nullopt`,
		double	eps_HF = `1.0e-14`,
		bool	print = `true`
	)

This version will first set_HF(), then solve_core()

◆ solve_core() [3/3]

void Wavefunction::solve_core	(	const std::string &	method,
		const std::string &	in_core = `""`,
		std::optional< HF::Breit::Params >	breit_params = `std::nullopt`,
		double	eps_HF = `1.0e-14`,
		bool	print = `true`
	)

This is an overloaded member function, provided for convenience. It differs from the above function only in what argument(s) it accepts.

◆ solve_valence()

void Wavefunction::solve_valence	(	const std::string &	in_valence_str = `""`,
		const bool	print = `true`
	)

Performs Hartree-Fock procedure for valence; populates valence.

Adds the states in in_valence_str (basis-string format, e.g., "7sp5d": up to n=7 for s and p, n=5 for d) that are not already in valence, and not in the core (except those after the ':' in the core string; see set_HF). Then solves for all valence states (existing ones are re-solved) in the HF potential, without Sigma, and copies them to hf_valence. If HF is not set, solves in the nuclear potential (H-like).

◆ solve_exotic()

void Wavefunction::solve_exotic	(	const std::string &	in_exotic_str,
		double	mass = `PhysConst::m_muon`,
		bool	screening = `true`,
		bool	print = `true`
	)

Solves for exotic atoms (e.g., muonic), including screening.

Resulting states are included in valence (flagged as exotic); the screening also updates core.

Each state is first solved in the nuclear (+QED) potential. With screening (and a non-empty core), the core HF equations are iterated including the direct potential of the first exotic state; the HF object (core) is then replaced by the screened one (Vnuc itself unchanged), and the exotic states are re-solved in the screened potential.

Note: The exotic states are just added to the valence list, so they can be used more simply with all the modules. However, be careful; for example, RPA will now be meaningless for these states!

Parameters

in_exotic_str	States, in basis-string format (e.g., "1s2p")
mass	Mass of the exotic particle, in units of m_e
screening	Include screening (core and exotic states)
print	Print energies and screening iterations

Note: With screening, any (electron) valence states already present are not re-solved, so are no longer consistent with the screened core.

◆ hartreeFockBrueckner()

void Wavefunction::hartreeFockBrueckner ( const bool print = true )

Forms Brueckner valence orbitals; replaces existing valence states.

\( (H_{\rm HF} + \Sigma)\ket{nk} = \varepsilon\ket{nk} \). Does nothing if Sigma has not been formed; prints a warning and does nothing if HF is not set.

Note: Sigma is only formed for the HartreeFock method (see formSigma).

◆ fitSigma_hfBrueckner()

void Wavefunction::fitSigma_hfBrueckner	(	const std::string &	valence_list,
		const std::vector< double > &	fit_energies,
		bool	separate_lambda = `false`
	)

First, fits Sigma to energies, then forms fitted Brueckner orbitals.

For each of the first fit_energies.size() valence states (in order), finds the scaling factor lambda such that the Brueckner energy, with lambda*Sigma, matches the given energy (damped iteration; lambda kept within [0.25, 2.5]). By default, fitted lambdas are averaged over each fine-structure pair (same n and l); separate_lambda = true keeps them separate. Then calls hartreeFockBrueckner.

Parameters

valence_list	Unused
fit_energies	Energies to fit to (au), in the same order as valence; sign ignored (binding energies)
separate_lambda	Keep lambdas of fine-structure pairs separate

Note: Each fitted state needs its own Sigma (same kappa and n; see each_valence in formSigma). Otherwise it is not fitted (lambda = 1 is printed): its energy is that with the lower-n Sigma.; Prints a warning and does nothing if Sigma has not been formed.

◆ radiativePotential() [1/2]

void Wavefunction::radiativePotential	(	QED::RadPot::Scale	s,
		double	rcut,
		double	scale_rN,
		const std::vector< double > &	x_spd,
		bool	do_readwrite = `true`,
		bool	print = `true`
	)

OLD: deprecated.

◆ radiativePotential() [2/2]

void Wavefunction::radiativePotential	(	const IO::InputBlock &	qed_input,
		bool	do_readwrite,
		bool	print
	)

Calculates radiative potential, adds to HF potential.

Replaces any existing radiative potential in HF. Included in subsequent HF solutions: call before solve_core to include QED in the core, after it for valence only. Prints a warning and does nothing if HF is not set.

Parameters

qed_input	QED options (see QED::ConstructRadPot)
do_readwrite	Read/write the potential from/to file
print	Print details

◆ formBasis()

void Wavefunction::formBasis ( const SplineBasis::Parameters & params )

Calculates + populates basis [see BSplineBasis].

◆ formSpectrum()

void Wavefunction::formSpectrum ( const SplineBasis::Parameters & params )

Calculates + populates Spectrum [see BSplineBasis].

◆ formSigma()

void Wavefunction::formSigma	(	int	nmin_core = `1`,
		double	r0 = `1.0e-4`,
		double	rmax = `30.0`,
		int	stride = `4`,
		bool	each_valence = `false`,
		bool	include_G = `false`,
		bool	include_Breit = `false`,
		int	n_max_breit = `0`,
		const std::vector< double > &	lambdas = `{}`,
		const std::vector< double > &	fk = `{}`,
		const std::vector< double > &	etak = `{}`,
		bool	read_write = `true`,
		const std::string &	fname = `""`,
		bool	FeynmanQ = `false`,
		bool	ScreeningQ = `false`,
		bool	hole_particleQ = `false`,
		int	lmax = `6`,
		double	omre = `-0.2`,
		double	w0 = `0.01`,
		double	wratio = `1.5`,
		bool	complex_green = `false`,
		const std::optional< IO::InputBlock > &	ek = `std::nullopt`,
		const std::string &	ladder_file = `""`,
		bool	derivative = `false`,
		bool	print_de = `false`,
		bool	fk_both_lines = `false`,
		bool	tail_extrapolation = `false`
	)

Forms + stores correlation potential Sigma.

Constructs Sigma (see MBPT::CorrelationPotential), reading it from file if it exists, and forms it for the valence states (see each_valence, ek). Then applies lambdas, optionally prints the energy shifts, and writes Sigma to file. Requires the basis (formBasis) for the Goldstone sums. Does nothing if the core is empty or HF is not set. Requires the HartreeFock method: prints a warning and does nothing otherwise.

Parameters

nmin_core	Minimum n of core states included in polarisation
r0	Sub-grid first point (aB)
rmax	Sub-grid last point (aB)
stride	Sub-grid stride: every stride-th point of the full grid
each_valence	Form Sigma for each valence state; else for the first valence state of each kappa. Ignored if `ek` is given
include_G	Include small-component (g) parts of Sigma
include_Breit	Include two-body Breit correction, B2 (only if HF includes Breit)
n_max_breit	Maximum n of excited states in B2 (0: all)
lambdas	Scaling factors, lambda, for each Sigma, in the order formed
fk	Effective screening factors for exchange, for each k (empty: calculated for Feynman with screening)
etak	Effective hole-particle factors, for each k (Goldstone only)
read_write	Read Sigma from file (if it exists) and write it; false: no file
fname	Sigma file name (empty: identity() plus an extension encoding the method)
FeynmanQ	Feynman method for direct diagrams; else Goldstone
ScreeningQ	Include all-order screening (Feynman only)
hole_particleQ	Include hole-particle interaction (Feynman only)
lmax	Maximum l of internal lines (Feynman)
omre	Real part of frequency for the Im(w) integration (Feynman)
w0	First non-zero point of the Im(w) grid (Feynman)
wratio	Ratio of the logarithmic Im(w) grid (Feynman)
complex_green	Solve Green's function directly at complex energy (see MBPT::FeynmanOptions)
ek	State/energy pairs (energies in au): form Sigma for each at the given energy, instead of the valence energies
ladder_file	File with ladder correction, Sigma_L (empty: none)
derivative	Also form dSigma/dE
print_de	Print MBPT(2) energy corrections for valence
fk_both_lines	Apply fk to both Coulomb lines of the exchange diagrams (else outer only)
tail_extrapolation	Extrapolate Sigma beyond the sub-grid as a_eff/r^4 (experimental)

Note: If read from file, the sub-grid and include_G are taken from the file (see MBPT::CorrelationPotential).

◆ copySigma()

void Wavefunction::copySigma ( const MBPT::CorrelationPotential *const Sigma )

inline

Copies correlation potential Sigma (nothing if Sigma is null).

E.g., to re-use Sigma from another Wavefunction (different isotope).

Note: The copy still refers to the Hartree-Fock object it was formed with, which must outlive it if it is used to form or write Sigma.

◆ update_Vnuc()

void Wavefunction::update_Vnuc ( const std::vector< double > & v_new )

inline

Replaces nuclear potential, in both Wavefunction and HartreeFock.

Allows extra potential to be added to Vnuc: pass vnuc() + extra. Updates both in Wavefunction and HartreeFock, which must be kept in sync. Orbitals are not re-solved, and nucleus() is not updated.

nb: two versions of Vnuc...

◆ lminmax_core_range()

std::tuple< double, double > Wavefunction::lminmax_core_range	(	int	l,
		double	eps = `0.0`
	)		const

Returns [min,max] r values for which the core density is > cutoff.

Core density for given l; cutoff = eps*max_value. Returns the r values (au) for which the value of \( \rho(r) = \sum|\psi|^2 \) drops below cutoff. Sum goes over all m for given l. Cut-off defined as eps*max, where max is maximum value for rho(r). Set l<0 to get for all l (entire core). Returns the full grid range if there are no core states with given l (or if eps = 0).

◆ H0ab() [1/2]

double Wavefunction::H0ab	(	const DiracSpinor &	Fa,
		const DiracSpinor &	Fb
	)		const

Returns <a|H|b> for Hamiltonian H (NOT Sigma, Breit, or exchange).

Includes the local potential, vlocal (incl. Rad.pot), and the QED magnetic form factor, Hmag. Uses the symmetric form of the derivative term (integrated by parts). Zero if Fa and Fb have different kappa.

◆ H0ab() [2/2]

double Wavefunction::H0ab	(	const DiracSpinor &	Fa,
		const DiracSpinor &	dFa,
		const DiracSpinor &	Fb,
		const DiracSpinor &	dFb
	)		const

As above, but with the radial derivatives given.

dFa and dFb are the derivatives (d/dr) of Fa and Fb; only their g components are used.

◆ Hab()

double Wavefunction::Hab	(	const DiracSpinor &	Fa,
		const DiracSpinor &	Fb
	)		const

Returns <a|H|b>, including exchange, Sigma, and Breit.

H0ab, plus HF exchange (if HF is non-local), Sigma (if formed, including lambda), and Breit (if included in HF). Zero if kappa differs.

◆ ConfigurationInteraction()

void Wavefunction::ConfigurationInteraction ( const IO::InputBlock & input )

Runs the CI+MBPT routines; stores wavefunctions.

◆ output_to_json()

nlohmann::json Wavefunction::output_to_json ( const std::string & out_name = "ampsci_output.json" )

Writes wavefunction information to json file.

Grid, nuclear parameters, and core, valence, and exotic (muon) orbitals, with a description of each field ("metadata").

Parameters

out_name If given (not empty), will print to that file

Returns: The json object.

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

Wavefunction/Wavefunction.hpp
Wavefunction/Wavefunction.cpp

ok

Public Member Functions

Constructor & Destructor Documentation

◆ Wavefunction() [1/3]

◆ Wavefunction() [2/3]

◆ Wavefunction() [3/3]

Member Function Documentation

◆ grid()

◆ grid_sptr()

◆ alpha()

◆ dalpha2()

◆ nucleus()

◆ Znuc()

◆ Anuc()

◆ get_rrms()

◆ core()

◆ valence()

◆ hf_valence()

◆ basis()

◆ spectrum()

◆ CIwfs()

◆ CIwf()

◆ CI_integrals()

◆ vnuc()

◆ vHF()

◆ vlocal()

◆ Hmag()

◆ vrad()

◆ Sigma()

◆ Ncore()

◆ getState() [1/2]

◆ getState() [2/2]

◆ FermiLevel()

◆ energy_gap()

◆ coreConfiguration()

◆ coreConfiguration_nice()

◆ atom()

◆ atomicSymbol()

◆ identity()

◆ run_label()

◆ ion_degree()

◆ ion_symbol()

◆ Zion()

◆ printCore()

◆ printValence()

◆ printBasis()

◆ isInCore()

◆ isInValence()

◆ coreDensity()

◆ coreEnergyHF()

◆ set_HF()

◆ solve_core() [1/3]

◆ solve_core() [2/3]

◆ solve_core() [3/3]

◆ solve_valence()

◆ solve_exotic()

◆ hartreeFockBrueckner()

◆ fitSigma_hfBrueckner()

◆ radiativePotential() [1/2]

◆ radiativePotential() [2/2]

◆ formBasis()

◆ formSpectrum()

◆ formSigma()

◆ copySigma()

◆ update_Vnuc()

◆ lminmax_core_range()

◆ H0ab() [1/2]

◆ H0ab() [2/2]

◆ Hab()

◆ ConfigurationInteraction()

◆ output_to_json()