TDHF.hpp

#pragma once
#include "CorePolarisation.hpp"
#include <string>
#include <vector>
class DiracSpinor;
namespace DiracOperator {
class TensorOperator;
}
namespace MBPT {
class CorrelationPotential;
}
 
namespace HF {
class HartreeFock;
class Breit;
} // namespace HF
 
namespace ExternalField {
 
/*!
  @brief Uses TDHF to include core-polarisation (RPA) corrections to matrix
  elements of an external field operator.
 
  @details
  Solves the TDHF equations for each core orbital \f$ \phi_a \f$,
  \f[
    (h_{\rm HF} - \en_a \mp \omega)\varphi^a_\pm
      = -(t_\pm + \delta V_\pm - \delta\en^a_\pm)\phi_a,
  \f]
  self-consistently to determine \f$ \delta V_\pm \f$.
  See the @ref ExternalField namespace documentation for full physics description.
 
  Each core orbital acquires both a forward (\f$ \varphi^a_+ \f$, stored as X)
  and backward (\f$ \varphi^a_- \f$, stored as Y) correction. Both contribute
  to \f$ \delta V_\pm \f$:
  \f[
    \delta V_\pm \phi_i =
    \sum_a^{\rm core} \left[
      \matel{\phi_a}{Q}{\varphi^a_+}\phi_i - \matel{\phi_a}{Q}{\phi_i}\varphi^a_+
    + \matel{\varphi^a_-}{Q}{\phi_a}\phi_i - \matel{\varphi^a_-}{Q}{\phi_i}\phi_a
    \right].
  \f]
  It is via \f$ \delta V_\pm \f$ that the \f$ e^{\pm i\omega t} \f$ terms are mixed.
 
  \par Construction
  Requires a pointer to the operator (\f$ t_+ \f$, at \f$ +|\omega| \f$)
  and a const @ref HF::HartreeFock object. The user is responsible for
  updating the operator frequency externally before calling solve_core().
  The backward operator \f$ t_- \f$ (at \f$ -|\omega| \f$) is
  \f$ t_+ \f$ times TensorOperator::h_minus_sign(): -1 for an operator
  odd in the frequency (E1v), +1 otherwise.
 
  \par Usage
  @ref solve_core (omega) solves the TDHF equations for a given frequency.
  @ref dV (Fa, Fb) then returns the RPA correction to the reduced matrix element
  \f$ \redmatel{a}{\delta V}{b} \f$.
*/
class TDHF : public CorePolarisation {
 
protected:
  // dPhi = X exp(-iwt) + Y exp(+iwt)
  // (H - e - w)X = -(h + dV - de)Phi
  // (H - e + w)Y = -(h* + dV* - de)Phi
  // X_c = sum_x X_x,
  // j(x)=j(c)-k,...,j(c)+k.  And: pi(x) = pi(c)*pi(h)
  std::vector<std::vector<DiracSpinor>> m_X{};
  std::vector<std::vector<DiracSpinor>> m_Y{};
  std::vector<std::vector<DiracSpinor>> m_hFcore{};
  std::vector<std::vector<DiracSpinor>> m_hFcore_minus{};
  // can just write these to disk! Read them in, continue as per normal
 
  const HF::HartreeFock *const p_hf;
  const std::vector<DiracSpinor> m_core;
  const double m_alpha;
  const HF::Breit *const p_VBr;
  // If true, eps is the relative change |dPsi| (sqrt of the ratio); if false,
  // the squared ratio. See eps_dPsi().
  bool m_eps_sqrt{false};
  // Method of solve_core(): Anderson mixing if true, damped iteration if
  // false; the other is the fallback
  bool m_use_anderson{true};
  // Method of the core mixed-state solves: Anderson mixing
  // (solveMixedState_Anderson) if true, damped iteration (solveMixedState)
  // if false
  bool m_use_anderson_ms{true};
 
public:
  /*!
    @brief Constructs TDHF for operator h.
 
    @param h_plus  The operator \f$ t_+ \f$; must be set to positive
                   frequency before each call to solve_core().
    @param hf      @ref HF::HartreeFock object defining the core.
  */
  TDHF(const DiracOperator::TensorOperator *const h_plus,
       const HF::HartreeFock *const hf);
 
