DiracSpinor.hpp

#pragma once
#include <memory>
#include <string>
#include <utility>
#include <vector>
class Grid;
 
class DiracSpinor {
 
public:
  DiracSpinor(int in_n, int in_kappa, std::shared_ptr<const Grid> in_rgrid);
  using Index = uint16_t;
 
private:
  // Radial Grid; links F[i] to F(r)
  std::shared_ptr<const Grid> m_rgrid;
  // Principal quantum number
  int m_n;
  // Dirac quantum number, kappa
  int m_kappa;
  // Single-particle energy, not including rest energy
  double m_en = 0.0;
  // Upper (large) radial component
  std::vector<double> m_f;
  // Lower (small) radial component
  std::vector<double> m_g;
  // `practical zero': p0 is first non-zero point for f(r) [usually p0=0]
  std::size_t m_p0 = 0;
  // `practical infinity': pinf is last non-zero point for f(r)
  std::size_t m_pinf;
  // Number of iterations until energy convergence (for latest routine only)
  int m_its = -1;
  // Fractional energy convergence: eps = |(en'-en)/en|
  double m_eps = -1.0;
  // Occupation fraction. =1 for closed shells. =1/(2j+1) for valence
  double m_occ_frac = 0.0;
 
  // 2j, l, pi, kappa_index (for convenience): these are defined by n and kappa
  int m_twoj;
  int m_l;
  int m_parity;
  int m_kappa_index;
  Index m_nkappa_index;
 
public:
  int n() const { return m_n; }
  int kappa() const { return m_kappa; }
  int l() const { return m_l; }
  double jjp1() const { return 0.25 * double(m_twoj * (m_twoj + 2)); }
  int twoj() const { return m_twoj; }
  int twojp1() const { return m_twoj + 1; }
  int parity() const { return m_parity; }
  int k_index() const { return m_kappa_index; }
  Index nk_index() const { return m_nkappa_index; }
  std::string symbol(bool gnuplot = false) const;
  std::string shortSymbol() const;
 
  void set_new_kappa(int new_kappa);
 
  const Grid &grid() const { return *m_rgrid; };
  std::shared_ptr<const Grid> grid_sptr() const { return m_rgrid; };
 
  double en() const { return m_en; }
  double &en() { return m_en; }
 
  const std::vector<double> &f() const { return m_f; }
  std::vector<double> &f() { return m_f; }
  double f(std::size_t i) const { return m_f.at(i); }
  double &f(std::size_t i) { return m_f.at(i); }
 
  const std::vector<double> &g() const { return m_g; }
  std::vector<double> &g() { return m_g; }
  double g(std::size_t i) const { return m_g.at(i); }
  double &g(std::size_t i) { return m_g.at(i); }
 
  // XXX Kill this?
  auto min_pt() const { return m_p0; }
  auto &min_pt() { return m_p0; }
 
  auto max_pt() const { return m_pinf; }
  auto &max_pt() { return m_pinf; }
 
  double r0() const;
  double rinf() const;
 
  auto its() const { return m_its; }
  auto &its() { return m_its; }
 
  auto occ_frac() const { return m_occ_frac; }
  auto &occ_frac() { return m_occ_frac; }
 
  auto eps() const { return m_eps; }
  auto &eps() { return m_eps; }
 
  const DiracSpinor &scale(const double factor);
  const DiracSpinor &scale(const std::vector<double> &v);
  void normalise(double norm_to = 1.0);
  void zero_boundaries();
 
  double norm() const;
  double norm2() const;
 
  std::pair<double, double> r0pinfratio() const;
 
  std::vector<double> rho() const;
 
  int num_electrons() const;
 
public:
  // Operator overloads
 
  friend double operator*(const DiracSpinor &Fa, const DiracSpinor &Fb);
 
