CoulombIntegrals.hpp

#pragma once
#include "Angular/include.hpp"
#include "Wavefunction/DiracSpinor.hpp"
#include <algorithm>
#include <cstdlib>
#include <optional>
#include <utility>
#include <vector>
 
//! Functions (+classes) for computing Coulomb integrals
namespace Coulomb {
 
//! Calculates Hartree Screening functions \f$y^k_{ab}(r)\f$
//! @details maxi is max point to calculate; blank or zero means all the way
std::vector<double> yk_ab(const int k, const DiracSpinor &Fa,
                          const DiracSpinor &Fb, const std::size_t maxi = 0);
 
//! Overload: does not allocate ykab
void yk_ab(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
           std::vector<double> &ykab, const std::size_t maxi = 0);
 
//==============================================================================
 
//! Calculates R^k_abcd for given k. From scratch (calculates y)
double Rk_abcd(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
               const DiracSpinor &Fc, const DiracSpinor &Fd);
 
//! Overload for when y^k_bd already exists [much faster]
double Rk_abcd(const DiracSpinor &Fa, const DiracSpinor &Fc,
               const std::vector<double> &ykbd);
 
//! "Right-hand-side" R^k{v}_bcd [i.e., without Fv integral]
DiracSpinor Rkv_bcd(const int k, const int kappa_v, const DiracSpinor &Fb,
                    const DiracSpinor &Fc, const DiracSpinor &Fd);
 
//! Overload for when y^k_bd already exists [much faster]
DiracSpinor Rkv_bcd(const int kappa_v, const DiracSpinor &Fc,
                    const std::vector<double> &ykbd);
 
//! Overload for when spinor exists. Rkv is overwritten
void Rkv_bcd(DiracSpinor *const Rkv, const DiracSpinor &Fc,
             const std::vector<double> &ykbd);
 
//==============================================================================
 
//! Calculates Q^k_abcd for given k. From scratch (calculates y) [see YkTable
//! version if already have YkTable]
double Qk_abcd(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
               const DiracSpinor &Fc, const DiracSpinor &Fd);
 
//! Calculates Q^k(v)_bcd for given k,kappa_v. From scratch (calculates y) [see
//! YkTable version if already have YkTable]
DiracSpinor Qkv_bcd(const int k, int kappa_v, const DiracSpinor &Fb,
                    const DiracSpinor &Fc, const DiracSpinor &Fd);
 
//==============================================================================
//! Calculates g from scratch - not used often
double g_abcd(const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c,
              const DiracSpinor &d, int tma, int tmb, int tmc, int tmd);
 
//==============================================================================
//! Just selection rule for Qk_abcd
bool Qk_abcd_SR(int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
                const DiracSpinor &Fc, const DiracSpinor &Fd);
//! Just selection rule for Qk_abcd
bool Qk_abcd_SR(int k, int ka, int kb, int kc, int kd);
 
//! Just selection rule for Pk_abcd
bool Pk_abcd_SR(int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
                const DiracSpinor &Fc, const DiracSpinor &Fd);
 
//! Just selection rule for Pk_abcd
bool Pk_abcd_SR(int k, int ka, int kb, int kc, int kd);
 
//==============================================================================
 
//! Exchange only version of W (W-Q): W = Q + P [see Qk above]
double Pk_abcd(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
               const DiracSpinor &Fc, const DiracSpinor &Fd);
 
//! Exchange only version of W (W-Q): W = Q + P [see Qk above]
DiracSpinor Pkv_bcd(const int k, int kappa_v, const DiracSpinor &Fb,
                    const DiracSpinor &Fc, const DiracSpinor &Fd);
 
//==============================================================================
 
//! Calculates W^k_abcd for given k. From scratch (calculates y)
/*! @details
  \f[ W^k_{abcd} = Q^k_{abcd} + \sum_l [k]
      \begin{Bmatrix}a&c&k\\b&d&l\end{Bmatrix} * Q^l_{abdc} \f]
 \f[ W^k_{abcd} = Q^k_{abcd} + P^k_{abcd}  \f]
 */
double Wk_abcd(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
               const DiracSpinor &Fc, const DiracSpinor &Fd);
 
DiracSpinor Wkv_bcd(const int k, int kappa_v, const DiracSpinor &Fb,
                    const DiracSpinor &Fc, const DiracSpinor &Fd);
 
//==============================================================================
 
//! Returns min and max k (multipolarity) allowed for Triangle(k,a,b),
//! NOT accounting for parity (2j only, not kappa/l) (used by k_minmax_P,W)
//! @details made inline - called very often
inline std::pair<int, int> k_minmax_tj(int tja, int tjb) {
  return std::make_pair(std::abs(tja - tjb) / 2, (tja + tjb) / 2);
}
 
//! Returns min and max k (multipolarity) allowed for C^k_ab, accounting for
//! parity (used by k_minmax_Q)
//! @details made inline - called very often
inline std::pair<int, int> k_minmax_Ck(const DiracSpinor &a,
                                       const DiracSpinor &b) {
  auto minmax = k_minmax_tj(a.twoj(), b.twoj());
  auto &min_k = minmax.first;
  auto &max_k = minmax.second;
  if ((a.l() + b.l() + min_k) % 2 != 0) {
    ++min_k;
  }
  if ((a.l() + b.l() + max_k) % 2 != 0) {
    --max_k;
  }
  return minmax;
}
//! Returns min and max k (multipolarity) allowed for C^k_ab, accounting for
//! parity (used by k_minmax_Q)
//! @details inline (header) - hot inner-loop helper.
inline std::pair<int, int> k_minmax_Ck(int kappa_a, int kappa_b) {
  auto minmax = k_minmax_tj(Angular::twoj_k(kappa_a), Angular::twoj_k(kappa_b));
  auto &min_k = minmax.first;
  auto &max_k = minmax.second;
  const auto la = Angular::l_k(kappa_a);
  const auto lb = Angular::l_k(kappa_b);
  if ((la + lb + min_k) % 2 != 0) {
    ++min_k;
  }
  if ((la + lb + max_k) % 2 != 0) {
    --max_k;
  }
  return minmax;
}
 
//! Returns min and max k (multipolarity) allowed for Q^k_abcd.
//! Parity rule is included, so you may safely call k+=2.
//! Guaranteed to be non-zero *at* min and max, and every 2nd inbetween.
//! @details made inline - called very often
inline std::pair<int, int> k_minmax_Q(const DiracSpinor &a,
                                      const DiracSpinor &b,
                                      const DiracSpinor &c,
                                      const DiracSpinor &d) {
  // Determine if K needs to be even/odd (parity selection rule)
  const auto k_even_ac = (a.l() + c.l()) % 2 == 0;
  const auto k_even_bd = (b.l() + d.l()) % 2 == 0;
  if (k_even_ac != k_even_bd) {
    // no K satisfies selection rule!
    return {1, 0};
  }
 
  // Find min/max k from triangle rule:
  const auto [l1, u1] = k_minmax_Ck(a, c);
  const auto [l2, u2] = k_minmax_Ck(b, d);
  return {std::max(l1, l2), std::min(u1, u2)};
}
 
//! Returns min and max k (multipolarity) allowed for Q^k_abcd.
//! Parity rule is included, so you may safely call k+=2.
//! Guaranteed to be non-zero *at* min and max, and every 2nd inbetween.
//! @details made inline - called very often
inline std::pair<int, int> k_minmax_Q(int kap_a, int kap_b, int kap_c,
                                      int kap_d) {
  // Determine if K needs to be even/odd (parity selection rule)
  const auto k_even_ac = (Angular::l_k(kap_a) + Angular::l_k(kap_c)) % 2 == 0;
  const auto k_even_bd = (Angular::l_k(kap_b) + Angular::l_k(kap_d)) % 2 == 0;
  if (k_even_ac != k_even_bd) {
    // no K satisfies selection rule!
    return {1, 0};
  }
 
  // Find min/max k from triangle rule:
  const auto [l1, u1] = k_minmax_Ck(kap_a, kap_c);
  const auto [l2, u2] = k_minmax_Ck(kap_b, kap_d);
  return {std::max(l1, l2), std::min(u1, u2)};
}
 
//! Returns min and max k (multipolarity) allowed for P^k_abcd.
//! DOES NOT contain parity rules (6j only) - so NOT safe to call k+=2.
//! Guaranteed to be non-zero *at* min and max.
std::pair<int, int> k_minmax_P(const DiracSpinor &a, const DiracSpinor &b,
                               const DiracSpinor &c, const DiracSpinor &d);
 
//! Returns min and max k (multipolarity) allowed for W^k_abcd.
//! DOES NOT contain parity rules (6j only) - so NOT safe to call k+=2.
//! Cannot guarantee non-zero, since when a=b or c=d, sometimes have P=-Q.
//! But when a!=b and c!=d, guaranteed to be non-zero *at* min and max.
std::pair<int, int> k_minmax_W(const DiracSpinor &a, const DiracSpinor &b,
                               const DiracSpinor &c, const DiracSpinor &d);
 
//==============================================================================
 
//! Returns number of orbitals that are below Fermi level. Used for Qk selection
int number_below_Fermi(const DiracSpinor &i, const DiracSpinor &j,
                       const DiracSpinor &k, const DiracSpinor &l,
                       double eFermi);
 
//==============================================================================
 
template <class A>
static int twojk(const A &a) {
  if constexpr (std::is_same_v<A, DiracSpinor>) {
    return a.twoj();
  } else {
    static_assert(std::is_same_v<A, int>);
    return 2 * a;
  }
}
 
template <class A>
static std::optional<int> twojknull(const A &a) {
  if constexpr (std::is_same_v<A, DiracSpinor>) {
    return a.twoj();
  } else if constexpr (std::is_same_v<A, int>) {
    static_assert(std::is_same_v<A, int>);
    return 2 * a;
  } else {
    return std::nullopt;
  }
}
 
template <class A, class B, class C, class D, class E, class F>
static double sixj(const A &a, const B &b, const C &c, const D &d, const E &e,
                   const F &f) {
  return Angular::sixj_2(twojk(a), twojk(b), twojk(c), twojk(d), twojk(e),
                         twojk(f));
}
 
template <class A = std::optional<int>, class B = std::optional<int>,
          class C = std::optional<int>, class D = std::optional<int>,
          class E = std::optional<int>, class F = std::optional<int>>
static bool sixjTriads(const A &a, const B &b, const C &c, const D &d,
                       const E &e, const F &f) {
  return Angular::sixjTriads(twojknull(a), twojknull(b), twojknull(c),
                             twojknull(d), twojknull(e), twojknull(f));
}
 
template <class A, class B, class C>
static bool triangle(const A &a, const B &b, const C &c) {
  return Angular::triangle(twojk(a), twojk(b), twojk(c));
}
 
} // namespace Coulomb