YkTable.hpp

#pragma once
#include "Angular/CkTable.hpp"
#include "Angular/SixJTable.hpp"
#include "Wavefunction/DiracSpinor.hpp"
#include <unordered_map>
#include <vector>
 
namespace Coulomb {
 
 
class YkTable {
 
private:
  std::vector<std::unordered_map<uint32_t, std::vector<double>>> m_Y{};
  Angular::CkTable m_Ck{};
  Angular::SixJTable m_6j{};
 
public:
  YkTable(const std::vector<DiracSpinor> &a_orbs,
          const std::vector<DiracSpinor> &b_orbs) {
    calculate(a_orbs, b_orbs);
  }
  YkTable(const std::vector<DiracSpinor> &a_orbs) { calculate(a_orbs); }
  YkTable() {}
 
  //----------------------------------------------------------------------------
  void calculate(const std::vector<DiracSpinor> &as,
                 const std::vector<DiracSpinor> &bs);
  void calculate(const std::vector<DiracSpinor> &as) { calculate(as, as); }
 
  void extend_angular(int new_max_2j);
 
  const Angular::CkTable &Ck() const { return m_Ck; }
  const Angular::SixJTable &SixJ() const { return m_6j; }
 
  const std::vector<double> *get(const int k, const DiracSpinor &Fa,
                                 const DiracSpinor &Fb) const;
 
  [[nodiscard]] double R(const int k, const DiracSpinor &Fa,
                         const DiracSpinor &Fb, const DiracSpinor &Fc,
                         const DiracSpinor &Fd) const;
 
  [[nodiscard]] double Q(const int k, const DiracSpinor &Fa,
                         const DiracSpinor &Fb, const DiracSpinor &Fc,
                         const DiracSpinor &Fd) const;
 
  [[nodiscard]] double P(const int k, const DiracSpinor &Fa,
                         const DiracSpinor &Fb, const DiracSpinor &Fc,
                         const DiracSpinor &Fd) const;
 
  [[nodiscard]] double P2(const int k, const DiracSpinor &Fa,
                          const DiracSpinor &Fb, const DiracSpinor &Fc,
                          const DiracSpinor &Fd,
                          const std::vector<double> &fk = {}) const;
 
  [[nodiscard]] double W(const int k, const DiracSpinor &Fa,
                         const DiracSpinor &Fb, const DiracSpinor &Fc,
                         const DiracSpinor &Fd) const;
 
  [[nodiscard]] DiracSpinor Qkv_bcd(const int k, int kappa,
                                    const DiracSpinor &Fb,
                                    const DiracSpinor &Fc,
                                    const DiracSpinor &Fd) const;
 
  [[nodiscard]] DiracSpinor
  Pkv_bcd(const int k, int kappa, const DiracSpinor &Fb, const DiracSpinor &Fc,
          const DiracSpinor &Fd, const std::vector<double> &f2k = {}) const;
 
private:
  // Allocates space for the Yk table, but does not calculate Yk. This is
  // because allocation cannot be done in parallel, but once allocation is done,
  // calculation can be done in //
  void allocate_space(const std::vector<DiracSpinor> &a_orbs,
                      const std::vector<DiracSpinor> &b_orbs);
 
  // Returns key used for look-up table (unordered_map), considers symmetry
  uint32_t ab_key(const DiracSpinor &Fa, const DiracSpinor &Fb) const;
 
  // Returns reference to yk vector at 'key'; if no such vector exists, creates
  // it first.
  std::vector<double> &get_or_insert(std::size_t k, uint32_t key);
 
  // Returns reference to yk vector at 'key'. This vector *must* exist already
  std::vector<double> &get_ref(const int k, const DiracSpinor &Fa,
                               const DiracSpinor &Fb);
};
 
} // namespace Coulomb