BSplineBasis.hpp

#pragma once
#include "IO/InputBlock.hpp"
#include "LinAlg/include.hpp"
#include <memory>
#include <string>
#include <utility>
class DiracSpinor;
class Wavefunction;
class Grid;
namespace MBPT {
class CorrelationPotential;
}
namespace IO {
class InputBlock;
}
 
namespace SplineBasis {
 
enum class SplineType { Derevianko, Johnson };
inline auto parseSplineType(std::string_view type) {
  return (type == "Johnson" || type == "johnson") ? SplineType::Johnson :
                                                    SplineType::Derevianko;
}
 
struct Parameters {
  Parameters() {}
  Parameters(const std::string &states, std::size_t n, std::size_t k, double r0,
             double reps, double rmax, const std::string &positron = "",
             SplineType itype = SplineType::Derevianko,
             bool in_orthogonalise = false, bool in_verbose = true);
  Parameters(IO::InputBlock input);
 
  std::string states{};
  std::size_t n{}, k{};
  double r0{}, reps{}, rmax{};
  std::string positron{};
  SplineType type{SplineType::Derevianko};
  bool orthogonalise{false};
  bool verbose{true};
};
 
 
std::vector<DiracSpinor> form_basis(const Parameters &params,
                                    const Wavefunction &wf,
                                    const bool correlationsQ = false);
 
double check(const std::vector<DiracSpinor> &basis,
             const std::vector<DiracSpinor> &orbs, bool print_warning = true);
 
std::pair<std::vector<DiracSpinor>, std::vector<DiracSpinor>>
form_spline_basis(const int kappa, const std::size_t n_states,
                  const std::size_t k_spl, const double r0_spl,
                  const double rmax_spl, std::shared_ptr<const Grid> rgrid,
                  const double alpha,
                  SplineType itype = SplineType::Derevianko);
 
std::pair<LinAlg::Matrix<double>, LinAlg::Matrix<double>>
fill_Hamiltonian_matrix(const std::vector<DiracSpinor> &spl_basis,
                        const std::vector<DiracSpinor> &d_basis,
                        const Wavefunction &wf,
                        const bool correlationsQ = false,
                        SplineType itype = SplineType::Derevianko);
 
void add_NotreDameBoundary(LinAlg::Matrix<double> *Aij, const int kappa,
                           const double alpha);
 
void expand_basis_orbitals(std::vector<DiracSpinor> *basis,
                           std::vector<DiracSpinor> *basis_positron,
                           const std::vector<DiracSpinor> &spl_basis,
                           const int kappa, const int max_n, int max_n_positron,
                           const LinAlg::Vector<double> &e_values,
                           const LinAlg::Matrix<double> &e_vectors,
                           const Wavefunction &wf);
 
std::vector<double> sumrule_TKR(const std::vector<DiracSpinor> &basis,
                                const std::vector<double> &r,
                                bool print = false);
 
std::vector<double> sumrule_DG(int nDG, const std::vector<DiracSpinor> &basis,
                               const Grid &gr, double alpha, bool print);
 
std::pair<double, double> r_completeness(const DiracSpinor &Fv,
                                         const std::vector<DiracSpinor> &basis,
                                         const Grid &gr, bool print = false);
} // namespace SplineBasis