numerical_radial.h

Go to the documentation of this file.

#ifndef NUMERICAL_RADIAL_H_
#define NUMERICAL_RADIAL_H_
 
#include <cassert>
#include <string>
#include <memory>
 
#include "source_base/spherical_bessel_transformer.h"
#include "source_basis/module_ao/ORB_atomic_lm.h"
 
class NumericalRadial
{
public:
    NumericalRadial() = default;
    NumericalRadial(NumericalRadial const&); 
 
    NumericalRadial& operator=(NumericalRadial const&);
 
    ~NumericalRadial();
 
    void build(const int l,
               const bool for_r_space,
               const int ngrid,
               const double* const grid,
               const double* const value,
               const int p = 0,
               const int izeta = 0,
               const std::string symbol = "",
               const int itype = 0,
               const bool init_sbt = true
    );
 
    void to_numerical_orbital_lm(Numerical_Orbital_Lm& orbital_lm,
                                 const int nk_legacy = 4005, // equivalent to lcao_ecut = 1600
                                 const double lcao_dk = 0.01,
                                 const bool out_element_info = false,
                                 const bool cal_force = false) const;
 
    void set_transformer(ModuleBase::SphericalBesselTransformer sbt, int update = 0);
 
    void set_grid(const bool for_r_space, const int ngrid, const double* const grid, const char mode = 'i');
 
    void set_uniform_grid(const bool for_r_space,
                          const int ngrid,
                          const double cutoff,
                          const char mode = 'i',
                          const bool enable_fft = false);
 
    void set_value(const bool for_r_space,
                   const double* const value,
                   const int p
    );
 
    void wipe(const bool r_space = true, const bool k_space = true);
 
 
    void radtab(const char op,              
                const NumericalRadial& ket, 
                const int l,                
                double* const table,        
                const int nr_tab,           
                const double rmax_tab,      
                const bool deriv = false    
    ) const;
 
    void normalize(bool for_r_space = true);
 
    std::string const& symbol() const { return symbol_; }
    int itype() const { return itype_; }
    int izeta() const { return izeta_; }
    int l() const { return l_; }
    int nr() const { return nr_; } // paired with rmax(), not rcut!
    int nk() const { return nk_; }
    double rcut() const { return rgrid_[std::min(ircut_, nr_-1)]; } 
    double kcut() const { return kgrid_[std::min(ikcut_, nk_-1)]; }
    double rmax() const { return rgrid_[nr_-1]; } 
    double kmax() const { return kgrid_[nk_-1]; }
    const double* rgrid() const { return rgrid_; }
    const double* kgrid() const { return kgrid_; }
    const double* rvalue() const { return rvalue_; }
    const double* kvalue() const { return kvalue_; }
    double pr() const { return pr_; }
    double pk() const { return pk_; }
    bool is_fft_compliant() const { return is_fft_compliant_; }
    ModuleBase::SphericalBesselTransformer sbt() const { return sbt_; }
 
    double rgrid(int ir) const { return rgrid_[ir]; }
    double kgrid(int ik) const { return kgrid_[ik]; }
    double rvalue(int ir) const { return rvalue_[ir]; }
    double kvalue(int ik) const { return kvalue_[ik]; }
 
private:
    std::string symbol_ = "";   
    int itype_ = 0;             
    int l_ = -1;                
    int izeta_ = 0;             
 
    int nr_ = 0;                
    int nk_ = 0;                
 
    double* rgrid_ = nullptr;   
    double* kgrid_ = nullptr;   
 
    int ircut_ = 0;
    int ikcut_ = 0;
 
    double* rvalue_ = nullptr;  
    double* kvalue_ = nullptr;  
 
    bool is_fft_compliant_ = false;
 
 
    int pr_ = 0; 
 
    int pk_ = 0; 
 
    ModuleBase::SphericalBesselTransformer sbt_;
 
    void transform(const bool forward);
 
    void set_icut(const bool for_r_space, const bool for_k_space, const double tol = 1e-15);
 
    // FIXME is_uniform and is_fft_compliant should be more robust for arrays whose elements
    // are all close to machine precision
 
    static bool is_uniform(const int n, const double* const grid, const double tol = 1e-15);
 
    static bool is_fft_compliant(const int nr,
                                 const double* const rgrid,
                                 const int nk,
                                 const double* const kgrid,
                                 const double tol = 1e-15
                                 );
};
 
#endif