ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
Loading...
Searching...
No Matches
jL.hpp
1#pragma once
2#include "DiracOperator/TensorOperator.hpp"
3#include "Maths/Grid.hpp"
4#include "Maths/SphericalBessel.hpp"
5#include "Wavefunction/Wavefunction.hpp"
6
7namespace DiracOperator {
8
10
19class jL : public TensorOperator {
20public:
22 jL(const Grid &r_grid, const Grid &q_grid, std::size_t max_l,
23 bool subtract_one = false)
24 : TensorOperator(-1, Parity::blank),
25 m_max_l(max_l),
26 m_j_lq_r(m_max_l + 1, q_grid.num_points()),
27 m_r_grid(&r_grid),
28 m_q_grid(&q_grid),
29 m_subtract_one(subtract_one) {
30 fill_table();
31 }
33 jL(const jL &other)
34 : TensorOperator(-1, Parity::blank),
35 m_max_l(other.m_max_l),
36 m_j_lq_r(other.m_j_lq_r),
37 m_r_grid(other.m_r_grid),
38 m_q_grid(other.m_q_grid),
39 m_subtract_one(other.m_subtract_one) {}
40
41 // jL(const jL &) = default;
42 jL &operator=(const jL &) = delete;
43
44protected:
45 std::size_t m_max_l;
47 const Grid *m_r_grid;
48 const Grid *m_q_grid;
49 bool m_subtract_one;
50
51private:
52 // fills lookup J_l_q_r table
53 void fill_table() {
54 std::cout << "Filling jL lookup table: " << std::flush;
55 for (std::size_t l = 0; l <= m_max_l; ++l) {
56#pragma omp parallel for
57 for (std::size_t iq = 0; iq < m_q_grid->num_points(); ++iq) {
58 m_j_lq_r.at(l, iq) = SphericalBessel::fillBesselVec_kr(
59 int(l), m_q_grid->r(iq), m_r_grid->r());
60 }
61 }
62 std::cout << "done.\n";
63 }
64
65public:
67 std::size_t L() const { return std::size_t(m_rank); }
69 std::size_t max_L() const { return m_max_l; }
70
71 const auto &q_grid() const { return *m_q_grid; }
72 const auto &r_grid() const { return *m_r_grid; }
73
75 virtual void set_L_q(std::size_t L, double q) {
76 assert(L <= m_max_l && "L must be <= max L");
77 m_rank = int(L);
78 m_parity = Angular::evenQ(m_rank) ? Parity::even : Parity::odd;
79 const auto iq = m_q_grid->getIndex(q);
80 m_vec = m_j_lq_r.at(L, iq);
81 }
82
83 // This is _not_ thread safe
84 DiracSpinor radial_rhs(const int kappa_a,
85 const DiracSpinor &Fb) const override {
86 // Fa * radial_rhs(Fa.kappa(),Fb) = h.radialIntegral(Fa, Fb)
87
88 // const auto &gr = Fb.grid();
89 DiracSpinor dF(0, kappa_a, Fb.grid_sptr());
90 dF.min_pt() = Fb.min_pt();
91 dF.max_pt() = Fb.max_pt();
92
93 if (isZero(kappa_a, Fb.kappa())) {
94 dF.min_pt() = Fb.min_pt();
95 dF.max_pt() = Fb.min_pt();
96 return dF;
97 }
98
99 const auto &df = Fb.f();
100 const auto &dg = Fb.g();
101
102 const auto cff = angularCff(kappa_a, Fb.kappa());
103 const auto cgg = angularCgg(kappa_a, Fb.kappa());
104 const auto cfg = angularCfg(kappa_a, Fb.kappa());
105 const auto cgf = angularCgf(kappa_a, Fb.kappa());
106
107 if (m_subtract_one && m_rank == 0 && Fb.kappa() == kappa_a) {
108 // jL -> (jL-1)
109 // Vector
110 // (ff + gg)*jL - (ff+gg)
111 // = (ff + gg)*(jL-1)
112
113 // Scalar:
114 // (ff - gg)*jL - (ff+gg)
115 // = ff*(jL-1) - gg*(jL+1)
116 if (cgg > 0.0) {
117 //vector
118 for (auto i = Fb.min_pt(); i < Fb.max_pt(); i++) {
119 dF.f(i) = m_constant * (m_vec[i] - 1.0) * cff * df[i];
120 dF.g(i) = m_constant * (m_vec[i] - 1.0) * cgg * dg[i];
121 }
122 } else if (cgg < 0.0) {
123 //scalar
124 for (auto i = Fb.min_pt(); i < Fb.max_pt(); i++) {
125 dF.f(i) = m_constant * (m_vec[i] - 1.0) * cff * df[i];
126 dF.g(i) = m_constant * (m_vec[i] + 1.0) * cgg * dg[i];
127 }
128 } else {
129 assert(false && "Error 111 - subtract_one for pseudo case?");
130 }
131 } else {
132 for (auto i = Fb.min_pt(); i < Fb.max_pt(); i++) {
133 dF.f(i) = m_constant * m_vec[i] * (cff * df[i] + cfg * dg[i]);
134 dF.g(i) = m_constant * m_vec[i] * (cgf * df[i] + cgg * dg[i]);
135 }
136 }
137
138 return dF;
139 }
140
141 virtual double radialIntegral(const DiracSpinor &Fa,
142 const DiracSpinor &Fb) const override final {
143 // nb: faster not to do this, but nicer this way
144 return Fa * radial_rhs(Fa.kappa(), Fb);
145 }
146
148 double rme(const DiracSpinor &a, const DiracSpinor &b, std::size_t L,
149 double q) const {
150 const auto iq = m_q_grid->getIndex(q);
151 const auto &jlqr = m_j_lq_r.at(L, iq);
152 const auto &gr = *m_r_grid;
153 const auto cff = angularCff(a.kappa(), b.kappa());
154 const auto cgg = angularCgg(a.kappa(), b.kappa());
155 const auto cfg = angularCfg(a.kappa(), b.kappa());
156 const auto cgf = angularCgf(a.kappa(), b.kappa());
157 auto max = std::min(a.max_pt(), b.max_pt());
158
159 double Rf{0.0}, Rg{0.0};
160 if (m_subtract_one && a.kappa() == b.kappa() && L == 0) {
161 auto jl_neg1 = jlqr;
162 for (auto &el : jl_neg1) {
163 el -= 1.0;
164 }
165 // jL -> (jL-1)
166 // Vector
167 // (ff + gg)*jL - (ff+gg)
168 // = (ff + gg)*(jL-1)
169
170 // Scalar:
171 // (ff - gg)*jL - (ff+gg)
172 // (ff - gg)*jL - (ff-gg)
173 // = ff*(jL-1) - gg*(jL+1)
174
175 Rf = NumCalc::integrate(1.0, 0, max, a.f(), b.f(), jl_neg1, gr.drdu());
176
177 if (cgg > 0.01) {
178 //vector
179 Rg = NumCalc::integrate(1.0, 0, max, a.g(), b.g(), jl_neg1, gr.drdu());
180 } else if (cgg < -0.01) {
181 //scalar (-ve sign is inside cgg)
182 auto jl_p1 = jlqr;
183 for (auto &el : jl_p1) {
184 el += 1.0;
185 }
186 Rg = NumCalc::integrate(1.0, 0, max, a.g(), b.g(), jl_p1, gr.drdu());
187 } else {
188 assert(false && "Error: subtract 1 in pseudo-case?");
189 }
190 } else {
191 Rf = cff == 0.0 ?
192 0.0 :
193 NumCalc::integrate(1.0, 0, max, a.f(), b.f(), jlqr, gr.drdu());
194 Rg = cgg == 0.0 ?
195 0.0 :
196 NumCalc::integrate(1.0, 0, max, a.g(), b.g(), jlqr, gr.drdu());
197 }
198 const auto Rfg = cfg == 0.0 ? 0.0 :
199 NumCalc::integrate(1.0, 0, max, a.f(), b.g(),
200 jlqr, gr.drdu());
201 const auto Rgf = cgf == 0.0 ? 0.0 :
202 NumCalc::integrate(1.0, 0, max, a.g(), b.f(),
203 jlqr, gr.drdu());
204
205 const auto kb = cfg == 0.0 ? b.kappa() : -b.kappa();
206 const auto s = cfg == 0.0 ? 1.0 : -1.0;
207 return s * Angular::Ck_kk(int(L), a.kappa(), kb) *
208 (cff * Rf + cgg * Rg + cfg * Rfg + cgf * Rgf) * gr.du();
209 }
210
212 bool is_zero(const DiracSpinor &a, const DiracSpinor &b,
213 std::size_t L) const {
214 const auto kb =
215 angularCfg(a.kappa(), b.kappa()) == 0.0 ? b.kappa() : -b.kappa();
216 return !Angular::Ck_kk_SR(int(L), a.kappa(), kb);
217 }
218
219public:
220 virtual double angularF(const int ka, const int kb) const override {
221 return Angular::Ck_kk(m_rank, ka, kb);
222 }
223 virtual double angularCff(int, int) const override { return 1.0; }
224 virtual double angularCgg(int, int) const override { return 1.0; }
225 virtual double angularCfg(int, int) const override { return 0.0; }
226 virtual double angularCgf(int, int) const override { return 0.0; }
227
228 virtual std::string name() const override { return std::string("jL"); }
229 std::string units() const override final { return "au"; }
230};
231
232//------------------------------------------------------------------------------
234class g0jL : public jL {
235public:
236 g0jL(const Grid &r_grid, const Grid &q_grid, std::size_t max_l,
237 bool subtract_one = false)
238 : jL(r_grid, q_grid, max_l, subtract_one) {}
239 g0jL(const jL &other) : jL(other) {}
240
241public:
242 double angularCff(int, int) const override final { return 1.0; }
243 double angularCgg(int, int) const override final { return -1.0; }
244 double angularCfg(int, int) const override final { return 0.0; }
245 double angularCgf(int, int) const override final { return 0.0; }
246
247 std::string name() const override final { return std::string("g0jL"); }
248};
249
250//------------------------------------------------------------------------------
252class ig5jL : public jL {
253public:
254 ig5jL(const Grid &r_grid, const Grid &q_grid, std::size_t max_l)
255 : jL(r_grid, q_grid, max_l, false) {}
256 ig5jL(const jL &other) : jL(other) {}
257
258public:
259 void set_L_q(std::size_t L, double q) override final {
260 assert(L <= m_max_l && "L must be <= max L");
261 m_rank = int(L);
262 // Note: opposite parity for 'pseudo' cases
263 m_parity = Angular::evenQ(m_rank) ? Parity::odd : Parity::even;
264 const auto iq = m_q_grid->getIndex(q);
265 m_vec = m_j_lq_r.at(L, iq);
266 }
267
268 double angularF(const int ka, const int kb) const override final {
269 return -1.0 * Angular::Ck_kk(m_rank, ka, -kb);
270 }
271
272 double angularCff(int, int) const override final { return 0.0; }
273 double angularCgg(int, int) const override final { return 0.0; }
274 double angularCfg(int, int) const override final { return 1.0; }
275 double angularCgf(int, int) const override final { return -1.0; }
276
277 std::string name() const override final { return std::string("ig5jL"); }
278};
279
280//------------------------------------------------------------------------------
282class ig0g5jL : public jL {
283public:
284 ig0g5jL(const Grid &r_grid, const Grid &q_grid, std::size_t max_l)
285 : jL(r_grid, q_grid, max_l, false) {}
286 ig0g5jL(const jL &other) : jL(other) {}
287
288public:
289 void set_L_q(std::size_t L, double q) override final {
290 assert(L <= m_max_l && "L must be <= max L");
291 m_rank = int(L);
292 // Note: opposite parity for 'pseudo' cases
293 m_parity = Angular::evenQ(m_rank) ? Parity::odd : Parity::even;
294 const auto iq = m_q_grid->getIndex(q);
295 m_vec = m_j_lq_r.at(L, iq);
296 }
297
298 virtual double angularF(const int ka, const int kb) const override {
299 return -1.0 * Angular::Ck_kk(m_rank, ka, -kb);
300 }
301
302 virtual double angularCff(int, int) const override { return 0.0; }
303 virtual double angularCgg(int, int) const override { return 0.0; }
304 virtual double angularCfg(int, int) const override { return 1.0; }
305 virtual double angularCgf(int, int) const override { return 1.0; }
306
307 std::string name() const override final { return std::string("ig0g5jL"); }
308};
309
310} // namespace DiracOperator
General operator (virtual base class); operators derive from this.
Definition TensorOperator.hpp:110
bool isZero(const int ka, int kb) const
If matrix element <a|h|b> is zero, returns true.
Definition TensorOperator.cpp:15
Matrix element of tensor operator: gamma^0 J_L(qr) C^L.
Definition jL.hpp:234
std::string name() const override final
Returns "name" of operator (e.g., 'E1')
Definition jL.hpp:247
Matrix element of tensor operator: i gamma^0gamma^5 J_L(qr) C^L. nb: i makes ME real.
Definition jL.hpp:282
virtual double angularF(const int ka, const int kb) const override
angularF: links radiation integral to RME. RME = <a||h||b> = angularF(a,b) * radial_int(a,...
Definition jL.hpp:298
void set_L_q(std::size_t L, double q) override final
Sets the current L and q values for use. Note: NOT thread safe!
Definition jL.hpp:289
std::string name() const override final
Returns "name" of operator (e.g., 'E1')
Definition jL.hpp:307
Matrix element of tensor operator: i gamma^5 J_L(qr) C^L. nb: i makes ME real.
Definition jL.hpp:252
std::string name() const override final
Returns "name" of operator (e.g., 'E1')
Definition jL.hpp:277
double angularF(const int ka, const int kb) const override final
angularF: links radiation integral to RME. RME = <a||h||b> = angularF(a,b) * radial_int(a,...
Definition jL.hpp:268
void set_L_q(std::size_t L, double q) override final
Sets the current L and q values for use. Note: NOT thread safe!
Definition jL.hpp:259
Matrix element of tensor operator: J_L(qr)*C^L.
Definition jL.hpp:19
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.
Definition jL.hpp:22
virtual std::string name() const override
Returns "name" of operator (e.g., 'E1')
Definition jL.hpp:228
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.
Definition jL.hpp:148
virtual double radialIntegral(const DiracSpinor &Fa, const DiracSpinor &Fb) const override final
Defined via <a||h||b> = angularF(a,b) * radialIntegral(a,b) (Note: if radial_rhs is overridden,...
Definition jL.hpp:141
std::string units() const override final
Returns units of operator (usually au, may be MHz, etc.)
Definition jL.hpp:229
std::size_t max_L() const
Maximum L value in table.
Definition jL.hpp:69
virtual double angularF(const int ka, const int kb) const override
angularF: links radiation integral to RME. RME = <a||h||b> = angularF(a,b) * radial_int(a,...
Definition jL.hpp:220
virtual void set_L_q(std::size_t L, double q)
Sets the current L and q values for use. Note: NOT thread safe!
Definition jL.hpp:75
DiracSpinor radial_rhs(const int kappa_a, const DiracSpinor &Fb) const override
radial_int = Fa * radial_rhs(a, Fb) (a needed for angular factor)
Definition jL.hpp:84
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)
Definition jL.hpp:212
std::size_t L() const
Current value of L (should = rank)
Definition jL.hpp:67
jL(const jL &other)
Constructing from existing operator: copies JL table - faster. Can copy from a jL of different type (...
Definition jL.hpp:33
l (orbital angular momentum) operator
Definition jls.hpp:27
s (spin) operator
Definition jls.hpp:50
Stores radial Dirac spinor: F_nk = (f, g)
Definition DiracSpinor.hpp:41
auto max_pt() const
Effective size(); index after last non-zero point (index for f[i])
Definition DiracSpinor.hpp:136
const std::vector< double > & f() const
Upper (large) radial component function, f.
Definition DiracSpinor.hpp:117
int kappa() const
Dirac quantum number, kappa.
Definition DiracSpinor.hpp:84
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:110
auto min_pt() const
First non-zero point (index for f[i])
Definition DiracSpinor.hpp:132
const std::vector< double > & g() const
Lower (small) radial component function, g.
Definition DiracSpinor.hpp:124
Holds grid, including type + Jacobian (dr/du)
Definition Grid.hpp:31
const std::vector< double > & r() const
Grid points, r.
Definition Grid.hpp:75
auto num_points() const
Number of grid points.
Definition Grid.hpp:64
std::size_t getIndex(double x, bool require_nearest=false) const
Given value, returns grid index. if not require_nearest, returns lower.
Definition Grid.cpp:76
Matrix class; row-major.
Definition Matrix.hpp:35
T & at(std::size_t row_i, std::size_t col_j)
() index access (with range checking). (i,j) returns ith row, jth col
Definition Matrix.hpp:109
bool Ck_kk_SR(int k, int ka, int kb)
Ck selection rule only. Returns false if C^k=0, true if non-zero.
Definition Wigner369j.hpp:324
constexpr bool evenQ(int a)
Returns true if a is even - for integer values.
Definition Wigner369j.hpp:146
double Ck_kk(int k, int ka, int kb)
Reduced (relativistic) angular ME: <ka||C^k||kb> [takes k and kappa].
Definition Wigner369j.hpp:298
Dirac Operators: General + derived.
Definition GenerateOperator.cpp:12