DiagramRPA.hpp

#pragma once
#include "CorePolarisation.hpp"
#include "Coulomb/QkTable.hpp"
#include "HF/Breit.hpp"
#include "IO/FRW_fileReadWrite.hpp"
#include "Wavefunction/DiracSpinor.hpp"
#include <vector>
class Wavefunction;
class DiracSpinor;
namespace DiracOperator {
class TensorOperator;
}
namespace HF {
class HartreeFock;
}
 
namespace ExternalField {
 
class DiagramRPA : public CorePolarisation {
 
private:
  const HF::HartreeFock *p_hf;
  std::vector<DiracSpinor> holes{};
  std::vector<DiracSpinor> excited{};
 
  std::optional<HF::Breit> m_Br{};
 
  // t0's never change
  // NO! They change if omega is updated (frequency dependent operator!)
  std::vector<std::vector<double>> t0am{};
  std::vector<std::vector<double>> t0ma{};
  // t's updated each solve_core itteration
  std::vector<std::vector<double>> tam{};
  std::vector<std::vector<double>> tma{};
 
  // Note: W's depend on rank (also parity)! Can re-use!?
  // These are probably an excellent candidate for unordered_map?
  std::vector<std::vector<std::vector<std::vector<double>>>> Wanmb{};
  std::vector<std::vector<std::vector<std::vector<double>>>> Wabmn{};
  std::vector<std::vector<std::vector<std::vector<double>>>> Wmnab{};
  std::vector<std::vector<std::vector<std::vector<double>>>> Wmban{};
 
public:
  DiagramRPA(const DiracOperator::TensorOperator *const h,
             const std::vector<DiracSpinor> &basis,
             const HF::HartreeFock *in_hf, const std::string &atom = "Atom");
 
  DiagramRPA(const DiracOperator::TensorOperator *const h,
             const DiagramRPA *const drpa);
 
public:
  virtual void solve_core(const double omega, int max_its = 200,
                          const bool print = true) override final;
 
  virtual Method method() const override final { return Method::diagram; }
 
  virtual double dV(const DiracSpinor &Fa,
                    const DiracSpinor &Fb) const override final;
 
  double dV_diagram(const DiracSpinor &Fa, const DiracSpinor &Fb) const;
 
  virtual void clear() override final;
 
  void update_t0s(const DiracOperator::TensorOperator *const h = nullptr);
 
  void grab_tam(const DiagramRPA *const drpa) {
    tam = drpa->tam;
    tma = drpa->tma;
  }
 
private:
  // Note: only writes W (depends on k/pi, and basis). Do not write t's, since
  // they depend on operator. This makes it very fast when making small changes
  // to operator (don't need to re-calc W)
  // Note: doesn't depend on grid!
  bool read_write(const std::string &fname, IO::FRW::RoW rw);
 
  void fill_W_matrix(const DiracOperator::TensorOperator *const h);
  void setup_ts(const DiracOperator::TensorOperator *const h);
 
public:
  DiagramRPA &operator=(const DiagramRPA &) = delete;
  DiagramRPA(const DiagramRPA &) = default;
  ~DiagramRPA() = default;
};
 
} // namespace ExternalField