ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
Loading...
Searching...
No Matches
hfs.hpp
1#pragma once
2#include "DiracOperator/TensorOperator.hpp"
3#include "IO/InputBlock.hpp"
4#include "Wavefunction/Wavefunction.hpp"
5#include "fmt/color.hpp"
6#include "qip/Maths.hpp"
7#include <functional>
8
9namespace DiracOperator {
10
11//==============================================================================
12
14namespace Hyperfine {
15using Func_R2_R = std::function<double(double, double)>; // save typing
16
18inline auto sphericalBall_F() -> Func_R2_R {
19 return [=](double r, double rN) {
20 return (r > rN) ? 1.0 : (r * r * r) / (rN * rN * rN);
21 };
22}
23
25inline auto sphericalShell_F() -> Func_R2_R {
26 return [=](double r, double rN) { return (r > rN) ? 1.0 : 0.0; };
27}
28
30inline auto pointlike_F() -> Func_R2_R {
31 return [=](double, double) { return 1.0; };
32}
33
34//------------------------------------------------------------------------------
36inline auto VolotkaSP_F(double mu, double I_nuc, double l_pn, int gl,
37 bool print = true) -> Func_R2_R
38// Function that returns generates + returns F_BW Bohr-Weiskopf
39// gl = 1 for proton, =0 for neutron. Double allowed for testing..
40// mu is in units of Nuclear Magneton!
41{
42 const auto two_I = int(2 * I_nuc + 0.0001);
43 const auto two_l = int(2 * l_pn + 0.0001);
44 const auto g_l = double(gl); // just safety
45 const auto gI = mu / I_nuc;
46
47 const auto K = (l_pn * (l_pn + 1.0) - (3. / 4.)) / (I_nuc * (I_nuc + 1.0));
48 const double g_s = (2.0 * gI - g_l * (K + 1.0)) / (1.0 - K);
49 if (print)
50 std::cout << "SingleParticle using: gl=" << g_l << ", gs=" << g_s
51 << ", l=" << l_pn << ", gI=" << gI << " (I=" << I_nuc << ")\n";
52 const double factor =
53 (two_I == two_l + 1) ?
54 g_s * (1 - two_I) / (4.0 * (two_I + 2)) + g_l * 0.5 * (two_I - 1) :
55 g_s * (3 + two_I) / (4.0 * (two_I + 2)) +
56 g_l * 0.5 * two_I * (two_I + 3) / (two_I + 2);
57 if (two_I != two_l + 1 && two_I != two_l - 1) {
58 std::cerr << "\nFAIL:59 in Hyperfine (VolotkaSP_F):\n "
59 "we must have I = l +/- 1/2, but we have: I,l = "
60 << I_nuc << "," << l_pn << "\n";
61 return [](double, double) { return 0.0; };
62 }
63 return [=](double r, double rN) {
64 return (r > rN) ? 1.0 :
65 ((r * r * r) / (rN * rN * rN)) *
66 (1.0 - (3.0 / mu) * std::log(r / rN) * factor);
67 };
68}
69
70//------------------------------------------------------------------------------
72inline auto uSP(double mu, double I_nuc, double l_pn, int gl, double n,
73 double R, bool u_option, bool print = true) -> Func_R2_R {
74 const auto two_I = int(2 * I_nuc + 0.0001);
75 const auto two_l = int(2 * l_pn + 0.0001);
76 const auto g_l = double(gl); // just safety
77 const auto gI = mu / I_nuc;
78
79 const auto K = (l_pn * (l_pn + 1.0) - (3. / 4.)) / (I_nuc * (I_nuc + 1.0));
80 const double g_s = (2.0 * gI - g_l * (K + 1.0)) / (1.0 - K);
81 if (print) {
82 std::cout << "uSP using: gl=" << g_l << ", gs=" << g_s << ", l=" << l_pn
83 << ", gI=" << gI << " (I=" << I_nuc << ")\n";
84 if (u_option) {
85 std::cout << "u(r) = u1(r) = u0 (R-r)^n ";
86 } else {
87 std::cout << "u(r) = u2(r) = u0 r^n ";
88 }
89 std::cout << "with n = " << n << "\n";
90 }
91 if (two_I != two_l + 1 && two_I != two_l - 1) {
92 std::cerr << "\nFAIL:59 in Hyperfine (uSP):\n "
93 "we must have I = l +/- 1/2, but we have: I,l = "
94 << I_nuc << "," << l_pn << "\n";
95 return [](double, double) { return 0.0; };
96 }
97
98 // nb: in theory, can read u(r) in from file!
99 // r^2 * u(r)^2
100 const auto r2u2 = [=](double r) {
101 return u_option ? r * r * std::pow(R - r, 2.0 * n) :
102 r * r * std::pow(r, 2.0 * n);
103 };
104
105 // normalisation
106 const auto u02 =
107 1.0 / NumCalc::num_integrate(r2u2, 0.0, R, 100, NumCalc::linear);
108
109 const auto F0r = [=](double r) {
110 const auto f = [=](double x) { return u02 * r2u2(x); };
111 return NumCalc::num_integrate(f, 0.0, r, 100, NumCalc::linear);
112 };
113
114 const auto FrR = [=](double r) {
115 const auto f = [=](double x) {
116 return u02 * r2u2(x) * r * r * r / x / x / x;
117 };
118 return NumCalc::num_integrate(f, r, R, 100, NumCalc::linear);
119 };
120
121 const auto two_Ip1 = 2.0 * (I_nuc + 1.0);
122 const auto twoI_p3 = 2.0 * I_nuc + 3.0;
123 const auto twoI_m1 = 2.0 * I_nuc - 1.0;
124
125 double f1 = (two_I == two_l + 1) ?
126 0.5 * g_s + (I_nuc - 0.5) * g_l :
127 -I_nuc / two_Ip1 * g_s + I_nuc * twoI_p3 / two_Ip1 * g_l;
128
129 double f2 =
130 (two_I == two_l + 1) ?
131 -twoI_m1 / (4.0 * two_Ip1) * g_s + (I_nuc - 0.5) * g_l :
132 twoI_p3 / (4.0 * two_Ip1) * g_s + I_nuc * twoI_p3 / two_Ip1 * g_l;
133
134 return [=](double r, double rN) {
135 return (r > rN) ? 1.0 : (1.0 / mu) * (f1 * F0r(r) + f2 * FrR(r));
136 };
137}
138
139//------------------------------------------------------------------------------
141inline auto doublyOddSP_F(double mut, double It, double mu1, double I1,
142 double l1, int gl1, double I2, double l2,
143 bool print = true) -> Func_R2_R
144// F(r) * g = 0.5 [ g1F1 + g2F2 + (g1F1 - g2F2) * K]
145// K = [I1(I1+1) - I2(I2+1)] / [I(I+1)]
146// return F(r) [divide by g]
147// VolotkaSP_F(mu, I, l, g_l); //gl is 1 or 0
148// g2 : from: g = 0.5 [ g1 + g2 + (g1 - g2) * K]
149{
150 const auto K = (I1 * (I1 + 1.0) - I2 * (I2 + 1.0)) / (It * (It + 1.0));
151 const auto gt = mut / It;
152 const auto g1 = mu1 / I1;
153 const auto g2 = (g1 * (K + 1.0) - 2.0 * gt) / (K - 1.0);
154 const auto mu2 = g2 * I2;
155 const auto gl2 = (gl1 == 0) ? 1 : 0;
156 const auto F1 = VolotkaSP_F(mu1, I1, l1, gl1, print);
157 const auto F2 = VolotkaSP_F(mu2, I2, l2, gl2, print);
158 return [=](double r, double rN) {
159 return (0.5 / gt) * (g1 * F1(r, rN) + g2 * F2(r, rN) +
160 K * (g1 * F1(r, rN) - g2 * F2(r, rN)));
161 };
162}
163
164//------------------------------------------------------------------------------
166inline double convert_RME_to_AB_2J(int k, int tja, int tjb) {
167 const auto tjz = std::min(tja, tjb); // arbitrary choice
168 const auto s = Angular::neg1pow_2(tja - tjb);
169 const auto tjs = Angular::threej_2(tja, 2 * k, tjb, -tjz, 0, tjz);
170 const auto f = k % 2 == 0 ? 2.0 : 1 / (0.5 * tjz);
171 return s * f * tjs;
172}
173
175inline double convert_RME_to_AB(int k, int ka, int kb) {
176 const auto tja = Angular::twoj_k(ka);
177 const auto tjb = Angular::twoj_k(kb);
178 return convert_RME_to_AB_2J(k, tja, tjb);
179}
180
181inline double hfsA(const TensorOperator *h, const DiracSpinor &Fa) {
182 auto Raa = h->radialIntegral(Fa, Fa); //?
183 return Raa * Fa.kappa() / (Fa.jjp1()) * PhysConst::muN_CGS_MHz;
184}
185
186inline double hfsB(const TensorOperator *h, const DiracSpinor &Fa) {
187 auto Raa = h->radialIntegral(Fa, Fa); //?
188 return Raa * double(Fa.twoj() - 1) / double(Fa.twoj() + 2) *
189 PhysConst::barn_MHz;
190}
191
192//==============================================================================
194inline std::vector<double> RadialFunc(int k, double rN, const Grid &rgrid,
195 const Func_R2_R &hfs_F) {
196 std::vector<double> rfunc;
197 rfunc.reserve(rgrid.num_points());
198 for (const auto r : rgrid)
199 rfunc.push_back(hfs_F(r, rN) / pow(r, k + 1));
200 return rfunc;
201}
202
203} // namespace Hyperfine
204
205//==============================================================================
208class hfs final : public TensorOperator {
209 // see Xiao, ..., Derevianko, Phys. Rev. A 102, 022810 (2020).
210 using Func_R2_R = std::function<double(double, double)>;
211
212public:
213 hfs(int in_k, double in_GQ, double rN_au, const Grid &rgrid,
214 const Func_R2_R &hfs_F = Hyperfine::pointlike_F(), bool MHzQ = true)
215 : TensorOperator(in_k, Parity::even, in_GQ,
216 Hyperfine::RadialFunc(in_k, rN_au, rgrid, hfs_F), 0),
217 k(in_k),
218 magnetic(k % 2 != 0),
219 cfg(magnetic ? 1.0 : 0.0),
220 cff(magnetic ? 0.0 : 1.0),
221 mMHzQ(MHzQ) {}
222
223 std::string name() const override final {
224 return "hfs" + std::to_string(k) + "";
225 }
226 std::string units() const override final {
227 return (k <= 2 && mMHzQ) ? "MHz" : "au";
228 }
229
230 double angularF(const int ka, const int kb) const override final {
231 // inludes unit: Assumes g in nuc. magneton units, and/or Q in barns
232 // This only for k=1 (mag dipole) and k=2 E. quad.
233 const auto unit = (mMHzQ && k == 1) ? PhysConst::muN_CGS_MHz :
234 (mMHzQ && k == 2) ? PhysConst::barn_MHz :
235 1.0;
236 return magnetic ? -(ka + kb) * Angular::Ck_kk(k, -ka, kb) * unit :
237 -Angular::Ck_kk(k, ka, kb) * unit;
238 }
239
240 double angularCff(int, int) const override final { return cff; }
241 double angularCgg(int, int) const override final { return cff; }
242 double angularCfg(int, int) const override final { return cfg; }
243 double angularCgf(int, int) const override final { return cfg; }
244
245private:
246 int k;
247 bool magnetic;
248 double cfg;
249 double cff;
250 bool mMHzQ;
251};
252
253//==============================================================================
254//==============================================================================
255inline std::unique_ptr<DiracOperator::TensorOperator>
256generate_hfs(const IO::InputBlock &input, const Wavefunction &wf) {
257 using namespace DiracOperator;
258
259 input.check(
260 {{"", "Most following will be taken from the default nucleus if "
261 "not explicitely given"},
262 {"mu", "Magnetic moment in mu_N"},
263 {"Q", "Nuclear quadrupole moment, in barns. Also used as overall "
264 "constant for any higher-order moments [1.0]"},
265 {"k", "Multipolarity. 1=mag. dipole, 2=elec. quad, etc. [1]"},
266 {"rrms",
267 "nuclear (magnetic) rms radius, in Fermi (fm) (defult is charge rms)"},
268 {"units", "Units for output (only for k=1,k=2). MHz or au [MHz]"},
269 {"F", "F(r): Nuclear moment distribution: ball, point, shell, "
270 "SingleParticle, or doublyOddSP [ball]"},
271 {"F(r)", "Obselete; use 'F' from now - will be removed"},
272 {"nuc_mag", "Obselete; use 'F' from now - will be removed"},
273 {"printF", "Writes F(r) to a text file [false]"},
274 {"print", "Write F(r) info to screen [true]"},
275 {"", "The following are only for F=SingleParticle or doublyOddSP"},
276 {"I", "Nuclear spin. Taken from nucleus"},
277 {"parity", "Nulcear parity: +/-1"},
278 {"l", "l for unpaired nucleon (automatically derived from I and "
279 "parity; best to leave as default)"},
280 {"gl", "=1 for proton, =0 for neutron"},
281 {"", "The following are only used if F=doublyOddSP"},
282 {"mu1", "mag moment of 'first' unpaired nucleon"},
283 {"gl1", "gl of 'first' unpaired nucleon"},
284 {"l1", "l of 'first' unpaired nucleon"},
285 {"l2", "l of 'second' unpaired nucleon"},
286 {"I1", "total spin (J) of 'first' unpaired nucleon"},
287 {"I2", "total spin (J) of 'second' unpaired nucleon"},
288 {"", "The following are only for u(r) function"},
289 {"u", "u1 or u2 : u1(r) = (R-r)^n, u2(r) = r^n [u1]"},
290 {"n", "n that appears above. Should be between 0 and 2 [0]"}});
291 if (input.has_option("help")) {
292 return nullptr;
293 }
294
295 const auto nuc = wf.nucleus();
296 const auto isotope = Nuclear::findIsotopeData(nuc.z(), nuc.a());
297 const auto mu = input.get("mu", isotope.mu ? *isotope.mu : 0.0);
298 const auto I_nuc = input.get("I", isotope.I_N ? *isotope.I_N : 0.0);
299 const auto print = input.get("print", true);
300 const auto k = input.get("k", 1);
301
302 const auto use_MHz =
303 (k <= 2 &&
304 qip::ci_compare(input.get<std::string>("units", "MHz"), "MHz"));
305
306 if (k <= 0) {
307 fmt2::styled_print(fg(fmt::color::red), "\nError 246:\n");
308 std::cout << "In hyperfine: invalid K=" << k << "! meaningless results\n";
309 }
310 if (I_nuc <= 0) {
311 fmt2::styled_print(fg(fmt::color::orange), "\nWarning 253:\n");
312 std::cout << "In hyperfine: invalid I_nuc=" << I_nuc
313 << "! meaningless results\n";
314 }
315
316 const auto g_or_Q = (k == 1) ? (mu / I_nuc) :
317 (k == 2) ? input.get("Q", isotope.q ? *isotope.q : 1.0) :
318 input.get("Q", 1.0);
319
320 enum class DistroType {
321 point,
322 ball,
323 shell,
324 SingleParticle,
325 doublyOddSP,
326 spu,
327 Error
328 };
329
330 std::string default_distribution = "ball";
331
332 // For compatability with old notation of 'F(r)' input option
333 const auto Fr_str =
334 input.has_option("F") ?
335 input.get<std::string>("F", default_distribution) :
336 input.has_option("nuc_mag") ?
337 input.get<std::string>("nuc_mag", default_distribution) :
338 input.get<std::string>("F(r)", default_distribution);
339
340 const auto distro_type =
341 (qip::ci_wc_compare(Fr_str, "point*") || qip::ci_compare(Fr_str, "1")) ?
342 DistroType::point :
343 qip::ci_compare(Fr_str, "ball") ? DistroType::ball :
344 qip::ci_compare(Fr_str, "shell") ? DistroType::shell :
345 qip::ci_wc_compare(Fr_str, "Single*") ? DistroType::SingleParticle :
346 qip::ci_compare(Fr_str, "doublyOdd*") ? DistroType::doublyOddSP :
347 qip::ci_compare(Fr_str, "spu") ? DistroType::spu :
348 DistroType::Error;
349 if (distro_type == DistroType::Error) {
350 fmt2::styled_print(fg(fmt::color::red), "\nError 271:\n");
351 std::cout << "\nIn hyperfine. Unkown F(r) - " << Fr_str << "\n";
352 std::cout << "Defaulting to pointlike!\n";
353 }
354
355 const auto r_rmsfm =
356 distro_type == DistroType::point ? 0.0 : input.get("rrms", nuc.r_rms());
357 const auto r_nucfm = std::sqrt(5.0 / 3) * r_rmsfm;
358 const auto r_nucau = r_nucfm / PhysConst::aB_fm;
359
360 if (print) {
361 std::cout << "\nHyperfine structure: " << wf.atom() << "\n";
362 std::cout << "K=" << k << " ("
363 << (k == 1 ? "magnetic dipole" :
364 k == 2 ? "electric quadrupole" :
365 k % 2 == 0 ? "electric multipole" :
366 "magnetic multipole")
367 << ")\n";
368 std::cout << "Using " << Fr_str << " nuclear distro for F(r)\n"
369 << "w/ r_N = " << r_nucfm << "fm = " << r_nucau
370 << "au (r_rms=" << r_rmsfm << "fm)\n";
371 std::cout << "Points inside nucleus: " << wf.grid().getIndex(r_nucau)
372 << "\n";
373 if (k == 1) {
374 std::cout << "mu = " << mu << ", I = " << I_nuc << ", g = " << g_or_Q
375 << "\n";
376 } else {
377 std::cout << "Q = " << g_or_Q << "\n";
378 }
379 }
380
381 // default is BALL:
382 auto Fr = Hyperfine::sphericalBall_F();
383 if (distro_type == DistroType::ball) {
385 } else if (distro_type == DistroType::shell) {
387 } else if (distro_type == DistroType::SingleParticle) {
388 const auto pi = input.get("parity", isotope.parity ? *isotope.parity : 0);
389 const auto l_tmp = int(I_nuc + 0.5 + 0.0001);
390 auto l = ((l_tmp % 2 == 0) == (pi == 1)) ? l_tmp : l_tmp - 1;
391 l = input.get("l", l); // can override derived 'l' (not recommended)
392 const auto gl_default = wf.Znuc() % 2 == 0 ? 0 : 1; // unparied proton?
393 const auto gl = input.get<int>("gl", gl_default);
394 if (print) {
395 std::cout << "Single-Particle (Volotka formula) for unpaired";
396 if (gl == 1)
397 std::cout << " proton ";
398 else if (gl == 0)
399 std::cout << " neturon ";
400 else
401 std::cout << " gl=" << gl << "??? program will run, but prob wrong!\n";
402 std::cout << "with l=" << l << " (pi=" << pi << ")\n";
403 }
404 Fr = Hyperfine::VolotkaSP_F(mu, I_nuc, l, gl, print);
405 } else if (distro_type == DistroType::spu) {
406 const auto pi = input.get("parity", isotope.parity ? *isotope.parity : 0);
407 const auto u_func = input.get("u", std::string{"u1"}); // u1=(R-r)^n, u2=r^n
408 const bool u_option = u_func == std::string{"u1"};
409 const auto n = input.get("n", 0.0); // u1=(R-r)^n, u2=r^n
410 const auto l_tmp = int(I_nuc + 0.5 + 0.0001);
411 auto l = ((l_tmp % 2 == 0) == (pi == 1)) ? l_tmp : l_tmp - 1;
412 l = input.get("l", l); // can override derived 'l' (not recommended)
413 const auto gl_default = wf.Znuc() % 2 == 0 ? 0 : 1; // unparied proton?
414 const auto gl = input.get<int>("gl", gl_default);
415 if (print) {
416 std::cout << "Single-Particle (Volotka formula) with u(r) for unpaired";
417 if (gl == 1)
418 std::cout << " proton ";
419 else if (gl == 0)
420 std::cout << " neturon ";
421 else
422 std::cout << " gl=" << gl << "??? program will run, but prob wrong!\n";
423 std::cout << "with l=" << l << " (pi=" << pi << ")\n";
424 }
425 Fr = Hyperfine::uSP(mu, I_nuc, l, gl, n, r_nucau, u_option, print);
426 } else if (distro_type == DistroType::doublyOddSP) {
427 const auto mu1 = input.get<double>("mu1", 1.0);
428 const auto gl1 = input.get<int>("gl1", -1); // 1 or 0 (p or n)
429 if (gl1 != 0 && gl1 != 1) {
430 fmt2::styled_print(fg(fmt::color::red), "\nError 324:\n");
431 std::cout << "In " << input.name() << " " << Fr_str
432 << "; have gl1=" << gl1 << " but need 1 or 0\n";
433 return std::make_unique<NullOperator>(NullOperator());
434 }
435 const auto l1 = input.get<double>("l1", -1.0);
436 const auto l2 = input.get<double>("l2", -1.0);
437 const auto I1 = input.get<double>("I1", -1.0);
438 const auto I2 = input.get<double>("I2", -1.0);
439
440 Fr = Hyperfine::doublyOddSP_F(mu, I_nuc, mu1, I1, l1, gl1, I2, l2, print);
441 }
442
443 // Optionally print F(r) function to file
444 if (input.get<bool>("printF", false)) {
445 std::ofstream of(wf.identity() + "_" + Fr_str + ".txt");
446 of << "r/fm F(r)\n";
447 for (auto r : wf.grid()) {
448 of << r * PhysConst::aB_fm << " "
449 << Fr(r * PhysConst::aB_fm, r_nucau * PhysConst::aB_fm) << "\n";
450 }
451 }
452
453 return std::make_unique<hfs>(k, g_or_Q, r_nucau, wf.grid(), Fr, use_MHz);
454}
455
456} // namespace DiracOperator
Speacial operator: 0.
Definition TensorOperator.hpp:279
General operator (virtual base class); operators derive from this.
Definition TensorOperator.hpp:110
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
Units: Assumes g in nuc. magneton units (magnetic), and Q in barns (electric), and rN is atomic units...
Definition hfs.hpp:208
std::string units() const override final
Returns units of operator (usually au, may be MHz, etc.)
Definition hfs.hpp:226
std::string name() const override final
Returns "name" of operator (e.g., 'E1')
Definition hfs.hpp:223
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 hfs.hpp:230
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
int kappa() const
Dirac quantum number, kappa.
Definition DiracSpinor.hpp:84
double jjp1() const
j(j+1)
Definition DiracSpinor.hpp:88
Holds grid, including type + Jacobian (dr/du)
Definition Grid.hpp:31
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
Holds list of Options, and a list of other InputBlocks. Can be initialised with a list of options,...
Definition InputBlock.hpp:142
bool check(std::initializer_list< std::string > blocks, const std::vector< std::pair< std::string, std::string > > &list, bool print=false) const
Check all the options and blocks in this; if any of them are not present in 'list',...
Definition InputBlock.hpp:594
bool has_option(std::string_view key) const
Check is option is present (even if not set) in current block.
Definition InputBlock.hpp:201
T get(std::string_view key, T default_value) const
If 'key' exists in the options, returns value. Else, returns default_value. Note: If two keys with sa...
Definition InputBlock.hpp:417
Stores Wavefunction (set of valence orbitals, grid, HF etc.)
Definition Wavefunction.hpp:36
const Grid & grid() const
Returns a const reference to the radial grid.
Definition Wavefunction.hpp:81
int Znuc() const
Nuclear charge, Z.
Definition Wavefunction.hpp:98
std::string identity() const
Atomic symbol, including core ionisation degree and run_label.
Definition Wavefunction.cpp:687
std::string atom() const
String of atom info (e.g., "Cs, Z=55, A=133")
Definition Wavefunction.hpp:189
const Nuclear::Nucleus & nucleus() const
Returns Nuclear::nucleus object (contains nuc. parameters)
Definition Wavefunction.hpp:96
double threej_2(int two_j1, int two_j2, int two_j3, int two_m1, int two_m2, int two_m3)
Calculates wigner 3j symbol. Takes in 2*j (or 2*l) - intput is integer.
Definition Wigner369j.hpp:177
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 twoj_k(int ka)
returns 2j given kappa
Definition Wigner369j.hpp:45
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
auto doublyOddSP_F(double mut, double It, double mu1, double I1, double l1, int gl1, double I2, double l2, bool print=true) -> Func_R2_R
'Volotka' SP model, for doubly-odd nuclei: Phys. Rev. Lett. 125, 063002 (2020)
Definition hfs.hpp:141
auto uSP(double mu, double I_nuc, double l_pn, int gl, double n, double R, bool u_option, bool print=true) -> Func_R2_R
Elizarov et al., Optics and Spectroscopy (2006): u(r) = u0(R-r)^n.
Definition hfs.hpp:72
auto pointlike_F() -> Func_R2_R
Pointlike F(r): 1.
Definition hfs.hpp:30
auto sphericalBall_F() -> Func_R2_R
Spherical ball F(r): (r/rN)^3 for r<rN, 1 for r>rN.
Definition hfs.hpp:18
double convert_RME_to_AB(int k, int ka, int kb)
Converts reduced matrix element to A/B coeficients (takes k, kappa, kappa)
Definition hfs.hpp:175
double convert_RME_to_AB_2J(int k, int tja, int tjb)
Converts reduced matrix element to A/B coeficients (takes k, 2J, 2J)
Definition hfs.hpp:166
auto VolotkaSP_F(double mu, double I_nuc, double l_pn, int gl, bool print=true) -> Func_R2_R
'Volotka' single-particle model: see Phys. Rev. Lett. 125, 063002 (2020)
Definition hfs.hpp:36
auto sphericalShell_F() -> Func_R2_R
Spherical shell F(r): 0 for r<rN, 1 for r>rN.
Definition hfs.hpp:25
std::vector< double > RadialFunc(int k, double rN, const Grid &rgrid, const Func_R2_R &hfs_F)
Takes in F(r) and k, and forms hyperfine radial function: F(r,rN)/r^{k+1}.
Definition hfs.hpp:194
Dirac Operators: General + derived.
Definition GenerateOperator.cpp:12
Isotope findIsotopeData(int z, int a)
Looks up + returns an isotope from the list. If not in list, partially blank.
Definition NuclearData.cpp:6
constexpr double muN_CGS_MHz
Nulcear magneton in MHz (via Gaussian CGS-derived atomic units):
Definition PhysConst_constants.hpp:103
qip library: A collection of useful functions
Definition Array.hpp:9
bool ci_wc_compare(std::string_view s1, std::string_view s2)
Compares two strings, s1 and s2. s2 may contain ONE wildcard ('*') which will match anything....
Definition String.hpp:139
bool ci_compare(std::string_view s1, std::string_view s2)
Case insensitive string compare. Essentially: LowerCase(s1)==LowerCase(s2)
Definition String.hpp:131