ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
Loading...
Searching...
No Matches
DiagramRPA0_jL.hpp
1#pragma once
2#include "CorePolarisation.hpp"
3#include "Coulomb/YkTable.hpp"
4#include "DiracOperator/include.hpp"
5#include "HF/HartreeFock.hpp"
6#include "IO/FRW_fileReadWrite.hpp"
7#include "Wavefunction/DiracSpinor.hpp"
8#include <vector>
9
10namespace ExternalField {
11
14
15private:
16 const HF::HartreeFock *p_hf;
17 std::vector<DiracSpinor> holes{};
18 std::vector<DiracSpinor> excited{};
19 Coulomb::YkTable m_yk{};
20
21public:
23 const std::vector<DiracSpinor> &basis,
24 const HF::HartreeFock *in_hf, int max_l)
25 : CorePolarisation(h), p_hf(in_hf) {
26
27 if (p_hf == nullptr || h == nullptr) {
28 std::cout << "\nFAIL:25 in DiagramRPA0_jL - hf cannot be null\n"
29 << std::flush;
30 }
31
32 // Set up basis:
33 const auto &core = p_hf->core();
34 for (const auto &Fi : basis) {
35 const bool inCore = std::find(cbegin(core), cend(core), Fi) != cend(core);
36 if (inCore) {
37 holes.push_back(Fi);
38 } else {
39 excited.push_back(Fi);
40 }
41 }
42
43 m_yk.calculate(basis);
44 // need Ck up to continuum states; depends on max_l
45 const auto max_tj = DiracSpinor::max_tj(basis) + 2 * max_l;
46 m_yk.extend_angular(max_tj); //??
47 }
48
49public:
51 virtual void solve_core(const double, int, const bool) override final {
52 std::cout << "ERROR: solve_core() not implemented\n";
53 }
54
56 virtual Method method() const override final { return Method::Error; }
57
59 virtual double dV(const DiracSpinor &,
60 const DiracSpinor &) const override final {
61 std::cout << "\nWARNING: not implemented. Use dV_diagram_jL()\n"
62 << std::flush;
63 std::abort();
64 return 0.0;
65 }
66
67 //----------------------------------------------------------------------------
68 // nb: Fv MUST appear within basis!
69 double dV_diagram_jL(const DiracSpinor &Fw, const DiracSpinor &Fv,
70 const DiracOperator::jL *jl, std::size_t L,
71 double q) const {
72
73 if (holes.empty() || excited.empty())
74 return 0.0;
75
76 if (Fv.en() > Fw.en()) {
77 const auto sj = ((Fv.twoj() - Fw.twoj()) % 4 == 0) ? 1 : -1;
78 const auto si = m_imag ? -1 : 1;
79 return (sj * si) * dV_diagram_jL(Fv, Fw, jl, L, q);
80 }
81
82 const auto ww = 0.0;
83 const auto iL = int(L);
84
85 const auto f = (1.0 / (2 * iL + 1));
86
87 std::vector<double> sum_a(holes.size());
88#pragma omp parallel for
89 for (std::size_t ia = 0; ia < holes.size(); ia++) {
90 const auto &Fa = holes[ia];
91 const auto s1 = ((Fa.twoj() - Fw.twoj() + 2 * iL) % 4 == 0) ? 1 : -1;
92 for (std::size_t im = 0; im < excited.size(); im++) {
93 const auto &Fm = excited[im];
94 // if (t0am[ia][im] == 0.0)
95 if (jl->is_zero(Fa, Fm, L))
96 continue;
97 const auto s2 = ((Fa.twoj() - Fm.twoj()) % 4 == 0) ? 1 : -1;
98 // const auto Wwmva = Coulomb::Wk_abcd(Fw, Fm, Fv, Fa, L);
99 // const auto Wwavm = Coulomb::Wk_abcd(Fw, Fa, Fv, Fm, L);
100 const auto Wwmva = m_yk.W(iL, Fw, Fm, Fv, Fa);
101 const auto Wwavm = m_yk.W(iL, Fw, Fa, Fv, Fm);
102 auto ttam = jl->rme(Fa, Fm, L, q);
103
104 const auto ttma = jl->symm_sign(Fa, Fm) * ttam;
105 const auto A = ttam * Wwmva / (Fa.en() - Fm.en() - ww);
106 const auto B = Wwavm * ttma / (Fa.en() - Fm.en() + ww);
107 sum_a[ia] += s1 * (A + s2 * B);
108 }
109 }
110 return f * std::accumulate(begin(sum_a), end(sum_a), 0.0);
111 }
112
113 void clear() override final { return; }
114
115public:
116 DiagramRPA0_jL &operator=(const DiagramRPA0_jL &) = delete;
117 DiagramRPA0_jL(const DiagramRPA0_jL &) = default;
118 ~DiagramRPA0_jL() = default;
119};
120
121} // namespace ExternalField
Calculates + stores Hartree Y functions + Angular (w/ look-up), taking advantage of symmetry.
Definition YkTable.hpp:26
void calculate(const std::vector< DiracSpinor > &as, const std::vector< DiracSpinor > &bs)
Re-calculates all y_ab functions (will over-ride existing ones); NOTE: only calculates for a in as,...
Definition YkTable.cpp:14
double W(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb, const DiracSpinor &Fc, const DiracSpinor &Fd) const
Calculates Wk=Qk+Pk using the existing yk integrals. Note: Yk and Ck tables must include all required...
Definition YkTable.cpp:195
void extend_angular(int new_max_2j)
Extends the Ck and 6J tables up to new maximum 2*j.
Definition YkTable.cpp:46
General operator (virtual base class); operators derive from this.
Definition TensorOperator.hpp:110
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
Matrix element of tensor operator: J_L(qr)*C^L.
Definition jL.hpp:19
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
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
Stores radial Dirac spinor: F_nk = (f, g)
Definition DiracSpinor.hpp:41
static int max_tj(const std::vector< DiracSpinor > &orbs)
Returns maximum (2j) found in {orbs}.
Definition DiracSpinor.cpp:266
double en() const
Single-particle energy, not including rest energy.
Definition DiracSpinor.hpp:113
Virtual Core Polarisation class, for <a||dV||b>. See TDHF, DiagramRPA, etc.
Definition CorePolarisation.hpp:35
RPA correction to matrix elements, using Diagram technique.
Definition DiagramRPA0_jL.hpp:13
virtual void solve_core(const double, int, const bool) override final
Itterates the RPA equations for core electrons.
Definition DiagramRPA0_jL.hpp:51
void clear() override final
Clears the dPsi orbitals (sets to zero)
Definition DiagramRPA0_jL.hpp:113
virtual Method method() const override final
Returns RPA method.
Definition DiagramRPA0_jL.hpp:56
virtual double dV(const DiracSpinor &, const DiracSpinor &) const override final
Calculates RPA correction to matrix element: <A||dV||B>
Definition DiagramRPA0_jL.hpp:59
Solves relativistic Hartree-Fock equations for core and valence. Optionally includes Breit and QED ef...
Definition HartreeFock.hpp:70
const std::vector< DiracSpinor > & core() const
vector of core orbitals
Definition HartreeFock.hpp:186
double f(RaB r, PrincipalQN n, DiracQN k, Zeff z, AlphaFS a)
Upper radial component.
Definition DiracHydrogen.cpp:71
Calculates many-body corrections (RPA) to matrix elements of external field.
Definition calcMatrixElements.cpp:14