ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
Loading...
Searching...
No Matches
TensorOperator.hpp
1#pragma once
2#include "Angular/Wigner369j.hpp"
3#include "Maths/Grid.hpp"
4#include "Maths/NumCalc_quadIntegrate.hpp"
5#include "Maths/SphericalBessel.hpp"
6#include "Physics/PhysConst_constants.hpp"
7#include "Potentials/NuclearPotentials.hpp"
8#include "Wavefunction/DiracSpinor.hpp"
9#include "qip/Vector.hpp"
10#include <cmath>
11#include <string>
12#include <vector>
13
15namespace DiracOperator {
16
17enum class Parity { even, odd, blank };
18enum class Realness { real, imaginary, blank };
19
20//==============================================================================
22struct IntM4x4
23// 4x4 Integer matrix.
24// Notation for elements:
25// (e00 e01)
26// (e10 e11)
27{
28 IntM4x4(int in00, int in01, int in10, int in11)
29 : e00(in00), e01(in01), e10(in10), e11(in11) {}
30
31 const int e00, e01, e10, e11;
32};
33
34//==============================================================================
35class SpinorMatrix
36// 2x2 matrix that acts in radial spinor space.
37// Notation for elements:
38// (ff fg)
39// (gf gg)
40{
41 double ff{0.0}, fg{0.0}, gf{0.0}, gg{0.0};
42
43public:
44 // SpinorMatrix() = default;
45 constexpr SpinorMatrix(double iff, double ifg, double igf, double igg)
46 : ff(iff), fg(ifg), gf(igf), gg(igg) {}
47 constexpr SpinorMatrix() {}
48
49 constexpr SpinorMatrix &operator+=(const SpinorMatrix &rhs) {
50 ff += rhs.ff;
51 fg += rhs.fg;
52 gf += rhs.gf;
53 gg += rhs.gg;
54 return *this;
55 }
56 constexpr SpinorMatrix &operator-=(const SpinorMatrix &rhs) {
57 ff -= rhs.ff;
58 fg -= rhs.fg;
59 gf -= rhs.gf;
60 gg -= rhs.gg;
61 return *this;
62 }
63 friend SpinorMatrix operator+(SpinorMatrix lhs, const SpinorMatrix &rhs) {
64 return lhs += rhs;
65 }
66 friend SpinorMatrix operator-(SpinorMatrix lhs, const SpinorMatrix &rhs) {
67 return lhs -= rhs;
68 }
69 constexpr SpinorMatrix &operator*=(double x) {
70 ff *= x;
71 fg *= x;
72 gf *= x;
73 gg *= x;
74 return *this;
75 }
76 friend SpinorMatrix operator*(SpinorMatrix lhs, double x) { return lhs *= x; }
77 friend SpinorMatrix operator*(double x, SpinorMatrix rhs) { return rhs *= x; }
78
79 DiracSpinor act(const DiracSpinor &Fa) const {
80 using namespace qip::overloads;
81 auto mFa = Fa;
82 // mFa.f() = ff * Fa.f() + fg * Fa.g();
83 // mFa.g() = gf * Fa.f() + gg * Fa.g();
84
85 // since we start with a copy:
86 if (ff != 1.0)
87 mFa.f() *= ff;
88 if (fg != 0.0)
89 mFa.f() += fg * Fa.g();
90 if (gg != 1.0)
91 mFa.g() *= gg;
92 if (gf != 0.0)
93 mFa.g() += gf * Fa.f();
94 return mFa;
95 }
96 DiracSpinor operator*(const DiracSpinor &Fa) const { return act(Fa); }
97};
98
99//==============================================================================
111protected:
112 TensorOperator(int rank_k, Parity pi, double constant = 1.0,
113 const std::vector<double> &inv = {}, int diff_order = 0,
114 Realness RorI = Realness::real, bool freq_dep = false)
115 : m_rank(rank_k),
116 m_parity(pi),
117 m_diff_order(diff_order),
118 opC(RorI),
119 m_freqDependantQ(freq_dep),
120 m_constant(constant),
121 m_vec(inv){};
122
123public:
124 virtual ~TensorOperator() = default;
125
126protected:
127 int m_rank;
128 Parity m_parity;
129 int m_diff_order;
130 Realness opC;
131 bool m_freqDependantQ{false};
132
133protected:
134 // these may be updated for frequency-dependant operators
135 double m_constant; // included in radial integral
136 std::vector<double> m_vec;
137
138public:
139 bool freqDependantQ() const { return m_freqDependantQ; }
140
141public:
143 bool isZero(const int ka, int kb) const;
144 bool isZero(const DiracSpinor &Fa, const DiracSpinor &Fb) const;
145
146 bool selectrion_rule(int twoJA, int piA, int twoJB, int piB) const {
147 if (twoJA == twoJB && twoJA == 0.0)
148 return false;
149
150 if (Angular::triangle(twoJA, twoJB, 2 * m_rank) == 0)
151 return false;
152
153 // if (2 * m_rank < std::abs(twoJA - twoJB))
154 // return false;
155 return (m_parity == Parity::even) == (piA == piB);
156 }
157
159 virtual void updateFrequency(const double){};
160
162 void scale(double lambda);
163
165 const std::vector<double> &getv() const { return m_vec; }
167 double getc() const { return m_constant; }
168 int get_d_order() const { return m_diff_order; }
169
171 bool imaginaryQ() const { return (opC == Realness::imaginary); }
172 int rank() const { return m_rank; }
174 int parity() const { return (m_parity == Parity::even) ? 1 : -1; }
175
177 int symm_sign(const DiracSpinor &Fa, const DiracSpinor &Fb) const {
178 const auto sra_i = imaginaryQ() ? -1 : 1;
179 const auto sra = Angular::neg1pow_2(Fa.twoj() - Fb.twoj());
180 return sra_i * sra;
181 }
182
184 virtual std::string name() const { return "Operator"; };
186 virtual std::string units() const { return "au"; };
187
188public:
189 // These are needed for radial integrals
190 // Usually just constants, but can also be functions of kappa
191 virtual double angularCff(int /*k_a*/, int /*k_b*/) const { return 1.0; }
192 virtual double angularCgg(int, int) const { return 1.0; }
193 virtual double angularCfg(int, int) const { return 0.0; }
194 virtual double angularCgf(int, int) const { return 0.0; }
195
196public:
199 virtual double angularF(const int, const int) const = 0;
201 virtual DiracSpinor radial_rhs(const int kappa_a,
202 const DiracSpinor &Fb) const;
203
206 virtual double radialIntegral(const DiracSpinor &Fa,
207 const DiracSpinor &Fb) const;
208
210 double rme3js(const int twoja, const int twojb, int two_mb = 1,
211 int two_q = 0) const;
212
214 DiracSpinor reduced_rhs(const int ka, const DiracSpinor &Fb) const;
215
217 DiracSpinor reduced_lhs(const int ka, const DiracSpinor &Fb) const;
218
220 double reducedME(const DiracSpinor &Fa, const DiracSpinor &Fb) const;
221
224 double fullME(const DiracSpinor &Fa, const DiracSpinor &Fb,
225 std::optional<int> two_ma = std::nullopt,
226 std::optional<int> two_mb = std::nullopt,
227 std::optional<int> two_q = std::nullopt) const;
228};
229
230//============================================================================
231//============================================================================
234public:
235 ScalarOperator(Parity pi, double in_coef,
236 const std::vector<double> &in_v = {},
237 const IntM4x4 &in_g = {1, 0, 0, 1}, int in_diff = 0,
238 Realness rori = Realness::real)
239 : TensorOperator(0, pi, in_coef, in_v, in_diff, rori),
240 c_ff(in_g.e00),
241 c_fg(in_g.e01),
242 c_gf(in_g.e10),
243 c_gg(in_g.e11) {}
244
245 ScalarOperator(const std::vector<double> &in_v)
246 : TensorOperator(0, Parity::even, 1.0, in_v, 0),
247 c_ff(1.0),
248 c_fg(0.0),
249 c_gf(0.0),
250 c_gg(1.0) {}
251
252 ScalarOperator(double in_coef, const std::vector<double> &in_v = {})
253 : TensorOperator(0, Parity::even, in_coef, in_v, 0),
254 c_ff(1.0),
255 c_fg(0.0),
256 c_gf(0.0),
257 c_gg(1.0) {}
258
259public:
260 virtual double angularF(const int ka, const int kb) const override {
261 // For scalar operators, <a||h||b> = RadInt / 3js
262 // 3js:= 1/(Sqrt[2j+1]) ... depends on m???
263 // |k| = j+1/2
264 return (std::abs(ka) == std::abs(kb)) ? std::sqrt(2.0 * std::abs(ka)) : 0.0;
265 }
266
267private:
268 const double c_ff, c_fg, c_gf, c_gg;
269
270protected:
271 double virtual angularCff(int, int) const override { return c_ff; }
272 double virtual angularCgg(int, int) const override { return c_gg; }
273 double virtual angularCfg(int, int) const override { return c_fg; }
274 double virtual angularCgf(int, int) const override { return c_gf; }
275};
276
277//------------------------------------------------------------------------------
279class NullOperator final : public ScalarOperator {
280public:
281 NullOperator() : ScalarOperator(Parity::even, 0, {}) {}
282
283protected:
284 double virtual angularCff(int, int) const override final { return 0.0; }
285 double virtual angularCgg(int, int) const override final { return 0.0; }
286 double virtual angularCfg(int, int) const override final { return 0.0; }
287 double virtual angularCgf(int, int) const override final { return 0.0; }
288};
289
290} // namespace DiracOperator
Speacial operator: 0.
Definition TensorOperator.hpp:279
Speacial case for scalar operator.
Definition TensorOperator.hpp:233
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 TensorOperator.hpp:260
General operator (virtual base class); operators derive from this.
Definition TensorOperator.hpp:110
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,...
Definition TensorOperator.cpp:67
int symm_sign(const DiracSpinor &Fa, const DiracSpinor &Fb) const
returns relative sign between <a||x||b> and <b||x||a>
Definition TensorOperator.hpp:177
bool imaginaryQ() const
returns true if operator is imaginary (has imag MEs)
Definition TensorOperator.hpp:171
virtual std::string units() const
Returns units of operator (usually au, may be MHz, etc.)
Definition TensorOperator.hpp:186
int parity() const
returns parity, as integer (+1 or -1)
Definition TensorOperator.hpp:174
virtual double angularF(const int, const int) const =0
angularF: links radiation integral to RME. RME = <a||h||b> = angularF(a,b) * radial_int(a,...
bool isZero(const int ka, int kb) const
If matrix element <a|h|b> is zero, returns true.
Definition TensorOperator.cpp:15
virtual std::string name() const
Returns "name" of operator (e.g., 'E1')
Definition TensorOperator.hpp:184
double getc() const
Returns a const ref to constant c.
Definition TensorOperator.hpp:167
double reducedME(const DiracSpinor &Fa, const DiracSpinor &Fb) const
The reduced matrix element.
Definition TensorOperator.cpp:62
DiracSpinor reduced_lhs(const int ka, const DiracSpinor &Fb) const
<b||h||a> = Fa * reduced_lhs(a, Fb) (a needed for angular factor)
Definition TensorOperator.cpp:55
void scale(double lambda)
Permanently re-scales the operator by constant, lambda.
Definition TensorOperator.cpp:84
double rme3js(const int twoja, const int twojb, int two_mb=1, int two_q=0) const
ME = rme3js * RME.
Definition TensorOperator.cpp:40
virtual DiracSpinor radial_rhs(const int kappa_a, const DiracSpinor &Fb) const
radial_int = Fa * radial_rhs(a, Fb) (a needed for angular factor)
Definition TensorOperator.cpp:88
virtual double radialIntegral(const DiracSpinor &Fa, const DiracSpinor &Fb) const
Defined via <a||h||b> = angularF(a,b) * radialIntegral(a,b) (Note: if radial_rhs is overridden,...
Definition TensorOperator.cpp:133
DiracSpinor reduced_rhs(const int ka, const DiracSpinor &Fb) const
<a||h||b> = Fa * reduced_rhs(a, Fb) (a needed for angular factor)
Definition TensorOperator.cpp:50
virtual void updateFrequency(const double)
Update frequency for frequency-dependant operators.
Definition TensorOperator.hpp:159
const std::vector< double > & getv() const
Returns a const ref to vector v.
Definition TensorOperator.hpp:165
Stores radial Dirac spinor: F_nk = (f, g)
Definition DiracSpinor.hpp:41
const std::vector< double > & f() const
Upper (large) radial component function, f.
Definition DiracSpinor.hpp:117
const std::vector< double > & g() const
Lower (small) radial component function, g.
Definition DiracSpinor.hpp:124
constexpr int neg1pow_2(int two_a)
Evaluates (-1)^{two_a/2} (for integer a; two_a is even)
Definition Wigner369j.hpp:152
constexpr int triangle(int j1, int j2, int J)
Returns 1 if triangle rule is satisfied. nb: works with j OR twoj!
Definition Wigner369j.hpp:162
Dirac Operators: General + derived.
Definition GenerateOperator.cpp:12
namespace qip::overloads provides operator overloads for std::vector
Definition Vector.hpp:450
std::vector< T > & operator+=(std::vector< T > &a, const std::vector< T > &b)
Provide addition of two vectors:
Definition Vector.hpp:454
4x4 Integer matrix (for Gamma/Pauli). Can be real or imag. Not mixed.
Definition TensorOperator.hpp:27