ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
Loading...
Searching...
No Matches
StructureRad.hpp
1#pragma once
2#include "Coulomb/include.hpp"
3#include "Coulomb/meTable.hpp"
4#include "IO/FRW_fileReadWrite.hpp"
5#include "Wavefunction/DiracSpinor.hpp"
6#include <memory>
7#include <vector>
8class Wavefunction;
9class DiracSpinor;
10namespace DiracOperator {
11class TensorOperator;
12}
13namespace ExternalField {
14class CorePolarisation;
15}
16
18namespace MBPT {
19
53
54public:
65 StructureRad(const std::vector<DiracSpinor> &basis, double en_core,
66 std::pair<int, int> nminmax = {0, 999},
67 const std::string &Qk_fname = "",
68 const std::vector<double> &fk = {},
69 const std::vector<double> &etak = {});
70
71private:
72 bool m_use_Qk;
73 std::vector<double> m_root_fk; // sqrt - add to both sides!
74 std::vector<double> m_etak;
76 std::optional<Coulomb::QkTable> mQ{std::nullopt}; //?
77 // nb: it seems conter-intuative, but this copy makes it FASTER!
78 std::vector<DiracSpinor> mCore{}, mExcited{};
79
80public:
81 const std::vector<DiracSpinor> &core() const { return mCore; }
82 const std::vector<DiracSpinor> &excited() const { return mExcited; }
83
85 const Coulomb::YkTable &Yk() const { return mY; }
86 const std::optional<Coulomb::QkTable> &Qk() const { return mQ; }
87
90 std::pair<double, double>
91 srTB(const DiracOperator::TensorOperator *const h, const DiracSpinor &w,
92 const DiracSpinor &v, double omega = 0.0,
93 const ExternalField::CorePolarisation *const dV = nullptr) const;
94
97 std::pair<double, double>
98 srC(const DiracOperator::TensorOperator *const h, const DiracSpinor &w,
99 const DiracSpinor &v,
100 const ExternalField::CorePolarisation *const dV = nullptr) const;
101
104 std::pair<double, double>
105 norm(const DiracOperator::TensorOperator *const h, const DiracSpinor &w,
106 const DiracSpinor &v,
107 const ExternalField::CorePolarisation *const dV = nullptr) const;
108
111 std::pair<double, double>
113 const DiracSpinor &v, double omega = 0.0,
114 const ExternalField::CorePolarisation *const dV = nullptr) const {
115 const auto [tb, dvtb] = srTB(h, w, v, omega, dV);
116 const auto [c, dvc] = srC(h, w, v, dV);
117 const auto [n, dvn] = norm(h, w, v, dV);
118 return {tb + c + n, dvtb + dvc + dvn};
119 }
120
122 double f_root_scr(int k) const {
123 const auto sk = std::size_t(k);
124 return sk < m_root_fk.size() ? m_root_fk[sk] : 1.0;
125 }
126
128 double eta_hp(int k) const {
129 const auto sk = std::size_t(k);
130 return sk < m_etak.size() ? m_etak[sk] : 1.0;
131 }
132
136 std::pair<double, double>
137 BO(const DiracOperator::TensorOperator *const h, const DiracSpinor &w,
138 const DiracSpinor &v,
139 const ExternalField::CorePolarisation *const dV = nullptr, double fw = 1.0,
140 double fv = 1.0) const;
141
143 double Sigma_vw(const DiracSpinor &v, const DiracSpinor &w) const;
144
148 const std::vector<DiracSpinor> &as,
149 const std::vector<DiracSpinor> &tbs = {}, double omega = 0.0,
150 const ExternalField::CorePolarisation *const dV = nullptr) const;
151
152private:
153 // "Top" diagrams
154 double t1234(int k, const DiracSpinor &w, const DiracSpinor &r,
155 const DiracSpinor &v, const DiracSpinor &c) const;
156 double t1(int k, const DiracSpinor &w, const DiracSpinor &r,
157 const DiracSpinor &v, const DiracSpinor &c) const;
158 double t2(int k, const DiracSpinor &w, const DiracSpinor &r,
159 const DiracSpinor &v, const DiracSpinor &c) const;
160 double t3(int k, const DiracSpinor &w, const DiracSpinor &r,
161 const DiracSpinor &v, const DiracSpinor &c) const;
162 double t4(int k, const DiracSpinor &w, const DiracSpinor &r,
163 const DiracSpinor &v, const DiracSpinor &c) const;
164
165 // "Bottom" diagrams
166 double b1234(int k, const DiracSpinor &w, const DiracSpinor &c,
167 const DiracSpinor &v, const DiracSpinor &r) const;
168
169 // "Centre" diagrams (most important)
170 double c1(int k, const DiracSpinor &w, const DiracSpinor &a,
171 const DiracSpinor &v, const DiracSpinor &c) const;
172 double c2(int k, const DiracSpinor &w, const DiracSpinor &a,
173 const DiracSpinor &v, const DiracSpinor &c) const;
174 double d1(int k, const DiracSpinor &w, const DiracSpinor &r,
175 const DiracSpinor &v, const DiracSpinor &m) const;
176 double d2(int k, const DiracSpinor &w, const DiracSpinor &r,
177 const DiracSpinor &v, const DiracSpinor &m) const;
178
179 // "Normalisation" terms (most important)
180 double n1(const DiracSpinor &v) const;
181 double n2(const DiracSpinor &v) const;
182 double dSigma_dE(const DiracSpinor &v, const DiracSpinor &i,
183 const DiracSpinor &j, const DiracSpinor &k) const;
184
185 std::pair<double, double>
186 z_bo(const DiracOperator::TensorOperator *const h, const DiracSpinor &w,
187 const DiracSpinor &v, bool transpose = false,
188 const ExternalField::CorePolarisation *const dV = nullptr) const;
189
190 //
191 double Q(int k, const DiracSpinor &a, const DiracSpinor &b,
192 const DiracSpinor &c, const DiracSpinor &d) const {
193 // inlcudes effective Coulomb screening
194 const auto fk = f_root_scr(k);
195 return m_use_Qk ? fk * mQ->Q(k, a, b, c, d) : fk * mY.Q(k, a, b, c, d);
196 }
197
198 double P(int k, const DiracSpinor &a, const DiracSpinor &b,
199 const DiracSpinor &c, const DiracSpinor &d) const {
200 // inlcudes effective Coulomb screening
201 return m_use_Qk ? mQ->P2(k, a, b, c, d, mY.SixJ(), m_root_fk) :
202 mY.P2(k, a, b, c, d, m_root_fk);
203 }
204
205 double W(int k, const DiracSpinor &a, const DiracSpinor &b,
206 const DiracSpinor &c, const DiracSpinor &d) const {
207 // inlcudes effective Coulomb screening
208 return Q(k, a, b, c, d) + P(k, a, b, c, d);
209 }
210};
211
212} // namespace MBPT
Calculates + stores Hartree Y functions + Angular (w/ look-up), taking advantage of symmetry.
Definition YkTable.hpp:26
double Q(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb, const DiracSpinor &Fc, const DiracSpinor &Fd) const
Calculates Qk using the existing yk integrals. Note: Yk and Ck tables must include all required value...
Definition YkTable.cpp:127
const Angular::SixJTable & SixJ() const
Returns a (const ref) to SixJ table [see Angular::SixJTable].
Definition YkTable.hpp:58
Look-up table for matrix elements. Note: does not assume any symmetry: (a,b) is stored independantly ...
Definition meTable.hpp:17
General operator (virtual base class); operators derive from this.
Definition TensorOperator.hpp:110
Stores radial Dirac spinor: F_nk = (f, g)
Definition DiracSpinor.hpp:41
Virtual Core Polarisation class, for <a||dV||b>. See TDHF, DiagramRPA, etc.
Definition CorePolarisation.hpp:35
Calculates Structure Radiation + Normalisation of states, using diagram method.
Definition StructureRad.hpp:52
std::pair< double, double > srTB(const DiracOperator::TensorOperator *const h, const DiracSpinor &w, const DiracSpinor &v, double omega=0.0, const ExternalField::CorePolarisation *const dV=nullptr) const
Returns sum of Top+Bottom (SR) diagrams, reduced ME: <w||T+B||v>. Returns a pair: {TB,...
Definition StructureRad.cpp:61
double f_root_scr(int k) const
Effective screening factor for Coulomb lines.
Definition StructureRad.hpp:122
std::pair< double, double > BO(const DiracOperator::TensorOperator *const h, const DiracSpinor &w, const DiracSpinor &v, const ExternalField::CorePolarisation *const dV=nullptr, double fw=1.0, double fv=1.0) const
Returns Brueckner orbital contribution to the, reduced ME: <w||h||v>_norm. Returns a pair: {N,...
Definition StructureRad.cpp:751
double eta_hp(int k) const
Effective hole-particle factor for polarisation loops.
Definition StructureRad.hpp:128
Coulomb::meTable< std::pair< double, double > > srn_table(const DiracOperator::TensorOperator *const h, const std::vector< DiracSpinor > &as, const std::vector< DiracSpinor > &tbs={}, double omega=0.0, const ExternalField::CorePolarisation *const dV=nullptr) const
constructs an me table of {srn, srn+dv} for each pair or {a,b}
Definition StructureRad.cpp:200
std::pair< double, double > srn(const DiracOperator::TensorOperator *const h, const DiracSpinor &w, const DiracSpinor &v, double omega=0.0, const ExternalField::CorePolarisation *const dV=nullptr) const
Returns sum of SR+Norm diagrams, reduced ME: <w||T+B+C+N||v>. Returns a pair: {SRN,...
Definition StructureRad.hpp:112
double Sigma_vw(const DiracSpinor &v, const DiracSpinor &w) const
Calculates <v|Sigma|w>
Definition StructureRad.cpp:768
const Coulomb::YkTable & Yk() const
Returns reference to Yk table. NOTE: may not be initialised!
Definition StructureRad.hpp:85
std::pair< double, double > norm(const DiracOperator::TensorOperator *const h, const DiracSpinor &w, const DiracSpinor &v, const ExternalField::CorePolarisation *const dV=nullptr) const
Returns Normalisation of states, reduced ME: <w||h||v>_norm. Returns a pair: {N, N+dV}: second includ...
Definition StructureRad.cpp:183
std::pair< double, double > srC(const DiracOperator::TensorOperator *const h, const DiracSpinor &w, const DiracSpinor &v, const ExternalField::CorePolarisation *const dV=nullptr) const
Returns Centre (SR) diagrams, reduced ME: <w||C||v>. Returns a pair: {C, C+dV}: second includes RPA (...
Definition StructureRad.cpp:111
Stores Wavefunction (set of valence orbitals, grid, HF etc.)
Definition Wavefunction.hpp:36
Dirac Operators: General + derived.
Definition GenerateOperator.cpp:12
Calculates many-body corrections (RPA) to matrix elements of external field.
Definition calcMatrixElements.cpp:14
Many-body perturbation theory.
Definition CI_Integrals.hpp:9