HartreeFock.hpp

#pragma once
#include "Coulomb/YkTable.hpp"
#include "HF/Breit.hpp"
#include "Physics/PhysConst_constants.hpp"
#include "Potentials/Parametric_potentials.hpp"
#include "Potentials/RadPot.hpp"
#include <memory>
#include <string>
#include <vector>
class Wavefunction;
class DiracSpinor;
class Grid;
namespace MBPT {
class CorrelationPotential;
}
 
namespace HF {
 
//==============================================================================
struct EpsIts {
  double eps{0.0};
  int its{0};
  std::string symbol{};
  friend bool operator<(const EpsIts &l, const EpsIts &r) {
    return l.eps < r.eps;
  }
};
 
//==============================================================================
 
enum class Method { HartreeFock, ApproxHF, Hartree, KohnSham, Local };
Method parseMethod(const std::string &in_method);
std::string parseMethod(const Method &in_method);
 
//==============================================================================
std::vector<double> vex_approx(const DiracSpinor &Fa,
                               const std::vector<DiracSpinor> &core,
                               int k_cut = 99, double lambda_cut = 0.003);
 
DiracSpinor vexFa(const DiracSpinor &Fa, const std::vector<DiracSpinor> &core,
                  int k_cut = 99);
 
//==============================================================================
//==============================================================================
//==============================================================================
 
class HartreeFock {
 
private:
  std::shared_ptr<const Grid> m_rgrid;
  std::vector<DiracSpinor> m_core;
  std::vector<double> m_vnuc;
  std::optional<QED::RadPot> m_vrad;
  std::optional<HF::Breit> m_VBr;
  double m_alpha;
  Method m_method;
  double m_eps_HF;
  std::vector<double> m_vdir;
  Coulomb::YkTable m_Yab;
  int m_max_hf_its = 128;
 
public:
 
  HartreeFock(std::shared_ptr<const Grid> rgrid, std::vector<double> vnuc,
              std::vector<DiracSpinor> core,
              std::optional<QED::RadPot> vrad = std::nullopt,
              double m_alpha = PhysConst::alpha,
              Method method = Method::HartreeFock, double x_Breit = 0.0,
              double eps_HF = 0.0,
              Parametric::Type potential = Parametric::Type::Green,
              double H_g = 0.0, double d_t = 0.0);
 
  EpsIts solve_core(bool print = true);
 
  void
  solve_valence(std::vector<DiracSpinor> *valence, bool print = true,
                const MBPT::CorrelationPotential *const Sigma = nullptr) const;
 
  EpsIts hf_valence(DiracSpinor &Fv,
                    const MBPT::CorrelationPotential *const Sigma = nullptr,
                    std::optional<double> eta = std::nullopt,
                    std::optional<int> prev_its = std::nullopt) const;
 
  EpsIts hf_valence_Green(DiracSpinor &Fa,
                          const MBPT::CorrelationPotential *const Sigma) const;
 
  double calculateCoreEnergy() const;
 
  DiracSpinor vexFa(const DiracSpinor &Fa) const {
    // calls static version with HF core
    return ::HF::vexFa(Fa, m_core, 99);
  }
 
  DiracSpinor VBr(const DiracSpinor &Fv) const;
 
  //---------------------------
 
  const Grid &grid() const { return *m_rgrid; };
  std::shared_ptr<const Grid> grid_sptr() const { return m_rgrid; };
 
  const std::vector<double> &vdir() const { return m_vdir; }
  std::vector<double> &vdir() { return m_vdir; }
 
  const std::vector<double> &vnuc() const { return m_vnuc; }
  std::vector<double> &vnuc() { return m_vnuc; }
 
  std::vector<double> Hrad_el(int l = 0) const;
  std::vector<double> Hmag(int l = 0) const;
 
  std::vector<double> vlocal(int l = 0) const;
 
  Method method() const { return m_method; }
 
  double zion() const;
 
  bool excludeExchangeQ() const {
    return m_core.empty() ||
           !(m_method == Method::HartreeFock || m_method == Method::ApproxHF);
  }
 
  const std::vector<DiracSpinor> &core() const { return m_core; }
 
  double alpha() const { return m_alpha; }
 
