potential_new.h

Go to the documentation of this file.

#ifndef POTENTIALNEW_H
#define POTENTIALNEW_H
 
#include "source_base/complexmatrix.h"
#include "source_hamilt/module_surchem/surchem.h"
#include "source_pw/module_pwdft/VNL_in_pw.h"
#include "source_pw/module_pwdft/VSep_in_pw.h"
#include "source_pw/module_pwdft/structure_factor.h"
#include "pot_base.h"
 
#include <vector>
 
namespace elecstate
{
class Potential : public PotBase
{
  public:
    // default constructor for UT
    Potential(){};
    // In constructor, size of every potential components should be allocated
    // rho_basis_in is the dense grids, rho_basis_smooth_in is the smooth grids in USPP
    // charge density and potential are defined on dense grids,
    // but effective potential needs to be interpolated on smooth grids in order to compute Veff|psi>
    // Note: rho_basis_in and rho_basis_smooth_in are the same in NCPP
    Potential(const ModulePW::PW_Basis* rho_basis_in,
              const ModulePW::PW_Basis* rho_basis_smooth_in,
              const UnitCell* ucell_in,
              const ModuleBase::matrix* vloc_in,
              Structure_Factor* structure_factors_in,
              surchem* solvent_in,
              double* etxc_in,
              double* vtxc_in,
              VSep* vsep_cell_in = nullptr);
    ~Potential();
 
    // initialize potential when SCF begin
    void init_pot(int istep, const Charge*const chg);
    // initialize potential components before SCF
    void pot_register(const std::vector<std::string>& components_list);
    // update potential from current charge
    void update_from_charge(const Charge*const chg, const UnitCell*const ucell);
    // interface for SCF-converged, etxc vtxc for Energy, vnew for force_scc
    void get_vnew(const Charge* chg, ModuleBase::matrix& vnew);
 
    PotBase* get_pot_type(const std::string& pot_type);
 
    // interfaces to get values
    ModuleBase::matrix& get_effective_v()
    {
        return this->v_effective;
    }
    const ModuleBase::matrix& get_effective_v() const
    {
        return this->v_effective;
    }
 
    double* get_effective_v(int is)
    {
        if (this->v_effective.nc > 0)
        {
            return &(this->v_effective(is, 0));
        }
        else
        {
            return nullptr;
        }
    }
    const double* get_effective_v(int is) const
    {
        if (this->v_effective.nc > 0)
        {
            return &(this->v_effective(is, 0));
        }
        else
        {
            return nullptr;
        }
    }
    ModuleBase::matrix& get_effective_vofk()
    {
        return this->vofk_effective;
    }
    const ModuleBase::matrix& get_effective_vofk() const
    {
        return this->vofk_effective;
    }
    double* get_effective_vofk(int is)
    {
        if (this->vofk_effective.nc > 0)
        {
            return &(this->vofk_effective(is, 0));
        }
        else
        {
            return nullptr;
        }
    }
    const double* get_effective_vofk(int is) const
    {
        if (this->vofk_effective.nc > 0)
        {
            return &(this->vofk_effective(is, 0));
        }
        else
        {
            return nullptr;
        }
    }
 
    ModuleBase::matrix& get_veff_smooth()
    {
        return this->veff_smooth;
    }
    const ModuleBase::matrix& get_veff_smooth() const
    {
        return this->veff_smooth;
    }
 
    ModuleBase::matrix& get_vofk_smooth()
    {
        return this->vofk_smooth;
    }
    const ModuleBase::matrix& get_vofk_smooth() const
    {
        return this->vofk_smooth;
    }
 
    template <typename FPTYPE>
    FPTYPE* get_veff_smooth_data();
 
    template <typename FPTYPE>
    FPTYPE* get_vofk_smooth_data();
 
    double* get_fixed_v()
    {
        return this->v_effective_fixed.data();
    }
    const double* get_fixed_v() const
    {
        return this->v_effective_fixed.data();
    }
    const ModulePW::PW_Basis *get_rho_basis() const
    {
        return this->rho_basis_;
    }
    // What about adding a function to get the wfc?
    // This is useful for the calculation of the exx energy
 
 
    double get_vl_of_0() const
    {
        return this->vl_of_0;
    }
 
  private:
    void cal_v_eff(const Charge*const chg, const UnitCell*const ucell, ModuleBase::matrix& v_eff) override;
    void cal_fixed_v(double* vl_pseudo) override;
    // interpolate potential on the smooth mesh if necessary
    void interpolate_vrs();
 
    void allocate();
 
    std::vector<double> v_effective_fixed;
    ModuleBase::matrix v_effective;
 
    ModuleBase::matrix veff_smooth; // used in uspp liuyu 2023-10-12
    ModuleBase::matrix vofk_smooth; // used in uspp liuyu 2023-10-12
 
    ModuleBase::matrix v_xc; // if PAW is used, vxc must be stored separately
 
    float *s_veff_smooth = nullptr, *s_vofk_smooth = nullptr;
    double *d_veff_smooth = nullptr, *d_vofk_smooth = nullptr;
 
    ModuleBase::matrix vofk_effective;
 
    bool fixed_done = false;
 
    // gather etxc and vtxc in Potential, will be used in ESolver
    double* etxc_ = nullptr;
    double* vtxc_ = nullptr;
 
    double vl_of_0 = 0.0;
 
    std::vector<PotBase*> components;
 
    const UnitCell* ucell_ = nullptr;
    const ModuleBase::matrix* vloc_ = nullptr;
    Structure_Factor* structure_factors_ = nullptr;
    surchem* solvent_ = nullptr;
    VSep* vsep_cell = nullptr;
    bool use_gpu_ = false;
};
 
} // namespace elecstate
 
#endif