Skip to content

Commit

Permalink
Feature: enable energy calculation of uspp
Browse files Browse the repository at this point in the history
  • Loading branch information
@YuLiu98
YuLiu98 committed Oct 24, 2023
1 parent 68b7e17 commit a4718f9
Show file tree
Hide file tree
Showing 9 changed files with 201 additions and 35 deletions.
19 changes: 12 additions & 7 deletions source/module_cell/read_pp_upf201.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -552,15 +552,20 @@ void Pseudopot_upf::read_pseudo_upf201_nonlocal(std::ifstream& ifs)
}

// PP_Q
if (ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_Q>"))
if (ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_Q", true, false))
{
this->qqq.create(nbeta, nbeta);
for (int i = 0; i < nbeta; i++)
ifs.ignore(150, '>');
}
else
{
ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_Q>");
}
this->qqq.create(nbeta, nbeta);
for (int i = 0; i < nbeta; i++)
{
for (int j = 0; j < nbeta; j++)
{
for (int j = 0; j < nbeta; j++)
{
ifs >> qqq(i, j);
}
ifs >> qqq(i, j);
}
}
ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_Q>");
Expand Down
43 changes: 35 additions & 8 deletions source/module_elecstate/elecstate_pw.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,6 @@ ElecStatePW<T, Device>::ElecStatePW(ModulePW::PW_Basis_K* wfc_basis_in,
this->rhopw_smooth = rhopw_in;
this->ppcell = ppcell_in;
this->ucell = ucell_in;
this->vkb = this->ppcell->template get_vkb_data<Real>();
this->init_ks(chg_in, pkv_in, pkv_in->nks, rhodpw_in, bigpw_in);
}

Expand Down Expand Up @@ -246,6 +245,7 @@ void ElecStatePW<T, Device>::add_usrho(const psi::Psi<T, Device>& psi)
const int npwx = psi.get_nbasis() / npol;
const int nbands = psi.get_nbands() * npol;
const int nkb = this->ppcell->nkb;
this->vkb = this->ppcell->template get_vkb_data<Real>();
T* becp = nullptr;
resmem_complex_op()(this->ctx, becp, nbands * nkb, "ElecState<PW>::becp");
Real* becsum = nullptr;
Expand Down Expand Up @@ -385,16 +385,39 @@ void ElecStatePW<T, Device>::add_usrho(const psi::Psi<T, Device>& psi)
}
}
}
delmem_complex_op()(this->ctx, auxk1);
delmem_complex_op()(this->ctx, auxk2);
delmem_complex_op()(this->ctx, aux_gk);
}
}
}
delmem_complex_op()(this->ctx, becp);

// transform soft charge to recip space using smooth grids
T* rhog = nullptr;
resmem_complex_op()(this->ctx, rhog, this->charge->rhopw->npw * GlobalV::NSPIN, "ElecState<PW>::rhog");
setmem_complex_op()(this->ctx, rhog, 0, this->charge->rhopw->npw * GlobalV::NSPIN);
for (int is = 0; is < GlobalV::NSPIN; is++)
{
this->rhopw_smooth->real2recip(this->rho[is], &rhog[is * this->charge->rhopw->npw]);
}

// \sum_lm Q_lm(r) \sum_i <psi_i|beta_l><beta_m|psi_i> w_i
this->addusdens_g(becsum);
// add to the charge density in reciprocal space the part which is due to the US augmentation.
this->addusdens_g(becsum, rhog);

// transform back to real space using dense grids
for (int is = 0; is < GlobalV::NSPIN; is++)
{
this->charge->rhopw->recip2real(&rhog[is * this->charge->rhopw->npw], this->rho[is]);
}

delmem_var_op()(this->ctx, becsum);
delmem_complex_op()(this->ctx, rhog);
}

