CorePolarisation.hpp

#pragma once
#include "DiracOperator/TensorOperator.hpp"
#include "qip/String.hpp"
#include <string>
#include <vector>
class DiracSpinor;
namespace DiracOperator {
class TensorOperator;
}
 
namespace ExternalField {
 
enum class Method { TDHF, basis, diagram, none, Error };
 
inline Method ParseMethod(std::string_view str) {
  return qip::ci_compare(str, "TDHF")       ? Method::TDHF :
         qip::ci_compare(str, "true")       ? Method::TDHF :
         qip::ci_compare(str, "basis")      ? Method::basis :
         qip::ci_compare(str, "tdhf_basis") ? Method::basis :
         qip::ci_compare(str, "tdhfbasis")  ? Method::basis :
         qip::ci_compare(str, "diagram")    ? Method::diagram :
         qip::ci_compare(str, "rpad")       ? Method::diagram :
         qip::ci_compare(str, "rpa(d)")     ? Method::diagram :
         qip::ci_compare(str, "none")       ? Method::none :
         qip::ci_compare(str, "false")      ? Method::none :
         qip::ci_compare(str, "")           ? Method::none :
                                              Method::Error;
}
 
enum class dPsiType { X, Y };
enum class StateType { bra, ket }; // lhs, rhs
 
class CorePolarisation {
 
protected:
  CorePolarisation(const DiracOperator::TensorOperator *const h)
      : m_h(h), m_rank(h->rank()), m_pi(h->parity()), m_imag(h->imaginaryQ()) {}
 
protected:
  const DiracOperator::TensorOperator *m_h;
  double m_core_eps{1.0};
  int m_core_its{0};
  double m_core_omega{0.0};
  int m_rank;
  int m_pi;
  bool m_imag;
 
  double m_eta{0.4};
  double m_eps{1.0e-10};
 
public:
  double last_eps() const { return m_core_eps; }
  double last_its() const { return m_core_its; }
  double last_omega() const { return m_core_omega; }
  int rank() const { return m_rank; }
  int parity() const { return m_pi; }
  bool imagQ() const { return m_imag; }
 
  double &eps_target() { return m_eps; }
  double eps_target() const { return m_eps; }
 
  double eta() const { return m_eta; }
  void set_eta(double eta) {
    assert(eta >= 0.0 && eta < 1 && "Must have 0 <= eta < 1");
    m_eta = eta;
  }
 
  virtual Method method() const = 0;
 
  virtual void solve_core(const double omega, int max_its = 100,
                          const bool print = true) = 0;
 
  virtual void clear() = 0;
 
  virtual double dV(const DiracSpinor &Fn, const DiracSpinor &Fm) const = 0;
 
public:
  CorePolarisation &operator=(const CorePolarisation &) = delete;
  CorePolarisation(const CorePolarisation &) = default;
  virtual ~CorePolarisation() = default;
};
 
} // namespace ExternalField