Wavefunction.hpp

#pragma once
#include "CI/CSF.hpp"
#include "HF/HartreeFock.hpp"
#include "MBPT/CorrelationPotential.hpp"
#include "Maths/Grid.hpp"
#include "Physics/AtomData.hpp"
#include "Physics/NuclearPotentials.hpp"
#include "Physics/PhysConst_constants.hpp"
#include "Physics/RadPot.hpp"
#include "Wavefunction/BSplineBasis.hpp"
#include "Wavefunction/DiracSpinor.hpp"
#include <iostream>
#include <memory>
#include <numeric>
#include <optional>
#include <string>
#include <utility>
#include <vector>
 
namespace SplineBasis {
struct Parameters;
}
 
//==============================================================================
class Wavefunction {
 
public:
  Wavefunction(std::shared_ptr<const Grid> grid,
               const Nuclear::Nucleus &nucleus, double var_alpha = 1.0,
               const std::string &run_label = "");
  Wavefunction(const GridParameters &gridparams,
               const Nuclear::Nucleus &nucleus, double var_alpha = 1.0,
               const std::string &run_label = "");
 
  Wavefunction() : Wavefunction(GridParameters{}, Nuclear::Nucleus{}) {}
 
private:
  // Radial grid
  std::shared_ptr<const Grid> rgrid;
  // Internal value for alpha (alpha = var_alpha * alpha_0, alpha_0=~1/137)
  double m_alpha;
  std::string m_run_label;
  // Holds nuclear parameters (isotope, charge distro etc.)
  Nuclear::Nucleus m_nucleus;
  // Valence (single-particle) orbitals
  std::vector<DiracSpinor> m_valence{};
  // Basis, daigonalised over HF core. Used for MBPT
  std::vector<DiracSpinor> m_basis{};
  // Sprectrum: like basis, but includes Sigma (correlations).
  std::vector<DiracSpinor> m_spectrum{};
  // Nuclear potential // here AND hf?
  std::vector<double> m_vnuc{};
  // Hartree-Fock potential
  std::optional<HF::HartreeFock> m_HF{std::nullopt};
  // Correlation potential; for now unique_ptr; prefer std::optional
  std::optional<MBPT::CorrelationPotential> m_Sigma{};
  // Core configuration (non-rel terms)
  std::string m_core_string = "";
  std::string m_aboveFermi_core_string = "";
 
  std::vector<CI::PsiJPi> m_CIwfs{};
 
public:
  const Grid &grid() const { return *rgrid; };
  std::shared_ptr<const Grid> grid_sptr() const { return rgrid; };
 
  double alpha() const { return m_alpha; }
 
  double dalpha2() const {
    // (alpha/alpha_0)^2 -1
    return (m_alpha * m_alpha / PhysConst::alpha2) - 1.0;
  }
 
  const Nuclear::Nucleus &nucleus() const { return m_nucleus; }
  int Znuc() const { return m_nucleus.z(); }
  int Anuc() const { return m_nucleus.a(); }
  double get_rrms() const { return m_nucleus.r_rms(); }
 
  const std::vector<DiracSpinor> &core() const {
    static const auto empty = std::vector<DiracSpinor>{}; //?
    return m_HF ? m_HF->core() : empty;
  }
 
  const std::vector<DiracSpinor> &valence() const { return m_valence; }
  std::vector<DiracSpinor> &valence() { return m_valence; }
 
  const std::vector<DiracSpinor> &basis() const { return m_basis; }
  std::vector<DiracSpinor> &basis() { return m_basis; }
 
  const std::vector<DiracSpinor> &spectrum() const { return m_spectrum; }
  std::vector<DiracSpinor> &spectrum() { return m_spectrum; }
 
  const std::vector<CI::PsiJPi> &CIwfs() const { return m_CIwfs; }
 
  const CI::PsiJPi *CIwf(int J, int parity) const {
    for (const auto &ci_wf : m_CIwfs) {
      if (ci_wf.twoJ() == 2 * J && ci_wf.parity() == parity)
        return &ci_wf;
    }
    return nullptr;
  }
 
  const std::vector<double> &vnuc() const { return m_vnuc; }
 
  const HF::HartreeFock *vHF() const { return m_HF ? &*m_HF : nullptr; }
  HF::HartreeFock *vHF() { return m_HF ? &*m_HF : nullptr; }
 
  std::vector<double> vlocal(int l = 0) const;
 
  std::vector<double> Hmag(int l = 0) const;
 
  const QED::RadPot *vrad() const { return m_HF ? m_HF->Vrad() : nullptr; }
  QED::RadPot *vrad() { return m_HF ? m_HF->Vrad() : nullptr; }
 