  DiracSpinor &operator+=(const DiracSpinor &rhs);
  DiracSpinor &operator-=(const DiracSpinor &rhs);
  friend DiracSpinor operator+(DiracSpinor lhs, const DiracSpinor &rhs);
  friend DiracSpinor operator-(DiracSpinor lhs, const DiracSpinor &rhs);
 
  DiracSpinor &operator*=(const double x);
  friend DiracSpinor operator*(DiracSpinor Fa, const double x);
  friend DiracSpinor operator*(const double x, DiracSpinor Fa);
 
  DiracSpinor &operator*=(const std::vector<double> &v);
  friend DiracSpinor operator*(const std::vector<double> &v, DiracSpinor Fa);
 
  friend bool operator==(const DiracSpinor &lhs, const DiracSpinor &rhs);
  friend bool operator!=(const DiracSpinor &lhs, const DiracSpinor &rhs);
  friend bool operator<(const DiracSpinor &lhs, const DiracSpinor &rhs);
  friend bool operator>(const DiracSpinor &lhs, const DiracSpinor &rhs);
  friend bool operator<=(const DiracSpinor &lhs, const DiracSpinor &rhs);
  friend bool operator>=(const DiracSpinor &lhs, const DiracSpinor &rhs);
 
  static bool comp_l(const DiracSpinor &lhs, const DiracSpinor &rhs) {
    return lhs.m_l < rhs.m_l;
  }
  static bool comp_j(const DiracSpinor &lhs, const DiracSpinor &rhs) {
    return lhs.m_twoj < rhs.m_twoj;
  }
  static bool comp_ki(const DiracSpinor &lhs, const DiracSpinor &rhs) {
    return lhs.m_kappa_index < rhs.m_kappa_index;
  }
  static bool comp_n(const DiracSpinor &lhs, const DiracSpinor &rhs) {
    return lhs.m_n < rhs.m_n;
  }
  static bool comp_en(const DiracSpinor &lhs, const DiracSpinor &rhs) {
    return lhs.en() < rhs.en();
  }
 
  friend std::ostream &operator<<(std::ostream &ostr, const DiracSpinor &Fa) {
    return ostr << Fa.shortSymbol();
  }
 
  // Static (helper) functions:
 
  static std::pair<double, std::string>
  check_ortho(const std::vector<DiracSpinor> &a,
              const std::vector<DiracSpinor> &b);
 
  static void orthonormaliseOrbitals(std::vector<DiracSpinor> &orbs,
                                     int num_its = 1);
 
  [[nodiscard]] static DiracSpinor
  orthogonaliseWrt(const DiracSpinor &Fin,
                   const std::vector<DiracSpinor> &orbs);
 
  [[nodiscard]] static DiracSpinor
  orthonormaliseWrt(const DiracSpinor &Fin,
                    const std::vector<DiracSpinor> &orbs);
 
  void orthog(const DiracSpinor &rhs);
 
  static std::string shortSymbol(int n, int kappa);
 
  static std::string state_config(const std::vector<DiracSpinor> &orbs);
 
  static DiracSpinor exactHlike(int n, int k, std::shared_ptr<const Grid> rgrid,
                                double zeff, double alpha = 0.0);
 
  static const DiracSpinor *find(int n, int k,
                                 const std::vector<DiracSpinor> &orbs);
 
  static int max_tj(const std::vector<DiracSpinor> &orbs);
  static int max_l(const std::vector<DiracSpinor> &orbs);
  static int max_n(const std::vector<DiracSpinor> &orbs);
  static int max_kindex(const std::vector<DiracSpinor> &orbs);
 
  static std::pair<std::vector<DiracSpinor>, std::vector<DiracSpinor>>
  split_by_energy(const std::vector<DiracSpinor> &orbitals, double energy,
                  int n_min_core = 1);
 
  static std::pair<std::vector<DiracSpinor>, std::vector<DiracSpinor>>
  split_by_core(const std::vector<DiracSpinor> &orbitals,
                const std::vector<DiracSpinor> &core, int n_min_core = 1);
};