diff --git a/source/module_cell/read_pp_upf201.cpp b/source/module_cell/read_pp_upf201.cpp index 30e4d47112..6cf0021142 100644 --- a/source/module_cell/read_pp_upf201.cpp +++ b/source/module_cell/read_pp_upf201.cpp @@ -552,15 +552,20 @@ void Pseudopot_upf::read_pseudo_upf201_nonlocal(std::ifstream& ifs) } // PP_Q - if (ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "")) + if (ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "qqq.create(nbeta, nbeta); - for (int i = 0; i < nbeta; i++) + ifs.ignore(150, '>'); + } + else + { + ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, ""); + } + 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, ""); diff --git a/source/module_elecstate/elecstate_pw.cpp b/source/module_elecstate/elecstate_pw.cpp index cfd00c7255..bdc6df14ed 100644 --- a/source/module_elecstate/elecstate_pw.cpp +++ b/source/module_elecstate/elecstate_pw.cpp @@ -25,7 +25,6 @@ ElecStatePW::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(); this->init_ks(chg_in, pkv_in, pkv_in->nks, rhodpw_in, bigpw_in); } @@ -246,6 +245,7 @@ void ElecStatePW::add_usrho(const psi::Psi& 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(); T* becp = nullptr; resmem_complex_op()(this->ctx, becp, nbands * nkb, "ElecState::becp"); Real* becsum = nullptr; @@ -385,16 +385,39 @@ void ElecStatePW::add_usrho(const psi::Psi& 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::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 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 -void ElecStatePW::addusdens_g(const Real* becsum) +void ElecStatePW::addusdens_g(const Real* becsum, T* rhog) { const T one{1, 0}; const T zero{0, 0}; @@ -404,9 +427,6 @@ void ElecStatePW::addusdens_g(const Real* becsum) Structure_Factor* psf = this->ppcell->psf; const std::complex ci_tpi = ModuleBase::NEG_IMAG_UNIT * ModuleBase::TWO_PI; - T* aux = nullptr; - resmem_complex_op()(this->ctx, aux, GlobalV::NSPIN * npw, "ElecState::aux"); - setmem_complex_op()(this->ctx, aux, 0, GlobalV::NSPIN * npw); Real* qmod = nullptr; resmem_var_op()(this->ctx, qmod, npw, "ElecState::qmod"); T* qgm = nullptr; @@ -478,17 +498,24 @@ void ElecStatePW::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, psi::DEVICE_CPU>; diff --git a/source/module_elecstate/elecstate_pw.h b/source/module_elecstate/elecstate_pw.h index 1ae34c3c54..c84b19975b 100644 --- a/source/module_elecstate/elecstate_pw.h +++ b/source/module_elecstate/elecstate_pw.h @@ -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& psi); // \sum_lm Q_lm(r) \sum_i w_i - void addusdens_g(const Real* becsum); + void addusdens_g(const Real* becsum, T* rhog); void init_rho_data(); diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 1ec363097a..edf338de56 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -628,7 +628,7 @@ void Charge_Mixing::high_freq_mix(std::complex* 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]); } } diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 532f08d5da..60c2182374 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -715,6 +715,12 @@ void ESolver_KS_PW::updatepot(const int istep, const int iter) template void ESolver_KS_PW::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. diff --git a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp index b62ba655d4..471deb8884 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp @@ -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* qg) + std::complex* qg) const { // computes the indices which correspond to ih,jh const int nb = indv(itype, ih); @@ -1015,6 +1015,97 @@ void pseudopot_cell_vnl::radial_fft_q(const int ng, } } +template +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* qg) const +{ + using setmem_complex_op = psi::memory::set_memory_op, 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(qnorm); + + // makes the sum over the non zero LM + int l = -1; + std::complex 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>(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(work) * ylm[lp * ng + ig]; + } + } +} + #ifdef __LCAO std::complex pseudopot_cell_vnl::Cal_C(int alpha, int lu, int mu, int L, int M) // pengfei Li 2018-3-23 { @@ -1713,3 +1804,38 @@ std::complex* pseudopot_cell_vnl::get_qq_so_data() const int const&, std::complex*) const; #endif + + template void pseudopot_cell_vnl::radial_fft_q(psi::DEVICE_CPU*, + const int, + const int, + const int, + const int, + const float*, + const float*, + std::complex*) const; + template void pseudopot_cell_vnl::radial_fft_q(psi::DEVICE_CPU*, + const int, + const int, + const int, + const int, + const double*, + const double*, + std::complex*) const; +#if defined(__CUDA) || defined(__ROCM) + template void pseudopot_cell_vnl::radial_fft_q(psi::DEVICE_GPU*, + const int, + const int, + const int, + const int, + const float*, + const float*, + std::complex*) const; + template void pseudopot_cell_vnl::radial_fft_q(psi::DEVICE_GPU*, + const int, + const int, + const int, + const int, + const double*, + const double*, + std::complex*) const; +#endif diff --git a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h index f3399cc897..023e9acac2 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h +++ b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h @@ -141,7 +141,16 @@ class pseudopot_cell_vnl: public pseudopot_cell_vl const int itype, const double* qnorm, const ModuleBase::matrix ylm, - std::complex* qg); + std::complex* qg) const; + template + 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* qg) const; /** * @brief calculate the effective coefficient matrix for non-local pseudopotential projectors diff --git a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp index 99bcd83b52..cc538556f6 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp @@ -134,7 +134,6 @@ HamiltPW::HamiltPW(const HamiltPW *hamilt) this->qq_nt = hamilt->qq_nt; this->qq_so = hamilt->qq_so; this->vkb = hamilt->vkb; - this->becp = hamilt->becp; OperatorPW, Device_in> * node = reinterpret_cast, Device_in> *>(hamilt->ops); @@ -208,10 +207,12 @@ void HamiltPW::sPsi(const T* psi_in, // psi syncmem_op()(this->ctx, this->ctx, spsi, psi_in, static_cast(nbands * nrow)); if (GlobalV::use_uspp) { + T* becp = nullptr; + T* ps = nullptr; // psi updated, thus update if (this->ppcell->nkb > 0) { - resmem_complex_op()(this->ctx, this->becp, nbands * this->ppcell->nkb, "Hamilt::becp"); + resmem_complex_op()(this->ctx, becp, nbands * this->ppcell->nkb, "Hamilt::becp"); char transa = 'C'; char transb = 'N'; if (nbands == 1) @@ -227,7 +228,7 @@ void HamiltPW::sPsi(const T* psi_in, // psi psi_in, inc, &this->zero, - this->becp, + becp, inc); } else @@ -244,18 +245,17 @@ void HamiltPW::sPsi(const T* psi_in, // psi psi_in, nrow, &this->zero, - this->becp, + becp, this->ppcell->nkb); } - Parallel_Reduce::reduce_pool(this->becp, this->ppcell->nkb * nbands); + Parallel_Reduce::reduce_pool(becp, this->ppcell->nkb * nbands); } resmem_complex_op()(this->ctx, ps, this->ppcell->nkb * nbands, "Hamilt::ps"); setmem_complex_op()(this->ctx, ps, 0, this->ppcell->nkb * nbands); // spsi = psi + sum qq |beta> - int iat = 0; if (GlobalV::NONCOLIN) { // spsi_nc @@ -286,6 +286,7 @@ void HamiltPW::sPsi(const T* psi_in, // psi } for (int ia = 0; ia < atoms->na; ia++) { + const int iat = GlobalC::ucell.itia2iat(it, ia); gemm_op()(this->ctx, transa, transb, @@ -300,14 +301,9 @@ void HamiltPW::sPsi(const T* psi_in, // psi &this->zero, &ps[this->ppcell->indv_ijkb0[iat]], this->ppcell->nkb); - iat++; } delmem_complex_op()(ctx, qqc); } - else - { - iat += atoms->na; - } } if (nbands == 1) @@ -344,11 +340,10 @@ void HamiltPW::sPsi(const T* psi_in, // psi nrow); } } + delmem_complex_op()(this->ctx, ps); + delmem_complex_op()(this->ctx, becp); } - delmem_complex_op()(ctx, ps); - delmem_complex_op()(ctx, becp); - ModuleBase::TITLE("HamiltPW", "sPsi"); } diff --git a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.h b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.h index 6320081754..ff7fafd2b6 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.h +++ b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.h @@ -39,8 +39,6 @@ class HamiltPW : public Hamilt // used in sPhi, which are calculated in hPsi or sPhi const pseudopot_cell_vnl* ppcell = nullptr; mutable T* vkb = nullptr; - mutable T* becp = nullptr; - mutable T* ps = nullptr; Real* qq_nt = nullptr; T* qq_so = nullptr; const T one{1, 0};