  /*!
    @brief Solves TDHF equations self-consistently for core electrons at frequency omega.
    @details
    Iterates the TDHF equations until \f$ \delta V_\pm \f$ converges to within
    eps_target(), or @p max_its iterations have been performed.
    Can be re-run at a different frequency without restarting from scratch.
 
    Solved by solve_core_anderson() (default) or solve_core_damped(). If
    that does not converge, the solve is started again from scratch with the
    other method, and the better result is kept.
 
    @param omega    External-field frequency \f$ \omega \f$ in atomic units.
    @param max_its  Maximum number of iterations.
                   Set to 1 to get the first-order correction (no damping on
                   first iteration).
    @param print    If true, write convergence progress to screen.
 
    @note Solved at \f$ |\omega| \f$: both \f$ \varphi^a_+ \f$ and
    \f$ \varphi^a_- \f$ are found simultaneously (see class docs), so a
    single call covers both components, \f$ \delta V(\pm|\omega|) \f$;
    see dV().
 
    --
 
    @note Does not update the frequency of the operator itself; for frequency-dependent
    operators, update the operator frequency externally before calling.
  */
  virtual void solve_core(double omega, int max_its = 100,
                          bool print = true) override;
 
  /*!
    @brief As solve_core(), by damped iteration only (no fallback).
    @details
    Each iteration mixes the new corrections with the previous ones, with
    damping factor eta().
 
    @note nb: if the iterations stall, the solve is restarted from scratch with
    different damping (weaker, then stronger), up to a few times, each with
    a further max_its iterations; the best attempt is kept. eta() itself is
    not changed.
  */
  void solve_core_damped(double omega, int max_its = 128, bool print = true);
 
  /*!
    @brief As solve_core(), by Anderson mixing only (no fallback).
    @details
    Self-consistency is driven by Anderson mixing (@ref anderson_coefficients)
    of all the corrections together, X and Y. Each iteration applies the
    TDHF map undamped, then extrapolates over the recent iterates. The TDHF
    equations are linear in the corrections, so this is equivalent to
    truncated GMRES: there is no damping parameter, and it converges
    wherever the equations are non-singular.
  */
  void solve_core_anderson(double omega, int max_its = 128, bool print = true);
 
  //! Method of solve_core(): Anderson mixing if true (default), damped
  //! iteration if false; the other is the fallback
  void set_anderson(bool andersonQ) { m_use_anderson = andersonQ; }
 
  virtual Method method() const override { return Method::TDHF; }
 
  virtual void clear() override;
 
  //! Returns reduced matrix element <a||dV||b>, or the conjugate
  //! <a||dV^dag||b> if conj=true.
  double dV(const DiracSpinor &Fa, const DiracSpinor &Fb,
            bool conj = false) const override;
 
  DiracSpinor dV_rhs(int kappa_n, const DiracSpinor &Fm,
                     bool conj = false) const override;
 
  //! Returns const ref to dPsi orbitals for given core orbital Fc.
  const std::vector<DiracSpinor> &get_dPsis(const DiracSpinor &Fc,
                                            dPsiType XorY) const;
  //! Returns const ref to dPsi orbital of given kappa.
  const DiracSpinor &get_dPsi_x(const DiracSpinor &Fc, dPsiType XorY,
                                const int kappa_x) const;
 
