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 <cmath>
#include <memory>
#include <optional>
#include <string>
#include <vector>
class Wavefunction;
class DiracSpinor;
class Grid;
namespace MBPT {
class CorrelationPotential;
}
 
//! Functions and classes for Hartree-Fock
namespace HF {
 
//==============================================================================
/*!
  @brief Convergence results of for each self-consistent-field solve.
  @details
  Holds the output of solving for an orbital (core or valence) via an
  iterative HF-like method: the convergence metric, iteration count, and
  orbital identification symbol for reporting.
  @param eps Convergence metric: relative change in orbital norm, or similar,
             from the previous iteration.
  @param its Number of iterations performed.
  @param symbol Orbital identification (e.g. "2p-", "1s+") for labelling
                output tables.
*/
struct EpsIts {
  double eps{0.0};
  int its{0};
  std::string symbol{};
  // nan (failed solve) is ordered as the largest, i.e., worst, eps
  friend bool operator<(const EpsIts &l, const EpsIts &r) {
    return !std::isnan(l.eps) &&
           (std::isnan(r.eps) || std::abs(l.eps) < std::abs(r.eps));
  }
};
 
//==============================================================================
/*! @brief Methods available for self-consistant field model
  @details
 - HartreeFock: Self-consistent Hartree-Fock method
 - ApproxHF   : Approximate (localised) Hartree-Fock method
 - Hartree    : Core-Hartree method. No exchange, Vdir includes self-interaction
 - KohnSham   : Kohn-Sham (Density functional), includes Latter correction
 - Local      : Uses a local parameteric potential. [NOT self-consistant field]
 */
enum class Method { HartreeFock, ApproxHF, Hartree, KohnSham, Local };
 
//! Convers string (name) of method (e.g., HartreeFock) to HF::Method enum
Method parseMethod(const std::string &in_method);
//! Convers HF::Method enum to string (name) of method (e.g., HartreeFock)
std::string parseMethod(const Method &in_method);
//! Convers HF::Method enum to short string (name) of method (e.g., HF)
std::string parseMethod_short(const Method &in_method);
 
//==============================================================================
//! Forms approx (localised) exchange potential, from scratch
//! @details Needs existing orbital Fa, and the core orbitals.
//! k_cut is max multipolarity to sum over for exchange term [can limit to ~1
//! (e.g.) for speed when high accuracy is not required]
std::vector<double> vex_approx(const DiracSpinor &Fa,
                               const std::vector<DiracSpinor> &core,
                               int k_cut = 99, double lambda_cut = 0.003);
 
//! @brief Calculates V_exch * Fa, for any orbital Fa (calculates Coulomb
//! integral from scratch).
//! @details  k_cut is max multipolarity to sum over for exchange term [can
//! limit to ~1 (e.g.) for speed when high accuracy is not required]
DiracSpinor vexFa(const DiracSpinor &Fa, const std::vector<DiracSpinor> &core,
                  int k_cut = 99);
 
//! @brief Density-based (Kohn-Sham/Slater) local exchange potential, ~rho^1/3.
//! @details Does not depend on the orbital, so it
//! conditions all symmetries uniformly -- useful as the local-exchange term on
//! the LHS of the mixed-states iteration (where vex_approx, which divides by
//! the perturbation, is poorly conditioned). Grid is taken from @p core.
std::vector<double> vex_KS(const std::vector<DiracSpinor> &core);
 
/*!
  @brief Exchange of a single electron in orbital Fa, acting on Fv: X_a Fv.
  @details
 
  \f[ X_a F_v = -\frac{1}{[j_v]\,[j_a]}\sum_k C^k_{va}{}^2\, y^k_{av}(r)\,F_a \f]
 
  - The b=Fa term of vexFa(), scaled to ONE electron (weight 1/[j_a] rather
    than the subshell occupancy): the spherically-averaged one-electron
    self-exchange.
  - Used for the V^{N-1} (residual ion) potential: removing one electron from
    subshell a removes y^0_aa from the direct potential and removes this term
    from the exchange, V^{N-1}_a = V_HF - (y^0_aa + X_a).
 
  @note All multipoles k allowed by the (v,a) triangle/parity rules are
        included; the occupation fraction of Fa is NOT applied (exactly one
        electron is removed).
*/
DiracSpinor vexFa_1el(const DiracSpinor &Fv, const DiracSpinor &Fa);
 
//==============================================================================
//==============================================================================
//==============================================================================
 
//! Solves relativistic Hartree-Fock equations for core and valence. Optionally
//! includes Breit and QED effects. Can include Sigma (correlations) for valence
//! states. Class stores nuc. and direct potentials, a set of yk integrals, and
//! QED potential. Stores the core orbitals.
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 = 256;
 
public:
  //! @brief Method is enum class, eps_HF is convergence goal.
  /*! @details
    Required:
      - rgrid: Radial grid (shared pointer)
        - (This is required to allow no core orbitals)
        - Assumed to be same grid as for core orbitals
      - vnuc - nuclear potential
        - Assumed to be same length as radial grid
        - A copy is stored. May be updated
    Optional:
      - alpha (fine structure constant). default = true value
      - method default = HartreeFock
      - breit_params - Breit scaling factors. If std::nullopt (default), no Breit.
        If present, a Breit object is constructed from the params.
      - eps_HF: convergence goal
      - potential: which parametric potential used for initial Potential
      - h and d (or g and t) are parameters for above (if left zero, default
        will be chosen)
      - Note:  Parametric::Type potential (and parameters H,d) are for
    initial approx. Usually doesn't matter at all, and defaults should be used.
    If using local potential [method=Local], these are the final parameters. If
    any are set to zero - will be looked up.
  */
  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,
              std::optional<Breit::Params> breit_params = std::nullopt,
              double eps_HF = 0.0,
              Parametric::Type potential = Parametric::Type::Green,
              double H_g = 0.0, double d_t = 0.0);
 
  //! Solves HF equations self-consitantly for core orbs. Returns epsilon.
  EpsIts solve_core(bool print = true);
 
  //! Solves HF for given valence list. They need not already be solutions.
  //! @details Note: If given energy is set to zero, states assumed to not be
  //! existing solutions; initial energy is guessed and solved from scratch. If
  //! initial energy is non-zero, that energy is used and states are assumed to
  //! already be (approximate) solutions.
  void
  solve_valence(std::vector<DiracSpinor> *valence, bool print = true,
                const MBPT::CorrelationPotential *const Sigma = nullptr) const;
 
  //! Solves HF equation (+ Sigma) for single valence state.
  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;
 
  //! Solves HF equation (+ Sigma) for single valence state, alternative method
  EpsIts hf_valence_Green(DiracSpinor &Fa,
                          const MBPT::CorrelationPotential *const Sigma) const;
 
  //! Calculates the HF core energy (including Breit, if present)
  double calculateCoreEnergy() const;
 
  //! Calculates exchange term Vex*Fa
  DiracSpinor vexFa(const DiracSpinor &Fa) const {
    // calls static version with HF core
    return ::HF::vexFa(Fa, m_core, 99);
  }
 
  //! Breit interaction V_Br*Fa
  DiracSpinor VBr(const DiracSpinor &Fv) const;
 
  //---------------------------
 
  //! Resturns a const reference to the radial grid
  const Grid &grid() const { return *m_rgrid; };
  //! Resturns copy of shared_ptr to grid [shared resource] - used when we want
  //! to construct a new object that shares this grid
  std::shared_ptr<const Grid> grid_sptr() const { return m_rgrid; };
 
  //! Returns reference to Vdir (direct HF potential)
  const std::vector<double> &vdir() const { return m_vdir; }
  std::vector<double> &vdir() { return m_vdir; }
 
  //! Returns reference to Vnuc (nuclear potential)
  const std::vector<double> &vnuc() const { return m_vnuc; }
  std::vector<double> &vnuc() { return m_vnuc; }
 
  //! Electric part of radiative potential
  std::vector<double> Hrad_el(int l = 0) const;
  //! Magnetic (off-diagonal) part of radiative potential. Doesn't currently
  //! depend on l
  std::vector<double> Hmag(int l = 0) const;
 
  //! vlocal = vnuc + vrad_el + vdir
  std::vector<double> vlocal(int l = 0) const;
 
  //! Which method used to solve HF
  Method method() const { return m_method; }
 
  //! Effective charge at large Z : zion = Z - num_core_electrons
  double zion() const;
 
  //! Returns true if exchange not included
  bool is_localQ() const {
    return m_core.empty() ||
           !(m_method == Method::HartreeFock || m_method == Method::ApproxHF);
  }
 
  //! vector of core orbitals
  const std::vector<DiracSpinor> &core() const { return m_core; }
 
  //! Value of fine-structure constant used
  double alpha() const { return m_alpha; }
 
  //! Update the Vrad used inside HF (only used if we want QED into valence but
  // not core, for testing)
  void set_Vrad(QED::RadPot in_vrad) { m_vrad = std::move(in_vrad); } // XXX
  //! Get (const) ptr to Vrad - may be null
  const QED::RadPot *Vrad() const { return m_vrad ? &*m_vrad : nullptr; }
  QED::RadPot *Vrad() { return m_vrad ? &*m_vrad : nullptr; }
 
  //! pointer to Breit - may be nullptr if no breit
  const HF::Breit *vBreit() const { return m_VBr ? &*m_VBr : nullptr; }
  //! Breit scale factor (usualy 0 or 1)
  double x_Breit() const { return m_VBr ? m_VBr->scale_factor() : 0.0; }
 
  //! Number of electrons in the core
  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(bool print);
  // Solves HF equation for valence state, assuming local potential (Hartree,
  // Local, Kohn-Sham or approxHF)
  EpsIts local_valence(DiracSpinor &Fa) const;
 
  // Prints a warning for each orbital that does not have n-l-1 nodes (i.e.,
  // converged to the wrong state)
  void check_nodes(const std::vector<DiracSpinor> &orbitals,
                   const std::string &label) 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