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 {
 
 
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{};
  // 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 std::vector<double> m_Hmag;
  const double m_alpha;
  const HF::Breit *const p_VBr;
 
public:
  TDHF(const DiracOperator::TensorOperator *const h,
       const HF::HartreeFock *const hf);
 
 
  virtual void solve_core(const double omega, int max_its = 100,
                          const bool print = true) override;
 
  virtual Method method() const override { return Method::TDHF; }
 
  virtual void clear() override final;
 
  double dV(const DiracSpinor &Fa, const DiracSpinor &Fb, bool conj) const;
 
  virtual double dV(const DiracSpinor &Fa,
                    const DiracSpinor &Fb) const override final;
 
  DiracSpinor dV_rhs(const int kappa_n, const DiracSpinor &Fm,
                     bool conj = false) const;
 
  const std::vector<DiracSpinor> &get_dPsis(const DiracSpinor &Fc,
                                            dPsiType XorY) const;
  const DiracSpinor &get_dPsi_x(const DiracSpinor &Fc, dPsiType XorY,
                                const int kappa_x) const;
 
 
  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;
  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;
 
private:
  void initialise_dPsi();
 
  // Single iteration of TDHF equations
  std::pair<double, std::string> tdhf_core_it(double omega, double eta_damp);
  // Forms set of h*Fc for all core orbitals and all projections
  std::vector<std::vector<DiracSpinor>> form_hFcore() const;
  // 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;
 
public:
  TDHF &operator=(const TDHF &) = delete;
  TDHF(const TDHF &) = default;
  ~TDHF() = default;
};
 
} // namespace ExternalField