template <typename T, typename Device>
void ElecStatePW<T, Device>::addusdens_g(const Real* becsum)
void ElecStatePW<T, Device>::addusdens_g(const Real* becsum, T* rhog)
{
const T one{1, 0};
const T zero{0, 0};
Expand All @@ -404,9 +427,6 @@ void ElecStatePW<T, Device>::addusdens_g(const Real* becsum)
Structure_Factor* psf = this->ppcell->psf;
const std::complex<double> ci_tpi = ModuleBase::NEG_IMAG_UNIT * ModuleBase::TWO_PI;

T* aux = nullptr;
resmem_complex_op()(this->ctx, aux, GlobalV::NSPIN * npw, "ElecState<PW>::aux");
setmem_complex_op()(this->ctx, aux, 0, GlobalV::NSPIN * npw);
Real* qmod = nullptr;
resmem_var_op()(this->ctx, qmod, npw, "ElecState<PW>::qmod");
T* qgm = nullptr;
Expand Down Expand Up @@ -478,17 +498,24 @@ void ElecStatePW<T, Device>::addusdens_g(const Real* becsum)
{
for (int jh = ih; jh < atom->ncpp.nh; jh++)
{
// this->ppcell->radial_fft_q(npw, ih, jh, it, qmod, ylmk0, qgm);
this->ppcell->radial_fft_q(this->ctx, npw, ih, jh, it, qmod, ylmk0, qgm);
for (int ig = 0; ig < npw; ig++)
{
aux[is * npw + ig] += qgm[ig] * aux2[ijh * npw + ig];
rhog[is * npw + ig] += qgm[ig] * aux2[ijh * npw + ig];
}
ijh++;
}
}
}
delmem_complex_op()(this->ctx, skk);
delmem_complex_op()(this->ctx, aux2);
delmem_complex_op()(this->ctx, tbecsum);
}
}

delmem_var_op()(this->ctx, qmod);
delmem_complex_op()(this->ctx, qgm);
delmem_var_op()(this->ctx, ylmk0);
}

template class ElecStatePW<std::complex<float>, psi::DEVICE_CPU>;
Expand Down
2 changes: 1 addition & 1 deletion source/module_elecstate/elecstate_pw.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ class ElecStatePW : public ElecState
// add to the charge density in reciprocal space the part which is due to the US augmentation.
void add_usrho(const psi::Psi<T, Device>& psi);
// \sum_lm Q_lm(r) \sum_i <psi_i|beta_l><beta_m|psi_i> w_i
void addusdens_g(const Real* becsum);
void addusdens_g(const Real* becsum, T* rhog);

void init_rho_data();