  /*!
    @brief Forms the mixed states phi^v_+- for valence state Fv (including
    core pol.): single kappa channel.
    @details
    Solves
    \f[ (h_{\rm HF} + \Sigma - \en_v - \omega)\varphi^v_+
      = -(t + \delta V - \delta\en^v)\phi_v \f]
    or
    \f[ (h_{\rm HF} + \Sigma - \en_v + \omega)\varphi^v_-
      = -(t^\dagger + \delta V^\dagger - \delta\en^v)\phi_v \f]
    Returns \f$ \chi_\beta \f$ for given kappa_beta, where
    \f[ X_{j,m} = (-1)^{j_\beta-m}tjs(j,k,j;-m,0,m)\chi_j \f]
 
    @param Fv         Valence state \f$\phi_v\f$.
    @param omega      Perturbation frequency \f$\omega\f$.
    @param XorY       Selects X or Y solution; see @ref dPsiType.
    @param kappa_beta Kappa quantum number of the target channel.
    @param Sigma      Optional correlation potential; see @ref MBPT::CorrelationPotential.
    @param st         Bra or ket convention; see @ref StateType.
    @param incl_dV    Include the induced potential \f$\delta V\f$ if true.
  */
  DiracSpinor
  solve_dPsi(const DiracSpinor &Fv, const double omega, dPsiType XorY,
             const int kappa_beta,
             const MBPT::CorrelationPotential *const Sigma = nullptr,
             StateType st = StateType::ket, bool incl_dV = true) const;
 
  //! Forms the mixed states phi^v_+- for all kappa channels; see solve_dPsi.
  std::vector<DiracSpinor>
  solve_dPsis(const DiracSpinor &Fv, const double omega, dPsiType XorY,
              const MBPT::CorrelationPotential *const Sigma = nullptr,
              StateType st = StateType::ket, bool incl_dV = true) const;
 
  // //! Writes dPsi (f-component) to textfile
  // void print(const std::string &ofname = "dPsi.txt") const;
 
protected:
  // Forms set of h*Fc for all core orbitals and all projections
  std::vector<std::vector<DiracSpinor>>
  form_hFcore(const DiracOperator::TensorOperator *h) const;
  // The same at -|omega| (for the Y equations), from the set at +|omega|
  std::vector<std::vector<DiracSpinor>>
  form_hFcore_minus(const std::vector<std::vector<DiracSpinor>> &hFcore) const;
 
  // dV_rhs built from explicit correction sets X, Y (indexed as m_X, m_Y)
  // in place of the members; dV_rhs delegates here with (m_X, m_Y). Lets a
  // derived class apply the same (real-linear) builder to a second pair of
  // sets, e.g. the imaginary parts of complex corrections.
  DiracSpinor dV_rhs_sets(int kappa_n, const DiracSpinor &Fa, bool conj,
                          const std::vector<std::vector<DiracSpinor>> &X,
                          const std::vector<std::vector<DiracSpinor>> &Y) const;
 
private:
  void initialise_dPsi();
 
  // Single iteration of TDHF equations
  std::pair<double, std::string> tdhf_core_it(double omega, double eta_damp);
  // Solves the MS equations for all projections, single core state
  void solve_ms_core(std::vector<DiracSpinor> &dFb, const DiracSpinor &Fb,
                     const std::vector<DiracSpinor> &hFbs, const double omega,
                     dPsiType XorY, double eps_ms = 1.0e-9) const;
  // As solve_ms_core(), but a single channel (one projection). Thread-safe;
  // used to parallelise tdhf_core_it() over (orbital x channel x X/Y).
  void solve_ms_core_b(DiracSpinor &dF_beta, const DiracSpinor &Fb,
                       const DiracSpinor &hFb, const double omega,
                       dPsiType XorY, double eps_ms = 1.0e-9) const;
 
  // Convergence (eps): the relative L2 change of the (undamped) X spinors,
  // Sum|dX|^2 / Sum|X_new|^2 (summed over all channels); returns its sqrt
  // -- the relative change -- if @p relative. Returns {eps, worst-channel}.
  std::pair<double, std::string>
  eps_dPsi(const std::vector<std::vector<DiracSpinor>> &Xnew,
           bool relative) const;
 
  // True if |omega| exceeds the ionisation threshold of the core, where the
  // core equations have no bound solution
  bool above_ionisation_threshold(double omega) const;
 
  // The corrections X and Y as one real vector, and back: the Anderson
  // iterate of solve_core_anderson()
  std::vector<double> dPsi_vector() const;
  void set_dPsi_vector(const std::vector<double> &state);
 
public:
  TDHF &operator=(const TDHF &) = delete;
  TDHF(const TDHF &) = default;
  ~TDHF() = default;
};
 
} // namespace ExternalField