ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
Loading...
Searching...
No Matches
Wavefunction.hpp
1#pragma once
2#include "CI/CSF.hpp"
3#include "HF/HartreeFock.hpp"
4#include "MBPT/CorrelationPotential.hpp"
5#include "Maths/Grid.hpp"
6#include "Physics/AtomData.hpp"
7#include "Physics/PhysConst_constants.hpp"
8#include "Potentials/NuclearPotentials.hpp"
9#include "Potentials/RadPot.hpp"
10#include "Wavefunction/BSplineBasis.hpp"
11#include "Wavefunction/DiracSpinor.hpp"
12#include "json/json.hpp"
13#include <iostream>
14#include <memory>
15#include <numeric>
16#include <optional>
17#include <string>
18#include <utility>
19#include <vector>
20
21namespace SplineBasis {
22struct Parameters;
23}
24
25//==============================================================================
38
39public:
43 Wavefunction(std::shared_ptr<const Grid> grid,
44 const Nuclear::Nucleus &nucleus, double var_alpha = 1.0,
45 const std::string &run_label = "");
47 Wavefunction(const GridParameters &gridparams,
48 const Nuclear::Nucleus &nucleus, double var_alpha = 1.0,
49 const std::string &run_label = "");
50
52
53private:
54 // Radial grid
55 std::shared_ptr<const Grid> rgrid;
56 // Internal value for alpha (alpha = var_alpha * alpha_0, alpha_0=~1/137)
57 double m_alpha;
58 std::string m_run_label;
59 // Holds nuclear parameters (isotope, charge distro etc.)
60 Nuclear::Nucleus m_nucleus;
61 // Valence (single-particle) orbitals
62 std::vector<DiracSpinor> m_valence{};
63 std::vector<DiracSpinor> m_hf_valence{};
64 // Basis, daigonalised over HF core. Used for MBPT
65 std::vector<DiracSpinor> m_basis{};
66 // Sprectrum: like basis, but includes Sigma (correlations).
67 std::vector<DiracSpinor> m_spectrum{};
68 // Nuclear potential // here AND hf?
69 std::vector<double> m_vnuc{};
70 // Hartree-Fock potential
71 std::optional<HF::HartreeFock> m_HF{std::nullopt};
72 // Correlation potential; for now unique_ptr; prefer std::optional
73 std::optional<MBPT::CorrelationPotential> m_Sigma{};
74 // Core configuration (non-rel terms)
75 std::string m_core_string = "";
76 std::string m_aboveFermi_core_string = "";
77
78 std::vector<CI::PsiJPi> m_CIwfs{};
79
80public:
82 const Grid &grid() const { return *rgrid; };
85 std::shared_ptr<const Grid> grid_sptr() const { return rgrid; };
86
88 double alpha() const { return m_alpha; }
89
91 double dalpha2() const {
92 // (alpha/alpha_0)^2 -1
93 return (m_alpha * m_alpha / PhysConst::alpha2) - 1.0;
94 }
95
97 const Nuclear::Nucleus &nucleus() const { return m_nucleus; }
99 int Znuc() const { return m_nucleus.z(); }
101 int Anuc() const { return m_nucleus.a(); }
103 double get_rrms() const { return m_nucleus.r_rms(); }
104
106 const std::vector<DiracSpinor> &core() const {
107 static const auto empty = std::vector<DiracSpinor>{}; //?
108 return m_HF ? m_HF->core() : empty;
109 }
110
112 const std::vector<DiracSpinor> &valence() const { return m_valence; }
113 std::vector<DiracSpinor> &valence() { return m_valence; }
114
115 const std::vector<DiracSpinor> &hf_valence() const { return m_hf_valence; }
116
119 const std::vector<DiracSpinor> &basis() const { return m_basis; }
120 std::vector<DiracSpinor> &basis() { return m_basis; }
121
123 const std::vector<DiracSpinor> &spectrum() const { return m_spectrum; }
124 std::vector<DiracSpinor> &spectrum() { return m_spectrum; }
125
126 const std::vector<CI::PsiJPi> &CIwfs() const { return m_CIwfs; }
127
128 const CI::PsiJPi *CIwf(int J, int parity) const {
129 for (const auto &ci_wf : m_CIwfs) {
130 if (ci_wf.twoJ() == 2 * J && ci_wf.parity() == parity)
131 return &ci_wf;
132 }
133 return nullptr;
134 }
135
138 const std::vector<double> &vnuc() const { return m_vnuc; }
139
141 const HF::HartreeFock *vHF() const { return m_HF ? &*m_HF : nullptr; }
142 HF::HartreeFock *vHF() { return m_HF ? &*m_HF : nullptr; }
143
146 std::vector<double> vlocal(int l = 0) const;
147
150 std::vector<double> Hmag(int l = 0) const;
151
153 const QED::RadPot *vrad() const { return m_HF ? m_HF->Vrad() : nullptr; }
154 QED::RadPot *vrad() { return m_HF ? m_HF->Vrad() : nullptr; }
155
157 const MBPT::CorrelationPotential *Sigma() const {
158 return m_Sigma ? &*m_Sigma : nullptr;
159 }
160 MBPT::CorrelationPotential *Sigma() { return m_Sigma ? &*m_Sigma : nullptr; }
161
162 //----------------------------------
163
165 int Ncore() const;
166
170 const DiracSpinor *getState(int n, int k) const;
172 const DiracSpinor *getState(std::string_view state) const;
173
175 double FermiLevel() const;
176
179 double energy_gap() const;
180
182 std::string coreConfiguration() const { return m_core_string; }
183
185 std::string coreConfiguration_nice() const {
186 return AtomData::niceCoreOutput(m_core_string);
187 }
188
190 std::string atom() const {
191 return AtomData::atomicSymbol(m_nucleus.z()) +
192 ", Z=" + std::to_string(m_nucleus.z()) +
193 " A=" + std::to_string(m_nucleus.a());
194 }
195
197 std::string atomicSymbol() const {
198 return AtomData::atomicSymbol(m_nucleus.z());
199 }
200
202 std::string identity() const;
203
205 int ion_degree(int num_val) const { return Zion() - num_val; }
206
208 std::string ion_symbol(int num_val) const {
209 return qip::int_to_roman(Zion() - num_val + 1);
210 }
211
213 int Zion() const { return Znuc() - Ncore(); }
214
216 void printCore() const;
219 void printValence(const std::vector<DiracSpinor> &tmp_orbitals = {}) const;
220
222 void printBasis(const std::vector<DiracSpinor> &the_basis) const;
223
225 bool isInCore(int n, int k) const;
226 bool isInValence(int n, int k) const;
227
229 std::vector<double> coreDensity() const;
230
232 double coreEnergyHF() const;
233
234 //------------------------------------------------------------------
235
238 void set_HF(const std::string &method = "HartreeFock",
239 const double x_Breit = 0.0, const std::string &in_core = "",
240 double eps_HF = 1.0e-13, bool print = true);
241
243 void solve_core(bool print = true);
244
246 void solve_core(const std::string &method, const double x_Breit = 0.0,
247 const std::string &in_core = "", double eps_HF = 1.0e-13,
248 bool print = true);
249
251 void solve_valence(const std::string &in_valence_str = "",
252 const bool print = true);
253
260 void solve_exotic(const std::string &in_exotic_str,
261 double mass = PhysConst::m_muon, bool print = true);
262
265 void hartreeFockBrueckner(const bool print = true);
266
268 void fitSigma_hfBrueckner(const std::string &valence_list,
269 const std::vector<double> &fit_energies);
270
272 void radiativePotential(QED::RadPot::Scale s, double rcut, double scale_rN,
273 const std::vector<double> &x_spd,
274 bool do_readwrite = true, bool print = true);
275
277 void radiativePotential(const IO::InputBlock &qed_input, bool do_readwrite,
278 bool print);
279
281 void formBasis(const SplineBasis::Parameters &params);
282
284 void formSpectrum(const SplineBasis::Parameters &params);
285
287 void formSigma(int nmin_core = 1, int nmin_core_F = 1, double r0 = 1.0e-4,
288 double rmax = 30.0, int stride = 4, bool each_valence = false,
289 bool include_G = false, bool include_Breit = false,
290 int n_max_breit = 0, const std::vector<double> &lambdas = {},
291 const std::vector<double> &fk = {},
292 const std::vector<double> &etak = {},
293 const std::string &in_fname = "",
294 const std::string &out_fname = "", bool FeynmanQ = false,
295 bool ScreeningQ = false, bool holeParticleQ = false,
296 int lmax = 6, double omre = -0.2, double w0 = 0.01,
297 double wratio = 1.5,
298 const std::optional<IO::InputBlock> &ek = std::nullopt);
299
300 // void correlations(const IO::InputBlock &input);
301
302 void copySigma(const MBPT::CorrelationPotential *const Sigma) {
303 if (Sigma != nullptr)
304 m_Sigma = *Sigma;
305 }
306
308 // _and_ HartreeFock)
309 void update_Vnuc(const std::vector<double> &v_new) {
311 m_vnuc = v_new;
312 if (m_HF) {
313 m_HF->vnuc() = v_new;
314 }
315 }
316
323 std::tuple<double, double> lminmax_core_range(int l, double eps = 0.0) const;
324
326 double H0ab(const DiracSpinor &Fa, const DiracSpinor &Fb) const;
328 double H0ab(const DiracSpinor &Fa, const DiracSpinor &dFa,
329 const DiracSpinor &Fb, const DiracSpinor &dFb) const;
330
331 double Hab(const DiracSpinor &Fa, const DiracSpinor &Fb) const;
332
334 void ConfigurationInteraction(const IO::InputBlock &input);
335
338 nlohmann::json
339 output_to_json(const std::string &out_name = "ampsci_output.json");
340
341private:
342 double H0ab_impl(const DiracSpinor &Fa, std::vector<double> dga,
343 const DiracSpinor &Fb, std::vector<double> dgb) const;
344
345 // Creates set of blank core orbitals
346 std::vector<DiracSpinor> determineCore(const std::string &str_core_in);
347 bool isInAboveFermiCore(int n, int k) const;
348};
Stores the CI Solutions for given J and parity (only two-electron).
Definition CSF.hpp:62
Stores radial Dirac spinor: F_nk = (f, g)
Definition DiracSpinor.hpp:41
Holds grid, including type + Jacobian (dr/du)
Definition Grid.hpp:31
Solves relativistic Hartree-Fock equations for core and valence. Optionally includes Breit and QED ef...
Definition HartreeFock.hpp:70
Holds list of Options, and a list of other InputBlocks. Can be initialised with a list of options,...
Definition InputBlock.hpp:142
Stores set of nuclear parameters (all radii in fm)
Definition NuclearPotentials.hpp:20
Class holds Flambaum-Ginges QED Radiative Potential.
Definition RadPot.hpp:15
Stores Wavefunction (set of valence orbitals, grid, HF etc.)
Definition Wavefunction.hpp:37
void fitSigma_hfBrueckner(const std::string &valence_list, const std::vector< double > &fit_energies)
First, fits Sigma to energies, then forms fitted Brueckner orbitals.
Definition Wavefunction.cpp:615
void solve_valence(const std::string &in_valence_str="", const bool print=true)
Performs hartree-Fock procedure for valence: note: poplulates valnece.
Definition Wavefunction.cpp:185
void formBasis(const SplineBasis::Parameters &params)
Calculates + populates basis [see BSplineBasis].
Definition Wavefunction.cpp:468
void solve_core(bool print=true)
Performs hartree-Fock procedure for core.
Definition Wavefunction.cpp:164
const std::vector< DiracSpinor > & core() const
Core orbitals (frozen HF core)
Definition Wavefunction.hpp:106
void printValence(const std::vector< DiracSpinor > &tmp_orbitals={}) const
Prints table of valence orbitals + energies etc.
Definition Wavefunction.cpp:407
const Grid & grid() const
Returns a const reference to the radial grid.
Definition Wavefunction.hpp:82
double coreEnergyHF() const
Calculates HF core energy (doesn't include magnetic QED?)
Definition Wavefunction.cpp:177
std::tuple< double, double > lminmax_core_range(int l, double eps=0.0) const
Returns [min,max] r values for which the core density (given l) is larger than cutoff (= eps*max_valu...
Definition Wavefunction.cpp:349
const DiracSpinor * getState(int n, int k) const
Finds requested state; returns nullptr if not found.
Definition Wavefunction.cpp:290
int Znuc() const
Nuclear charge, Z.
Definition Wavefunction.hpp:99
double energy_gap() const
Energy gap between lowest valence + highest core state: e(v) - e(c) [should be positive].
Definition Wavefunction.cpp:336
std::string identity() const
Atomic symbol, including core ionisation degree and run_label.
Definition Wavefunction.cpp:687
std::string coreConfiguration() const
Returns full core configuration.
Definition Wavefunction.hpp:182
nlohmann::json output_to_json(const std::string &out_name="ampsci_output.json")
Writes wavefunction information to json file; if out_name given, will print to that file.
Definition Wavefunction.cpp:918
int Zion() const
Effective charge (for core) = Z-N_core.
Definition Wavefunction.hpp:213
std::vector< double > vlocal(int l=0) const
Local part of potential, e.g., Vl = Vnuc + Vdir + Vrad_el(l) - can be l-dependent....
Definition Wavefunction.cpp:668
std::vector< double > Hmag(int l=0) const
QED Magnetic form factor. May return empty vector. Not typically l-dependent, but may be in future....
Definition Wavefunction.cpp:673
std::vector< double > coreDensity() const
Calculates rho(r) = sum_c psi^2(r) for core states, c={n,k,m}.
Definition Wavefunction.cpp:456
int ion_degree(int num_val) const
0 for neutral, 1 for singly-ionised etc.
Definition Wavefunction.hpp:205
int Anuc() const
Nuclear mass number, A.
Definition Wavefunction.hpp:101
std::string atomicSymbol() const
e.g., "Cs"
Definition Wavefunction.hpp:197
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.
Definition Wavefunction.cpp:220
std::string ion_symbol(int num_val) const
I for neutral, II for singly-ionised etc.
Definition Wavefunction.hpp:208
double get_rrms() const
Nuclear rms charge radii, in fm (femptometres)
Definition Wavefunction.hpp:103
const std::vector< DiracSpinor > & spectrum() const
Sprectrum: like basis, but includes correlations.
Definition Wavefunction.hpp:123
const MBPT::CorrelationPotential * Sigma() const
Returns ptr to (const) Correlation Potential, Sigma.
Definition Wavefunction.hpp:157
std::string coreConfiguration_nice() const
Returns core configuration, in nice output notation.
Definition Wavefunction.hpp:185
void formSpectrum(const SplineBasis::Parameters &params)
Calculates + populates Spectrum [see BSplineBasis].
Definition Wavefunction.cpp:487
void formSigma(int nmin_core=1, int nmin_core_F=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={}, const std::string &in_fname="", const std::string &out_fname="", bool FeynmanQ=false, bool ScreeningQ=false, bool holeParticleQ=false, int lmax=6, double omre=-0.2, double w0=0.01, double wratio=1.5, const std::optional< IO::InputBlock > &ek=std::nullopt)
Forms + stores correlation potential Sigma.
Definition Wavefunction.cpp:507
const std::vector< DiracSpinor > & valence() const
Valence orbitals (HF or Brueckner orbitals)
Definition Wavefunction.hpp:112
std::shared_ptr< const Grid > grid_sptr() const
Copy of shared_ptr to grid [shared resource] - used when we want to construct a new object that share...
Definition Wavefunction.hpp:85
void printCore() const
Prints table of core orbitals + energies etc.
Definition Wavefunction.cpp:375
bool isInCore(int n, int k) const
Check if a state is in the core (or valence) list.
Definition Wavefunction.cpp:261
std::string atom() const
String of atom info (e.g., "Cs, Z=55, A=133")
Definition Wavefunction.hpp:190
double alpha() const
Local value of fine-structure constant.
Definition Wavefunction.hpp:88
double H0ab(const DiracSpinor &Fa, const DiracSpinor &Fb) const
Returns <a|H|b> for Hamiltonian H (inludes Rad.pot, NOT sigma, Breit, or exchange!...
Definition Wavefunction.cpp:718
const HF::HartreeFock * vHF() const
Returns ptr to Hartree Fock (class)
Definition Wavefunction.hpp:141
void ConfigurationInteraction(const IO::InputBlock &input)
Runs the CI+MBPT routines; stores wavefunctions.
Definition Wavefunction.cpp:778
void update_Vnuc(const std::vector< double > &v_new)
Allows extra potential to be added to Vnuc (updates both in Wavefunction.
Definition Wavefunction.hpp:309
void set_HF(const std::string &method="HartreeFock", const double x_Breit=0.0, const std::string &in_core="", double eps_HF=1.0e-13, bool print=true)
Initialises HF object and populates core orbitals (does not solve HF equations)
Definition Wavefunction.cpp:122
const std::vector< double > & vnuc() const
Nuclear potential. Only provide const version, since HF and WF version of vnuc must be kept in sync.
Definition Wavefunction.hpp:138
void solve_exotic(const std::string &in_exotic_str, double mass=PhysConst::m_muon, bool print=true)
Solves for exotic atoms (e.g., muonic), including screening. Resulting states are included in valence...
Definition Wavefunction.cpp:787
void hartreeFockBrueckner(const bool print=true)
Forms Bruckner valence orbitals: (H_hf + Sigma)|nk> = e|nk>. Replaces existing valence states.
Definition Wavefunction.cpp:605
int Ncore() const
Number of electrons in the core.
Definition Wavefunction.cpp:679
const Nuclear::Nucleus & nucleus() const
Returns Nuclear::nucleus object (contains nuc. parameters)
Definition Wavefunction.hpp:97
const std::vector< DiracSpinor > & basis() const
Basis, eigenstates of HF potential. Used for MBPT. Includes Breit and QED (if they are included),...
Definition Wavefunction.hpp:119
const QED::RadPot * vrad() const
Pointer to QED radiative potnential. May be nullptr.
Definition Wavefunction.hpp:153
double FermiLevel() const
Returns energy location of the "Fermi Level", - energy half way between core/valence....
Definition Wavefunction.cpp:317
void printBasis(const std::vector< DiracSpinor > &the_basis) const
Prints table of Basis/Spectrum orbitals, compares to HF orbitals.
Definition Wavefunction.cpp:433
double dalpha2() const
Variation in alpha^2 : x = (alpha/alpha_0)^2 - 1.
Definition Wavefunction.hpp:91
std::string niceCoreOutput(const std::string &full_core)
Given a full electron config., returns nicer format by recognising nobel gas.
Definition AtomData.cpp:198
std::string atomicSymbol(int Z)
e.g., 55 -> "Cs"
Definition AtomData.cpp:64
constexpr double m_muon
muon mass in atomic units (m_mu/m_e). Codata 2022: 206.768 2827(46)
Definition PhysConst_constants.hpp:39
Constucts of spinor/orbital basis using B-splines (DKB/Reno/Derevianko-Beloy method)
Definition BSplineBasis.cpp:20
std::string int_to_roman(int a)
Converts integer, a, to Roman Numerals. Assumed that |a|<=4000.
Definition String.hpp:321
Parmaters used to create Grid.
Definition Grid.hpp:9
Scale factors for Uehling, high, low, magnetic, Wickman-Kroll.
Definition RadPot.hpp:20