  // not core, for testing)
  void set_Vrad(QED::RadPot in_vrad) { m_vrad = std::move(in_vrad); } // XXX
  const QED::RadPot *Vrad() const { return m_vrad ? &*m_vrad : nullptr; }
  QED::RadPot *Vrad() { return m_vrad ? &*m_vrad : nullptr; }
 
  const HF::Breit *vBreit() const { return m_VBr ? &*m_VBr : nullptr; }
  double x_Breit() const { return m_VBr ? m_VBr->scale_factor() : 0.0; }
 
  int num_core_electrons() const;
 
private:
  // Solve Dirac equation for core states (just once) using existing vdir
  // (usually, vdir set to parametric potential beforehand)
  EpsIts solve_initial_core(const double eps);
  // Solve equations self-consistantly for core, using local method (either
  // core-Hartree or Kohn-Sham)
  EpsIts selfcon_local_core(const double eps_target_HF);
  // Solve HF equations self-consistantly for core, using approximate HF method
  EpsIts hf_approx_core(const double eps_target_HF);
  // Solve HF equations self-consistantly for core, using Hartree-Fock method
  EpsIts hartree_fock_core();
  // Solves HF equation for valence state, assuming local potential (Hartree,
  // Local, Kohn-Sham or approxHF)
  EpsIts local_valence(DiracSpinor &Fa) const;
 
  /*
    // same as hf_valence, but uses Green method
    EpsIts hf_valence_Green(
        DiracSpinor &Fv,
        const MBPT::CorrelationPotential *const Sigma = nullptr) const;
  */
 
  // Solves HF equation for given state, using non-local Green's method for
  // inhomogeneous ODE (used for hartree_fock_core()).
  // Solve Dirac Equation (Eigenvalue):
  //  (H0 + Vl + Vx)Fa = 0
  //  (H0 + Vl)Fa = -VxFa
  // Vl is local (e.g., Vnuc + fVdir), Vx is non-local (e.g., (1-f)Vdir + Vex)
  // where v0 = (1-f)Vdir  [f=1 for valence states!, so v0 may be empty]
  // Vx also includes Breit, and Sigma
  // Small energy adjustmenets (and wfs), solve:
  // (Hl - e) dF = de * F -VxFa
  // e -> e+de, F->F+dF
  // Core is input so can call in a thread-safe way! (with a 'old_core' copy)
  // Only used in dE from dF
  void hf_orbital_green(
      DiracSpinor &Fa, double en, const std::vector<double> &vl,
      const std::vector<double> &H_mag, const DiracSpinor &VxF,
      const std::vector<DiracSpinor> &static_core,
      const std::vector<double> &dv0 = {}, const HF::Breit *const VBr = nullptr,
      const MBPT::CorrelationPotential *const Sigma = nullptr) const;
 
  // Calc's Vex*Fa, for Fa in the core. Fa must be in the core
  void vex_Fa_core(const DiracSpinor &Fa, DiracSpinor &vexFa) const;
 
  // Option to re-scale diract potential so that V(r)~-zion/r at large r
  enum class ReScale { yes = true, no = false };
  // Forms direct potential
  void update_vdir(ReScale re_scale = ReScale::no);
  // Adds the additional Kohn-Sham parts to Vdir
  void add_KohnSham_vdir_addition();
 
  // Sets Vdir to be parametric potential. By default, Greens potential
  void
  set_parametric_potential(bool print = true,
                           Parametric::Type potential = Parametric::Type::Green,
                           double H_g = 0.0, double d_t = 0.0);
 
  // Forms approximate vex for all core states
  void form_approx_vex_core(std::vector<std::vector<double>> &vex) const;
  std::vector<std::vector<double>> form_approx_vex_core() const;
  // Forms approximate vex for given core states
  void form_approx_vex_core_a(const DiracSpinor &Fa,
                              std::vector<double> &vex_a) const;
  std::vector<double> form_approx_vex_core_a(const DiracSpinor &Fa) const;
 
  // Energy guess for core states
  double enGuessCore(int n, int ka) const;
  // Energy guess for valence states
  double enGuessVal(int n, int ka) const;
};
 
} // namespace HF