  const MBPT::CorrelationPotential *Sigma() const {
    return m_Sigma ? &*m_Sigma : nullptr;
  }
  MBPT::CorrelationPotential *Sigma() { return m_Sigma ? &*m_Sigma : nullptr; }
 
  //----------------------------------
 
  int Ncore() const;
 
  const DiracSpinor *getState(int n, int k) const;
  const DiracSpinor *getState(std::string_view state) const;
 
  double FermiLevel() const;
 
  double energy_gap() const;
 
  std::string coreConfiguration() const { return m_core_string; }
 
  std::string coreConfiguration_nice() const {
    return AtomData::niceCoreOutput(m_core_string);
  }
 
  std::string atom() const {
    return AtomData::atomicSymbol(m_nucleus.z()) +
           ", Z=" + std::to_string(m_nucleus.z()) +
           " A=" + std::to_string(m_nucleus.a());
  }
 
  std::string atomicSymbol() const {
    return AtomData::atomicSymbol(m_nucleus.z());
  }
 
  std::string identity() const;
 
  int ion_degree(int num_val) const { return Zion() - num_val; }
 
  std::string ion_symbol(int num_val) const {
    return qip::int_to_roman(Zion() - num_val + 1);
  }
 
  int Zion() const { return Znuc() - Ncore(); }
 
  void printCore() const;
  void printValence(const std::vector<DiracSpinor> &tmp_orbitals = {}) const;
 
  void printBasis(const std::vector<DiracSpinor> &the_basis) const;
 
  bool isInCore(int n, int k) const;
  bool isInValence(int n, int k) const;
 
  std::vector<double> coreDensity() const;
 
  double coreEnergyHF() const;
 
  //------------------------------------------------------------------
 
  void set_HF(const std::string &method = "HartreeFock",
              const double x_Breit = 0.0, const std::string &in_core = "",
              double eps_HF = 1.0e-13, bool print = true);
 
  void solve_core(bool print = true);
 
  void solve_core(const std::string &method, const double x_Breit = 0.0,
                  const std::string &in_core = "", double eps_HF = 1.0e-13,
                  bool print = true);
 
  void solve_valence(const std::string &in_valence_str = "",
                     const bool print = true);
 
  void hartreeFockBrueckner(const bool print = true);
 
  void fitSigma_hfBrueckner(const std::string &valence_list,
                            const std::vector<double> &fit_energies);
 
  void SOEnergyShift();
 
  void radiativePotential(QED::RadPot::Scale s, double rcut, double scale_rN,
                          const std::vector<double> &x_spd,
                          bool do_readwrite = true, bool print = true);
 
  void radiativePotential(const IO::InputBlock &qed_input, bool do_readwrite,
                          bool print);
 
  void formBasis(const SplineBasis::Parameters &params);
 
  void formSpectrum(const SplineBasis::Parameters &params);
 
  void formSigma(int nmin_core = 1, int nmin_core_F = 1, double r0 = 1.0e-4,
                 double rmax = 30.0, int stride = 4, bool each_valence = false,
                 bool include_G = false, bool include_Breit = false,
                 const std::vector<double> &lambdas = {},
                 const std::vector<double> &fk = {},
                 const std::vector<double> &etak = {},
                 const std::string &in_fname = "",
                 const std::string &out_fname = "", bool FeynmanQ = false,
                 bool ScreeningQ = false, bool holeParticleQ = false,
                 int lmax = 6, double omre = -0.2, double w0 = 0.01,
                 double wratio = 1.5,
                 const std::optional<IO::InputBlock> &ek = std::nullopt);
 
  // void correlations(const IO::InputBlock &input);
 
  void copySigma(const MBPT::CorrelationPotential *const Sigma) {
    if (Sigma != nullptr)
      m_Sigma = *Sigma;
  }
 
  // _and_ HartreeFock)
  void update_Vnuc(const std::vector<double> &v_new) {
    m_vnuc = v_new;
    if (m_HF) {
      m_HF->vnuc() = v_new;
    }
  }
 
  std::tuple<double, double> lminmax_core_range(int l, double eps = 0.0) const;
 
  double H0ab(const DiracSpinor &Fa, const DiracSpinor &Fb) const;
  double H0ab(const DiracSpinor &Fa, const DiracSpinor &dFa,
              const DiracSpinor &Fb, const DiracSpinor &dFb) const;
 
  double Hab(const DiracSpinor &Fa, const DiracSpinor &Fb) const;
 
  void ConfigurationInteraction(const IO::InputBlock &input);
 
private:
  double H0ab_impl(const DiracSpinor &Fa, std::vector<double> dga,
                   const DiracSpinor &Fb, std::vector<double> dgb) const;
 
  // Creates set of blank core orbitals
  std::vector<DiracSpinor> determineCore(const std::string &str_core_in);
  bool isInAboveFermiCore(int n, int k) const;
};