QkTable.hpp

#pragma once
#include "Angular/SixJTable.hpp"
#include "Wavefunction/DiracSpinor.hpp"
#include "YkTable.hpp"
#include <array>
#include <cstdint>
#include <functional>
#include <unordered_map>
 
namespace Coulomb {
 
enum class Symmetry { Qk, Wk, Lk, none };
using Real = double;
using nk4Index = uint64_t;
using nkIndex = uint16_t;
 
static_assert(sizeof(nkIndex) == sizeof(DiracSpinor::Index));
 
using CoulombFunction =
    std::function<double(int, const DiracSpinor &a, const DiracSpinor &b,
                         const DiracSpinor &c, const DiracSpinor &d)>;
 
using SelectionRules =
    std::function<bool(int, const DiracSpinor &a, const DiracSpinor &b,
                       const DiracSpinor &c, const DiracSpinor &d)>;
 
//==============================================================================
template <Symmetry S>
class CoulombTable {
 
private:
  // each vector element corresponds to a 'k'
  std::vector<std::unordered_map<nk4Index, Real>> m_data{};
 
public:
  // 'Rule of zero' (except virtual destructor)
  // virtual ~CoulombTable() = default;
 
  void fill(const std::vector<DiracSpinor> &basis, const YkTable &yk,
            int k_cut = -1, bool print = true);
 
  void fill_if(const std::vector<DiracSpinor> &basis, const YkTable &yk,
               const SelectionRules &SelectionFunction, int k_cut = -1,
               bool print = true);
 
  void fill(const std::vector<DiracSpinor> &basis, const CoulombFunction &Fk,
            const SelectionRules &Fk_SR, int k_cut = -1, bool print = true);
 
  auto operator->() { return &m_data; }
 
  void summary() const;
  std::size_t count() const;
 
  void add(int k, const DiracSpinor &a, const DiracSpinor &b,
           const DiracSpinor &c, const DiracSpinor &d, Real value);
  void add(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d, Real value);
  void add(int k, nk4Index index, Real value);
 
  void update(int k, const DiracSpinor &a, const DiracSpinor &b,
              const DiracSpinor &c, const DiracSpinor &d, Real value);
  void update(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d, Real value);
  void update(int k, nk4Index index, Real value);
 
  bool contains(int k, const DiracSpinor &a, const DiracSpinor &b,
                const DiracSpinor &c, const DiracSpinor &d) const;
  bool contains(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d) const;
  bool contains(int k, nk4Index index) const;
 
  bool emptyQ() const {
    return (m_data.size() == 0);
    // nb: possible to be empty if size > 0... (but doesn't matter..)
  }
 
  Real Q(int k, const DiracSpinor &a, const DiracSpinor &b,
         const DiracSpinor &c, const DiracSpinor &d) const;
 
  Real Q(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d) const;
  Real Q(int k, nk4Index index) const;
 
  Real R(int k, const DiracSpinor &a, const DiracSpinor &b,
         const DiracSpinor &c, const DiracSpinor &d) const;
  Real R(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d) const;
 
  Real P(int k, const DiracSpinor &a, const DiracSpinor &b,
         const DiracSpinor &c, const DiracSpinor &d,
         const Angular::SixJTable *const sj = nullptr) const;
  Real P(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d,
         const Angular::SixJTable *const sj = nullptr) const;
 
  Real P2(int k, const DiracSpinor &a, const DiracSpinor &b,
          const DiracSpinor &c, const DiracSpinor &d,
          const Angular::SixJTable &sj, const std::vector<double> &fk) const;
 
  Real W(int k, const DiracSpinor &a, const DiracSpinor &b,
         const DiracSpinor &c, const DiracSpinor &d,
         const Angular::SixJTable *const sj = nullptr) const;
  Real W(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d,
         const Angular::SixJTable *const sj = nullptr) const;
 
  Real g(const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c,
         const DiracSpinor &d, int tma, int tmb, int tmc, int tmd) const;
 
  nk4Index NormalOrder(const DiracSpinor &a, const DiracSpinor &b,
                       const DiracSpinor &c, const DiracSpinor &d) const;
  nk4Index NormalOrder(nkIndex a, nkIndex b, nkIndex c, nkIndex d) const;
 
  bool is_NormalOrdered(nkIndex a, nkIndex b, nkIndex c, nkIndex d) const {
    return NormalOrder(a, b, c, d) == CurrentOrder(a, b, c, d);
  }
  bool is_NormalOrdered(const DiracSpinor &a, const DiracSpinor &b,
                        const DiracSpinor &c, const DiracSpinor &d) const {
    return NormalOrder(a, b, c, d) == CurrentOrder(a, b, c, d);
  }
 
  void write(const std::string &fname) const;
  bool read(const std::string &fname);
 
  double *get(int k, nk4Index index);
 
  nk4Index CurrentOrder(const DiracSpinor &a, const DiracSpinor &b,
                        const DiracSpinor &c, const DiracSpinor &d) const;
 
  nk4Index CurrentOrder(nkIndex a, nkIndex b, nkIndex c, nkIndex d) const;
 
  static nk4Index FormIndex(nkIndex a, nkIndex b, nkIndex c, nkIndex d);
 
  std::array<nkIndex, 4> UnFormIndex(const nk4Index &index) const;
 
private:
  // Returns the "NormalOrder" nk4Index for given set. Implemented by
  // specific symmetries
  static inline nk4Index NormalOrder_impl(nkIndex a, nkIndex b, nkIndex c,
                                          nkIndex d);
};
 
using QkTable = CoulombTable<Symmetry::Qk>;
using WkTable = CoulombTable<Symmetry::Wk>;
using LkTable = CoulombTable<Symmetry::Lk>;
using NkTable = CoulombTable<Symmetry::none>;
 
} // namespace Coulomb
 
#include "QkTable.ipp"