CSF.hpp

#pragma once
#include "LinAlg/include.hpp"
#include "Wavefunction/DiracSpinor.hpp"
#include <array>
#include <optional>
#include <utility>
#include <vector>
 
namespace CI {
 
//==============================================================================
class CSF2 {
  int m_parity;
 
public:
  // nb: array of states is always sorted
  std::array<DiracSpinor::Index, 2> states;
 
  CSF2(const DiracSpinor &a, const DiracSpinor &b);
 
  DiracSpinor::Index state(std::size_t i) const;
 
  friend bool operator==(const CSF2 &A, const CSF2 &B);
  friend bool operator!=(const CSF2 &A, const CSF2 &B);
 
  static int num_different(const CSF2 &A, const CSF2 &B);
 
  static std::array<DiracSpinor::Index, 2> diff_1_na(const CSF2 &V,
                                                     const CSF2 &X);
  static DiracSpinor::Index same_1_j(const CSF2 &A, const CSF2 &B);
 
  int parity() const;
 
  std::string config(bool relativistic = false) const;
};
 
//==============================================================================
std::vector<CSF2> form_CSFs(int twoJ, int parity,
                            const std::vector<DiracSpinor> &cisp_basis);
 
//==============================================================================
struct ConfigInfo {
  // configuration (typically uses non-rel notation)
  std::string config{};
  // square of config coeficient (if non-rel, sum of all)
  double ci2{0.0};
  double gJ{0.0};
  double L{-1.0};
  double twoS{-1.0};
};
 
//==============================================================================
class PsiJPi {
 
  int m_twoj{-1};
  int m_pi{0};
 
  // Number of solutions stored:
  std::size_t m_num_solutions{0};
  // List of CSFs
  std::vector<CSF2> m_CSFs{};
  // Energy, and CI expansion coeficients
  std::pair<LinAlg::Vector<double>, LinAlg::Matrix<double>> m_Solution{};
  std::vector<ConfigInfo> m_Info{};
 
public:
  PsiJPi(int twoJ, int pi, const std::vector<DiracSpinor> &cisp_basis)
      : m_twoj(twoJ), m_pi(pi), m_CSFs(form_CSFs(twoJ, pi, cisp_basis)) {}
 
  PsiJPi() {}
 
  void solve(const LinAlg::Matrix<double> &Hci, int num_solutions = 0,
             std::optional<double> all_below = {});
 
  void update_config_info(std::size_t i, const ConfigInfo &info);
 
  const std::vector<CSF2> &CSFs() const;
 
  const CSF2 &CSF(std::size_t i) const;
 
  double energy(std::size_t i) const;
 
  LinAlg::View<const double> coefs(std::size_t i) const;
 
  double coef(std::size_t i, std::size_t j) const;
 
  int parity() const;
 
  int twoJ() const;
 
  std::size_t num_solutions() const;
 
  const ConfigInfo &info(std::size_t i) const;
};
 
} // namespace CI