ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
Loading...
Searching...
No Matches
DiracSpinor.hpp
1#pragma once
2#include <memory>
3#include <string>
4#include <utility>
5#include <vector>
6class Grid;
7
42
43public:
46 DiracSpinor(int in_n, int in_kappa, std::shared_ptr<const Grid> in_rgrid);
47 using Index = uint16_t;
48
49private:
50 // Radial Grid; links F[i] to F(r)
51 std::shared_ptr<const Grid> m_rgrid;
52 // Principal quantum number
53 int m_n;
54 // Dirac quantum number, kappa
55 int m_kappa;
56 // Single-particle energy, not including rest energy
57 double m_en = 0.0;
58 // Upper (large) radial component
59 std::vector<double> m_f;
60 // Lower (small) radial component
61 std::vector<double> m_g;
62 // `practical zero': p0 is first non-zero point for f(r) [usually p0=0]
63 std::size_t m_p0 = 0;
64 // `practical infinity': pinf is last non-zero point for f(r)
65 std::size_t m_pinf;
66 // Number of iterations until energy convergence (for latest routine only)
67 int m_its = -1;
68 // Fractional energy convergence: eps = |(en'-en)/en|
69 double m_eps = -1.0;
70 // Occupation fraction. =1 for closed shells. =1/(2j+1) for valence
71 double m_occ_frac = 0.0;
72
73 // 2j, l, pi, kappa_index (for convenience): these are defined by n and kappa
74 int m_twoj;
75 int m_l;
76 int m_parity;
77 int m_kappa_index;
78 Index m_nkappa_index;
79
80 // flag for regular electron, or "exotic"
81 bool m_exotic{false};
82
83public:
85 int n() const { return m_n; }
87 int kappa() const { return m_kappa; }
89 int l() const { return m_l; }
91 double jjp1() const { return 0.25 * double(m_twoj * (m_twoj + 2)); }
92 int twoj() const { return m_twoj; }
94 int twojp1() const { return m_twoj + 1; }
96 int parity() const { return m_parity; }
98 int k_index() const { return m_kappa_index; }
100 Index nk_index() const { return m_nkappa_index; }
102 std::string symbol(bool gnuplot = false) const;
104 std::string shortSymbol() const;
105
107 bool &exotic() { return m_exotic; }
109 bool exotic() const { return m_exotic; }
110
112 void set_new_kappa(int new_kappa);
113
115 const Grid &grid() const { return *m_rgrid; };
118 std::shared_ptr<const Grid> grid_sptr() const { return m_rgrid; };
119
121 double en() const { return m_en; }
122 double &en() { return m_en; }
123
125 const std::vector<double> &f() const { return m_f; }
126 std::vector<double> &f() { return m_f; }
128 double f(std::size_t i) const { return m_f.at(i); }
129 double &f(std::size_t i) { return m_f.at(i); }
130
132 const std::vector<double> &g() const { return m_g; }
133 std::vector<double> &g() { return m_g; }
135 double g(std::size_t i) const { return m_g.at(i); }
136 double &g(std::size_t i) { return m_g.at(i); }
137
139 // XXX Kill this?
140 auto min_pt() const { return m_p0; }
141 auto &min_pt() { return m_p0; }
142
144 auto max_pt() const { return m_pinf; }
145 auto &max_pt() { return m_pinf; }
146
148 double r0() const;
150 double rinf() const;
151
153 auto its() const { return m_its; }
154 auto &its() { return m_its; }
155
157 auto occ_frac() const { return m_occ_frac; }
158 auto &occ_frac() { return m_occ_frac; }
159
161 auto eps() const { return m_eps; }
162 auto &eps() { return m_eps; }
163
165 const DiracSpinor &scale(const double factor);
167 const DiracSpinor &scale(const std::vector<double> &v);
169 void normalise(double norm_to = 1.0);
171 void zero_boundaries();
172
174 double norm() const;
176 double norm2() const;
177
179 std::pair<double, double> r0pinfratio() const;
180
182 std::vector<double> rho() const;
183
185 int num_electrons() const;
186
187public:
188 // Operator overloads
189
191 friend double operator*(const DiracSpinor &Fa, const DiracSpinor &Fb);
192
195 DiracSpinor &operator-=(const DiracSpinor &rhs);
196 friend DiracSpinor operator+(DiracSpinor lhs, const DiracSpinor &rhs);
197 friend DiracSpinor operator-(DiracSpinor lhs, const DiracSpinor &rhs);
198
200 DiracSpinor &operator*=(const double x);
201 friend DiracSpinor operator*(DiracSpinor Fa, const double x);
202 friend DiracSpinor operator*(const double x, DiracSpinor Fa);
203
205 DiracSpinor &operator*=(const std::vector<double> &v);
206 friend DiracSpinor operator*(const std::vector<double> &v, DiracSpinor Fa);
207
209 friend bool operator==(const DiracSpinor &lhs, const DiracSpinor &rhs);
210 friend bool operator!=(const DiracSpinor &lhs, const DiracSpinor &rhs);
211 friend bool operator<(const DiracSpinor &lhs, const DiracSpinor &rhs);
212 friend bool operator>(const DiracSpinor &lhs, const DiracSpinor &rhs);
213 friend bool operator<=(const DiracSpinor &lhs, const DiracSpinor &rhs);
214 friend bool operator>=(const DiracSpinor &lhs, const DiracSpinor &rhs);
215
217 static bool comp_l(const DiracSpinor &lhs, const DiracSpinor &rhs) {
218 return lhs.m_l < rhs.m_l;
219 }
220 static bool comp_j(const DiracSpinor &lhs, const DiracSpinor &rhs) {
221 return lhs.m_twoj < rhs.m_twoj;
222 }
223 static bool comp_ki(const DiracSpinor &lhs, const DiracSpinor &rhs) {
224 return lhs.m_kappa_index < rhs.m_kappa_index;
225 }
226 static bool comp_n(const DiracSpinor &lhs, const DiracSpinor &rhs) {
227 return lhs.m_n < rhs.m_n;
228 }
229 static bool comp_en(const DiracSpinor &lhs, const DiracSpinor &rhs) {
230 return lhs.en() < rhs.en();
231 }
232
234 friend std::ostream &operator<<(std::ostream &ostr, const DiracSpinor &Fa) {
235 return ostr << Fa.shortSymbol();
236 }
237
238 // Static (helper) functions:
239
241 static std::pair<double, std::string>
242 check_ortho(const std::vector<DiracSpinor> &a,
243 const std::vector<DiracSpinor> &b);
244
247 static void orthonormaliseOrbitals(std::vector<DiracSpinor> &orbs,
248 int num_its = 1);
249
251 [[nodiscard]] static DiracSpinor
252 orthogonaliseWrt(const DiracSpinor &Fin,
253 const std::vector<DiracSpinor> &orbs);
254
256 [[nodiscard]] static DiracSpinor
258 const std::vector<DiracSpinor> &orbs);
259
261 void orthog(const DiracSpinor &rhs);
262
264 static std::string shortSymbol(int n, int kappa);
265
267 static std::string state_config(const std::vector<DiracSpinor> &orbs);
268
270 static DiracSpinor exactHlike(int n, int k, std::shared_ptr<const Grid> rgrid,
271 double zeff, double alpha = 0.0);
272
274 static const DiracSpinor *find(int n, int k,
275 const std::vector<DiracSpinor> &orbs);
276
278 static int max_tj(const std::vector<DiracSpinor> &orbs);
280 static int max_l(const std::vector<DiracSpinor> &orbs);
282 static int max_n(const std::vector<DiracSpinor> &orbs);
284 static int max_kindex(const std::vector<DiracSpinor> &orbs);
285
292 static std::pair<std::vector<DiracSpinor>, std::vector<DiracSpinor>>
293 split_by_energy(const std::vector<DiracSpinor> &orbitals, double energy,
294 int n_min_core = 1);
295
302 static std::pair<std::vector<DiracSpinor>, std::vector<DiracSpinor>>
303 split_by_core(const std::vector<DiracSpinor> &orbitals,
304 const std::vector<DiracSpinor> &core, int n_min_core = 1);
305
309 static std::vector<DiracSpinor> subset(const std::vector<DiracSpinor> &basis,
310 const std::string &subset_string);
311};
Stores radial Dirac spinor: F_nk = (f, g)
Definition DiracSpinor.hpp:41
std::pair< double, double > r0pinfratio() const
Returns [f[p0]/f_max , f[pinf]/f_max] - for tests.
Definition DiracSpinor.cpp:106
void normalise(double norm_to=1.0)
By default normalises to 1, but can normalise to other number.
Definition DiracSpinor.cpp:91
static DiracSpinor orthonormaliseWrt(const DiracSpinor &Fin, const std::vector< DiracSpinor > &orbs)
As orthogonaliseWrt(), but normalises Fin afterwards.
Definition DiracSpinor.cpp:410
static int max_tj(const std::vector< DiracSpinor > &orbs)
Returns maximum (2j) found in {orbs}.
Definition DiracSpinor.cpp:266
int parity() const
(-1)^l, returns +/- 1
Definition DiracSpinor.hpp:96
double norm() const
Returns the norm, defined: Norm = Sqrt[<a|a>].
Definition DiracSpinor.cpp:65
static std::pair< double, std::string > check_ortho(const std::vector< DiracSpinor > &a, const std::vector< DiracSpinor > &b)
Returns worst |<a|b>| (or |<a|b>-1| for a=b) {val, state_names}.
Definition DiracSpinor.cpp:304
auto max_pt() const
Effective size(); index after last non-zero point (index for f[i])
Definition DiracSpinor.hpp:144
auto eps() const
Fractional energy convergence: eps = |(en'-en)/en|.
Definition DiracSpinor.hpp:161
const std::vector< double > & f() const
Upper (large) radial component function, f.
Definition DiracSpinor.hpp:125
int twojp1() const
2j+1
Definition DiracSpinor.hpp:94
int kappa() const
Dirac quantum number, kappa.
Definition DiracSpinor.hpp:87
static std::pair< std::vector< DiracSpinor >, std::vector< DiracSpinor > > split_by_core(const std::vector< DiracSpinor > &orbitals, const std::vector< DiracSpinor > &core, int n_min_core=1)
Splits orbitals into two groups (i.e., core, excited).
Definition DiracSpinor.cpp:455
double norm2() const
Returns the norm^2, defined: Norm2 = <a|a>
Definition DiracSpinor.cpp:66
auto occ_frac() const
Occupation fraction. =1 for closed shells. =1/(2j+1) for valence.
Definition DiracSpinor.hpp:157
const Grid & grid() const
Resturns a const reference to the radial grid.
Definition DiracSpinor.hpp:115
friend bool operator==(const DiracSpinor &lhs, const DiracSpinor &rhs)
Comparitor overloads (compares n, then kappa):
Definition DiracSpinor.cpp:206
int n() const
Principal quantum number, n.
Definition DiracSpinor.hpp:85
static DiracSpinor orthogonaliseWrt(const DiracSpinor &Fin, const std::vector< DiracSpinor > &orbs)
Forces Fin to be orthogonal to orbs. If Fn (i.e., {n,k}) in in orbs, replaces Fin with that from orbs...
Definition DiracSpinor.cpp:393
std::shared_ptr< const Grid > grid_sptr() const
Resturns copy of shared_ptr to grid [shared resource] - used when we want to construct a new DiracSpi...
Definition DiracSpinor.hpp:118
static int max_l(const std::vector< DiracSpinor > &orbs)
Returns maximum l found in {orbs}.
Definition DiracSpinor.cpp:272
std::string symbol(bool gnuplot=false) const
Single-electron term symbol (e.g., 6s_1/2). Gnuplot=true => 6s_{1/2}.
Definition DiracSpinor.cpp:32
DiracSpinor & operator*=(const double x)
Scalar multiplication.
Definition DiracSpinor.cpp:187
static int max_n(const std::vector< DiracSpinor > &orbs)
Returns maximum n found in {orbs}.
Definition DiracSpinor.cpp:277
double en() const
Single-particle energy, not including rest energy.
Definition DiracSpinor.hpp:121
Index nk_index() const
(n,kappa) index (see AtomData)
Definition DiracSpinor.hpp:100
void zero_boundaries()
Forces f(r) and g(r) to be zero outside of [p0,pinf)
Definition DiracSpinor.cpp:94
static const DiracSpinor * find(int n, int k, const std::vector< DiracSpinor > &orbs)
Searches for {n,k} in list of orbitals, returns pointer (may be null)
Definition DiracSpinor.cpp:289
double f(std::size_t i) const
Upper (large) radial component, f(r_i)
Definition DiracSpinor.hpp:128
double r0() const
r0 = r[min_pt] (in atomic units) XXX Kill this?
Definition DiracSpinor.cpp:135
double g(std::size_t i) const
Lower (small) radial component, g(r_i)
Definition DiracSpinor.hpp:135
void orthog(const DiracSpinor &rhs)
A more correct way to force orthogonality than normal.
Definition DiracSpinor.cpp:325
static std::vector< DiracSpinor > subset(const std::vector< DiracSpinor > &basis, const std::string &subset_string)
Takes a subset of an input basis (by copy), according to subset_string.
Definition DiracSpinor.cpp:475
double jjp1() const
j(j+1)
Definition DiracSpinor.hpp:91
std::string shortSymbol() const
e.g., 6p_1/2 => 6p-, 6p_3/2 => 6p+
Definition DiracSpinor.cpp:43
static DiracSpinor exactHlike(int n, int k, std::shared_ptr< const Grid > rgrid, double zeff, double alpha=0.0)
Constructs H-like (pointlike) DiracSpinor - mainly for testing.
Definition DiracSpinor.cpp:418
static bool comp_l(const DiracSpinor &lhs, const DiracSpinor &rhs)
Custom comparitors (for sorting): l, j, kappa_index, energy.
Definition DiracSpinor.hpp:217
bool & exotic()
Checks if spinor is for "exotic" lepton, or regular electron.
Definition DiracSpinor.hpp:107
auto min_pt() const
First non-zero point (index for f[i])
Definition DiracSpinor.hpp:140
static void orthonormaliseOrbitals(std::vector< DiracSpinor > &orbs, int num_its=1)
(approximately) OrthoNormalises a set of any orbitals.
Definition DiracSpinor.cpp:338
int k_index() const
kappa index (see AtomData)
Definition DiracSpinor.hpp:98
int num_electrons() const
Number of occupied electrons: (2j+1)*occ_frac.
Definition DiracSpinor.cpp:130
friend std::ostream & operator<<(std::ostream &ostr, const DiracSpinor &Fa)
Writes short symbol to ostream.
Definition DiracSpinor.hpp:234
const DiracSpinor & scale(const double factor)
Scales DiracSpinor (f and g) by constant.
Definition DiracSpinor.cpp:69
int l() const
Orbital angular momentum Q number.
Definition DiracSpinor.hpp:89
static int max_kindex(const std::vector< DiracSpinor > &orbs)
Returns maximum kappa_index found in {orbs}.
Definition DiracSpinor.cpp:282
DiracSpinor & operator+=(const DiracSpinor &rhs)
Addition of two DiracSpinors - must have same kappa.
Definition DiracSpinor.cpp:150
friend double operator*(const DiracSpinor &Fa, const DiracSpinor &Fb)
Returns radial integral (Fa,Fb) = Int(fa*fb + ga*gb)
Definition DiracSpinor.cpp:140
bool exotic() const
Checks if spinor is for "exotic" lepton, or regular electron.
Definition DiracSpinor.hpp:109
void set_new_kappa(int new_kappa)
Changes 'kappa' angular quantum number. Use with caution!
Definition DiracSpinor.cpp:55
static std::pair< std::vector< DiracSpinor >, std::vector< DiracSpinor > > split_by_energy(const std::vector< DiracSpinor > &orbitals, double energy, int n_min_core=1)
Splits orbitals into two groups (i.e., core, excited) by energy.
Definition DiracSpinor.cpp:437
double rinf() const
rinf = r[max_pt]
Definition DiracSpinor.cpp:136
std::vector< double > rho() const
rho(r) = sum_m |Psi^2|(r) = (2j+1) * x_occ * |F^2|(r)
Definition DiracSpinor.cpp:119
auto its() const
Number of iterations until energy convergence (for latest routine only)
Definition DiracSpinor.hpp:153
static std::string state_config(const std::vector< DiracSpinor > &orbs)
Returns formatted states string (e.g., '7sp5d') given list of orbs.
Definition DiracSpinor.cpp:233
const std::vector< double > & g() const
Lower (small) radial component function, g.
Definition DiracSpinor.hpp:132
Holds grid, including type + Jacobian (dr/du)
Definition Grid.hpp:31