Expand Down
2 changes: 1 addition & 1 deletion source/module_elecstate/module_charge/charge_mixing.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -628,7 +628,7 @@ void Charge_Mixing::high_freq_mix(std::complex<double>* data,
{
for (int ig = 0; ig < number; ig++)
{
data[ig] += mixing_beta * (data_save[ig] - data[ig]);
data[ig] = data_save[ig] + mixing_beta * (data[ig] - data_save[ig]);
}
}

Expand Down
6 changes: 6 additions & 0 deletions source/module_esolver/esolver_ks_pw.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -715,6 +715,12 @@ void ESolver_KS_PW<T, Device>::updatepot(const int istep, const int iter)
template <typename T, typename Device>
void ESolver_KS_PW<T, Device>::eachiterfinish(const int iter)
{
// liuyu 2023-10-24
// D in uspp need vloc, thus needs update when veff updated
// calculate the effective coefficient matrix for non-local pseudopotential projectors
ModuleBase::matrix veff = this->pelec->pot->get_effective_v();
GlobalC::ppcell.cal_effective_D(veff, this->pw_rhod, GlobalC::ucell);

// print_eigenvalue(GlobalV::ofs_running);
this->pelec->cal_energies(2);
// We output it for restarting the scf.
Expand Down
128 changes: 127 additions & 1 deletion source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -934,7 +934,7 @@ void pseudopot_cell_vnl::radial_fft_q(const int ng,
const int itype,
const double* qnorm,
const ModuleBase::matrix ylm,
std::complex<double>* qg)
std::complex<double>* qg) const
{
// computes the indices which correspond to ih,jh
const int nb = indv(itype, ih);
Expand Down Expand Up @@ -1015,6 +1015,97 @@ void pseudopot_cell_vnl::radial_fft_q(const int ng,
}
}

template <typename FPTYPE, typename Device>
void pseudopot_cell_vnl::radial_fft_q(Device* ctx,
const int ng,
const int ih,
const int jh,
const int itype,
const FPTYPE* qnorm,
const FPTYPE* ylm,
std::complex<FPTYPE>* qg) const
{
using setmem_complex_op = psi::memory::set_memory_op<std::complex<FPTYPE>, Device>;

// computes the indices which correspond to ih,jh
const int nb = indv(itype, ih);
const int mb = indv(itype, jh);
assert(nb < nbetam);
assert(mb < nbetam);
int ijv = 0;
if (nb >= mb)
{
ijv = nb * (nb + 1) / 2 + mb;
}
else
{
ijv = mb * (mb + 1) / 2 + nb;
}
const int ivl = nhtolm(itype, ih);
const int jvl = nhtolm(itype, jh);

setmem_complex_op()(ctx, qg, 0, ng);

const double* qnorm_double = reinterpret_cast<const double*>(qnorm);

// makes the sum over the non zero LM
int l = -1;
std::complex<FPTYPE> pref(0.0, 0.0);
for (int lm = 0; lm < this->lpx(ivl, jvl); lm++)
{
int lp = this->lpl(ivl, jvl, lm);
assert(lp >= 0);
assert(lp < 49);
if (lp == 0)
{
l = 0;
}
else if (lp < 4)
{
l = 1;
}
else if (lp < 9)
{
l = 2;
}
else if (lp < 16)
{
l = 3;
}
else if (lp < 25)
{
l = 4;
}
else if (lp < 36)
{
l = 5;
}
else
{
l = 6;
}
pref = static_cast<std::complex<FPTYPE>>(pow(ModuleBase::NEG_IMAG_UNIT, l) * this->ap(lp, ivl, jvl));

double qm1 = -1.0; // any number smaller than qnorm
double work = 0.0;
for (int ig = 0; ig < ng; ig++)
{
if (std::abs(qnorm_double[ig] - qm1) > 1e-6)
{
work = ModuleBase::PolyInt::Polynomial_Interpolation(this->qrad,
itype,
l,
ijv,
GlobalV::NQXQ,
GlobalV::DQ,
qnorm_double[ig]);
qm1 = qnorm_double[ig];
}
qg[ig] += pref * static_cast<FPTYPE>(work) * ylm[lp * ng + ig];
}
}
}

#ifdef __LCAO
std::complex<double> pseudopot_cell_vnl::Cal_C(int alpha, int lu, int mu, int L, int M) // pengfei Li 2018-3-23
{
Expand Down Expand Up @@ -1713,3 +1804,38 @@ std::complex<double>* pseudopot_cell_vnl::get_qq_so_data() const
int const&,
std::complex<double>*) const;
#endif

template void pseudopot_cell_vnl::radial_fft_q<float, psi::DEVICE_CPU>(psi::DEVICE_CPU*,
const int,
const int,
const int,
const int,
const float*,
const float*,
std::complex<float>*) const;
template void pseudopot_cell_vnl::radial_fft_q<double, psi::DEVICE_CPU>(psi::DEVICE_CPU*,
const int,
const int,
const int,
const int,
const double*,
const double*,
std::complex<double>*) const;
#if defined(__CUDA) || defined(__ROCM)
template void pseudopot_cell_vnl::radial_fft_q<float, psi::DEVICE_GPU>(psi::DEVICE_GPU*,
const int,
const int,
const int,
const int,
const float*,
const float*,
std::complex<float>*) const;
template void pseudopot_cell_vnl::radial_fft_q<double, psi::DEVICE_GPU>(psi::DEVICE_GPU*,
const int,
const int,
const int,
const int,
const double*,
const double*,
std::complex<double>*) const;
#endif
11 changes: 10 additions & 1 deletion source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,7 +141,16 @@ class pseudopot_cell_vnl: public pseudopot_cell_vl
const int itype,
const double* qnorm,
const ModuleBase::matrix ylm,
std::complex<double>* qg);
std::complex<double>* qg) const;
template <typename FPTYPE, typename Device>
void radial_fft_q(Device* ctx,
const int ng,
const int ih,
const int jh,
const int itype,
const FPTYPE* qnorm,
const FPTYPE* ylm,
std::complex<FPTYPE>* qg) const;

/**
* @brief calculate the effective coefficient matrix for non-local pseudopotential projectors
Expand Down
Loading

0 comments on commit a4718f9

Please sign in to comment.