From 178c24c74f04bc144c9e6503bdd8f793aee63b2b Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Mon, 23 Dec 2024 17:19:06 +0800 Subject: [PATCH 1/3] Feature: enable cal_force and cal_stress in nscf --- source/module_io/read_input_item_system.cpp | 2 +- source/module_relax/relax_driver.cpp | 128 ++++++++++---------- 2 files changed, 63 insertions(+), 67 deletions(-) diff --git a/source/module_io/read_input_item_system.cpp b/source/module_io/read_input_item_system.cpp index 7398dc2e6a..58839c7b09 100644 --- a/source/module_io/read_input_item_system.cpp +++ b/source/module_io/read_input_item_system.cpp @@ -207,7 +207,7 @@ void ReadInput::item_system() item.annotation = "if calculate the force at the end of the electronic iteration"; item.reset_value = [](const Input_Item& item, Parameter& para) { std::vector use_force = {"cell-relax", "relax", "md"}; - std::vector not_use_force = {"get_wf", "get_pchg", "nscf", "get_S"}; + std::vector not_use_force = {"get_wf", "get_pchg", "get_S"}; if (std::find(use_force.begin(), use_force.end(), para.input.calculation) != use_force.end()) { if (!para.input.cal_force) diff --git a/source/module_relax/relax_driver.cpp b/source/module_relax/relax_driver.cpp index 752ce8687e..19fff5cc08 100644 --- a/source/module_relax/relax_driver.cpp +++ b/source/module_relax/relax_driver.cpp @@ -54,55 +54,67 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, UnitCell& uce time_t fstart = time(nullptr); ModuleBase::matrix force; ModuleBase::matrix stress; - if (PARAM.inp.calculation == "scf" || PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax") - { - // I'm considering putting force and stress - // as part of ucell and use ucell to pass information - // back and forth between esolver and relaxation - // but I'll use force and stress explicitly here for now - // calculate the total energy - this->etot = p_esolver->cal_energy(); + // I'm considering putting force and stress + // as part of ucell and use ucell to pass information + // back and forth between esolver and relaxation + // but I'll use force and stress explicitly here for now + + // calculate the total energy + this->etot = p_esolver->cal_energy(); + + // calculate and gather all parts of total ionic forces + if (PARAM.inp.cal_force) + { + p_esolver->cal_force(ucell, force); + } + // calculate and gather all parts of stress + if (PARAM.inp.cal_stress) + { + p_esolver->cal_stress(ucell, stress); + } - // calculate and gather all parts of total ionic forces - if (PARAM.inp.cal_force) + if (PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax") + { + if (PARAM.inp.relax_new) { - p_esolver->cal_force(ucell, force); + stop = rl.relax_step(ucell, force, stress, this->etot); } - // calculate and gather all parts of stress - if (PARAM.inp.cal_stress) + else { - p_esolver->cal_stress(ucell, stress); + stop = rl_old.relax_step(istep, + this->etot, + ucell, + force, + stress, + force_step, + stress_step); // pengfei Li 2018-05-14 } - - if (PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax") + // print structure + // changelog 20240509 + // because I move out the dependence on GlobalV from UnitCell::print_stru_file + // so its parameter is calculated here + bool need_orb = PARAM.inp.basis_type == "pw"; + need_orb = need_orb && PARAM.inp.psi_initializer; + need_orb = need_orb && PARAM.inp.init_wfc.substr(0, 3) == "nao"; + need_orb = need_orb || PARAM.inp.basis_type == "lcao"; + need_orb = need_orb || PARAM.inp.basis_type == "lcao_in_pw"; + std::stringstream ss, ss1; + ss << PARAM.globalv.global_out_dir << "STRU_ION_D"; + ucell.print_stru_file(ss.str(), + PARAM.inp.nspin, + true, + PARAM.inp.calculation == "md", + PARAM.inp.out_mul, + need_orb, + PARAM.globalv.deepks_setorb, + GlobalV::MY_RANK); + + if (Ions_Move_Basic::out_stru) { - if (PARAM.inp.relax_new) - { - stop = rl.relax_step(ucell,force, stress, this->etot); - } - else - { - stop = rl_old.relax_step(istep, - this->etot, - ucell, - force, - stress, - force_step, - stress_step); // pengfei Li 2018-05-14 - } - // print structure - // changelog 20240509 - // because I move out the dependence on GlobalV from UnitCell::print_stru_file - // so its parameter is calculated here - bool need_orb = PARAM.inp.basis_type == "pw"; - need_orb = need_orb && PARAM.inp.psi_initializer; - need_orb = need_orb && PARAM.inp.init_wfc.substr(0, 3) == "nao"; - need_orb = need_orb || PARAM.inp.basis_type == "lcao"; - need_orb = need_orb || PARAM.inp.basis_type == "lcao_in_pw"; - std::stringstream ss, ss1; - ss << PARAM.globalv.global_out_dir << "STRU_ION_D"; - ucell.print_stru_file(ss.str(), + ss1 << PARAM.globalv.global_out_dir << "STRU_ION"; + ss1 << istep << "_D"; + ucell.print_stru_file(ss1.str(), PARAM.inp.nspin, true, PARAM.inp.calculation == "md", @@ -110,34 +122,18 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, UnitCell& uce need_orb, PARAM.globalv.deepks_setorb, GlobalV::MY_RANK); - - if (Ions_Move_Basic::out_stru) - { - ss1 << PARAM.globalv.global_out_dir << "STRU_ION"; - ss1 << istep << "_D"; - ucell.print_stru_file(ss1.str(), - PARAM.inp.nspin, - true, - PARAM.inp.calculation == "md", - PARAM.inp.out_mul, - need_orb, - PARAM.globalv.deepks_setorb, - GlobalV::MY_RANK); - ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU_NOW.cif", - ucell, - "# Generated by ABACUS ModuleIO::CifParser", - "data_?"); - } - - ModuleIO::output_after_relax(stop, p_esolver->conv_esolver, GlobalV::ofs_running); + ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU_NOW.cif", + ucell, + "# Generated by ABACUS ModuleIO::CifParser", + "data_?"); } -#ifdef __RAPIDJSON - // add the energy to outout - Json::add_output_energy(p_esolver->cal_energy() * ModuleBase::Ry_to_eV); -#endif + ModuleIO::output_after_relax(stop, p_esolver->conv_esolver, GlobalV::ofs_running); } + #ifdef __RAPIDJSON + // add the energy to outout + Json::add_output_energy(p_esolver->cal_energy() * ModuleBase::Ry_to_eV); // add Json of cell coo stress force double unit_transform = ModuleBase::RYDBERG_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; double fac = ModuleBase::Ry_to_eV / 0.529177; From 9c46b772934240f30414b8e4af48674e30f3797a Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Mon, 23 Dec 2024 21:15:58 +0800 Subject: [PATCH 2/3] Fix: enable force and stress in uspp nscf --- source/module_elecstate/elecstate_pw.cpp | 8 +++++++- source/module_elecstate/elecstate_pw.h | 5 ++++- source/module_hsolver/hsolver_lcaopw.cpp | 4 ++++ source/module_hsolver/hsolver_pw.cpp | 4 ++++ 4 files changed, 19 insertions(+), 2 deletions(-) diff --git a/source/module_elecstate/elecstate_pw.cpp b/source/module_elecstate/elecstate_pw.cpp index 05fb2b57f3..a22c87bcf7 100644 --- a/source/module_elecstate/elecstate_pw.cpp +++ b/source/module_elecstate/elecstate_pw.cpp @@ -238,7 +238,7 @@ void ElecStatePW::rhoBandK(const psi::Psi& psi) } template -void ElecStatePW::add_usrho(const psi::Psi& psi) +void ElecStatePW::cal_becsum(const psi::Psi& psi) { const T one{1, 0}; const T zero{0, 0}; @@ -392,6 +392,12 @@ void ElecStatePW::add_usrho(const psi::Psi& psi) } } delmem_complex_op()(this->ctx, becp); +} + +template +void ElecStatePW::add_usrho(const psi::Psi& psi) +{ + this->cal_becsum(psi); // transform soft charge to recip space using smooth grids T* rhog = nullptr; diff --git a/source/module_elecstate/elecstate_pw.h b/source/module_elecstate/elecstate_pw.h index d879d872ac..dc992fa37a 100644 --- a/source/module_elecstate/elecstate_pw.h +++ b/source/module_elecstate/elecstate_pw.h @@ -35,6 +35,9 @@ class ElecStatePW : public ElecState virtual void cal_tau(const psi::Psi& psi); + //! calculate becsum for uspp + void cal_becsum(const psi::Psi& psi); + Real* becsum = nullptr; //! init rho_data and kin_r_data @@ -61,7 +64,7 @@ class ElecStatePW : public ElecState //! calcualte rho for each k void rhoBandK(const psi::Psi& psi); - + //! add to the charge density in reciprocal space the part which is due to the US augmentation. void add_usrho(const psi::Psi& psi); diff --git a/source/module_hsolver/hsolver_lcaopw.cpp b/source/module_hsolver/hsolver_lcaopw.cpp index f87a15664a..1f3bb2c17b 100644 --- a/source/module_hsolver/hsolver_lcaopw.cpp +++ b/source/module_hsolver/hsolver_lcaopw.cpp @@ -281,6 +281,10 @@ void HSolverLIP::solve(hamilt::Hamilt* pHamilt, // ESolver_KS_PW::p_hamilt reinterpret_cast*>(pes)->calEBand(); if (skip_charge) { + if (PARAM.globalv.use_uspp) + { + reinterpret_cast*>(pes)->cal_becsum(psi); + } ModuleBase::timer::tick("HSolverLIP", "solve"); return; } diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp index 1a8e7292f6..c217444a1a 100644 --- a/source/module_hsolver/hsolver_pw.cpp +++ b/source/module_hsolver/hsolver_pw.cpp @@ -337,6 +337,10 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, reinterpret_cast*>(pes)->calEBand(); if (skip_charge) { + if (PARAM.globalv.use_uspp) + { + reinterpret_cast*>(pes)->cal_becsum(psi); + } ModuleBase::timer::tick("HSolverPW", "solve"); return; } From 744cc333a9e71e6875bf50edc17cb7798151dc56 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Tue, 24 Dec 2024 14:43:16 +0800 Subject: [PATCH 3/3] update unit tests --- .../test/hsolver_supplementary_mock.h | 42 ++++++++++++++++++- .../module_hsolver/test/test_hsolver_sdft.cpp | 29 ------------- 2 files changed, 41 insertions(+), 30 deletions(-) diff --git a/source/module_hsolver/test/hsolver_supplementary_mock.h b/source/module_hsolver/test/hsolver_supplementary_mock.h index 66e49fb894..436ecd6565 100644 --- a/source/module_hsolver/test/hsolver_supplementary_mock.h +++ b/source/module_hsolver/test/hsolver_supplementary_mock.h @@ -1,5 +1,5 @@ #pragma once -#include "module_elecstate/elecstate.h" +#include "module_elecstate/elecstate_pw.h" #include "module_psi/wavefunc.h" namespace elecstate @@ -62,6 +62,46 @@ void ElecState::init_ks(Charge* chg_in, // pointer for class Charge return; } +template +ElecStatePW::ElecStatePW(ModulePW::PW_Basis_K* wfc_basis_in, + Charge* chg_in, + K_Vectors* pkv_in, + UnitCell* ucell_in, + pseudopot_cell_vnl* ppcell_in, + ModulePW::PW_Basis* rhodpw_in, + ModulePW::PW_Basis* rhopw_in, + ModulePW::PW_Basis_Big* bigpw_in) + : basis(wfc_basis_in) +{ +} + +template +ElecStatePW::~ElecStatePW() +{ +} + +template +void ElecStatePW::psiToRho(const psi::Psi& psi) +{ +} + +template +void ElecStatePW::cal_tau(const psi::Psi& psi) +{ +} + +template +void ElecStatePW::cal_becsum(const psi::Psi& psi) +{ +} + +template class ElecStatePW, base_device::DEVICE_CPU>; +template class ElecStatePW, base_device::DEVICE_CPU>; +#if ((defined __CUDA) || (defined __ROCM)) +template class ElecStatePW, base_device::DEVICE_GPU>; +template class ElecStatePW, base_device::DEVICE_GPU>; +#endif + Potential::~Potential() { } diff --git a/source/module_hsolver/test/test_hsolver_sdft.cpp b/source/module_hsolver/test/test_hsolver_sdft.cpp index 4d3db33096..642ee9daae 100644 --- a/source/module_hsolver/test/test_hsolver_sdft.cpp +++ b/source/module_hsolver/test/test_hsolver_sdft.cpp @@ -23,40 +23,11 @@ Sto_Func::Sto_Func() } template class Sto_Func; - -template <> -elecstate::ElecStatePW, base_device::DEVICE_CPU>::ElecStatePW(ModulePW::PW_Basis_K* wfc_basis_in, - Charge* chg_in, - K_Vectors* pkv_in, - UnitCell* ucell_in, - pseudopot_cell_vnl* ppcell_in, - ModulePW::PW_Basis* rhodpw_in, - ModulePW::PW_Basis* rhopw_in, - ModulePW::PW_Basis_Big* bigpw_in) - : basis(wfc_basis_in) -{ -} - -template<> -elecstate::ElecStatePW, base_device::DEVICE_CPU>::~ElecStatePW() -{ -} - template<> void elecstate::ElecStatePW, base_device::DEVICE_CPU>::init_rho_data() { } -template<> -void elecstate::ElecStatePW, base_device::DEVICE_CPU>::psiToRho(const psi::Psi, base_device::DEVICE_CPU>& psi) -{ -} - -template<> -void elecstate::ElecStatePW, base_device::DEVICE_CPU>::cal_tau(const psi::Psi, base_device::DEVICE_CPU>& psi) -{ -} - template StoChe::StoChe(const int& nche, const int& method, const REAL& emax_sto, const REAL& emin_sto) {