From 7056823227c43fc3772238559ec2dc8b965d3d6a Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Wed, 24 Jan 2024 15:47:16 +0800 Subject: [PATCH 01/11] add a new input parameter mixing_dmr and corresponding UnitTests --- source/module_base/global_variable.cpp | 1 + source/module_base/global_variable.h | 1 + source/module_io/input.cpp | 6 ++++++ source/module_io/input.h | 1 + source/module_io/input_conv.cpp | 1 + source/module_io/parameter_pool.cpp | 4 ++++ source/module_io/test/input_conv_test.cpp | 2 ++ source/module_io/test/input_test_para.cpp | 1 + source/module_io/test/write_input_test.cpp | 1 + source/module_io/write_input.cpp | 1 + 10 files changed, 19 insertions(+) diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp index eb0dc636e2..e76f45d7c0 100644 --- a/source/module_base/global_variable.cpp +++ b/source/module_base/global_variable.cpp @@ -255,6 +255,7 @@ double MIXING_GG0_MAG = 1.00; double MIXING_GG0_MIN = 0.1; double MIXING_ANGLE = 0.0; bool MIXING_TAU = 0; +bool MIXING_DMR = 0; //========================================================== // device flags added by denghui diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h index 9fa7dc1c8a..1bbe1edb91 100644 --- a/source/module_base/global_variable.h +++ b/source/module_base/global_variable.h @@ -284,6 +284,7 @@ extern double MIXING_BETA_MAG; extern double MIXING_GG0_MAG; extern double MIXING_GG0_MIN; extern double MIXING_ANGLE; +extern bool MIXING_DMR; //========================================================== // device flags added by denghui diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index 2f6b304a1c..32a5a358e3 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -311,6 +311,7 @@ void Input::Default(void) mixing_angle = -10.0; // defaultly close for npsin = 4 mixing_tau = false; mixing_dftu = false; + mixing_dmr = false; // whether to mixing real space density matrix //---------------------------------------------------------- // potential / charge / wavefunction / energy //---------------------------------------------------------- @@ -1291,6 +1292,10 @@ bool Input::Read(const std::string& fn) { read_bool(ifs, mixing_dftu); } + else if (strcmp("mixing_dmr", word) == 0) + { + read_bool(ifs, mixing_dmr); + } //---------------------------------------------------------- // charge / potential / wavefunction //---------------------------------------------------------- @@ -3338,6 +3343,7 @@ void Input::Bcast() Parallel_Common::bcast_double(mixing_angle); Parallel_Common::bcast_bool(mixing_tau); Parallel_Common::bcast_bool(mixing_dftu); + Parallel_Common::bcast_bool(mixing_dmr); Parallel_Common::bcast_string(read_file_dir); Parallel_Common::bcast_string(init_wfc); diff --git a/source/module_io/input.h b/source/module_io/input.h index ce723bab30..f732aed897 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -241,6 +241,7 @@ class Input bool mixing_tau; // whether to mix tau in mgga bool mixing_dftu; //whether to mix locale in DFT+U + bool mixing_dmr; // whether to mix real space density matrix //========================================================== // potential / charge / wavefunction / energy diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 8e00b1f774..f8d10df8e8 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -757,6 +757,7 @@ void Input_Conv::Convert(void) GlobalV::MIXING_GG0_MIN = INPUT.mixing_gg0_min; GlobalV::MIXING_ANGLE = INPUT.mixing_angle; GlobalV::MIXING_TAU = INPUT.mixing_tau; + GlobalV::MIXING_DMR = INPUT.mixing_dmr; //----------------------------------------------- // Quasiatomic Orbital analysis diff --git a/source/module_io/parameter_pool.cpp b/source/module_io/parameter_pool.cpp index 524df9de87..1763323054 100644 --- a/source/module_io/parameter_pool.cpp +++ b/source/module_io/parameter_pool.cpp @@ -850,6 +850,10 @@ void input_parameters_set(std::map input_parameters { INPUT.mixing_dftu = *static_cast(input_parameters["mixing_dftu"].get()); } + else if (input_parameters.count("mixing_dmr") != 0) + { + INPUT.mixing_dmr = *static_cast(input_parameters["mixing_dmr"].get()); + } else if (input_parameters.count("init_wfc") != 0) { INPUT.init_wfc = static_cast(input_parameters["init_wfc"].get())->c_str(); diff --git a/source/module_io/test/input_conv_test.cpp b/source/module_io/test/input_conv_test.cpp index a566827792..5b9c93dd77 100644 --- a/source/module_io/test/input_conv_test.cpp +++ b/source/module_io/test/input_conv_test.cpp @@ -181,9 +181,11 @@ TEST_F(InputConvTest, Conv) EXPECT_EQ(GlobalV::of_kernel_file,"WTkernel.txt"); EXPECT_EQ(GlobalV::global_readin_dir,GlobalV::global_out_dir); EXPECT_EQ(GlobalV::sc_mag_switch,0); + EXPECT_TRUE(GlobalV::decay_grad_switch); EXPECT_EQ(GlobalV::sc_file, "sc.json"); EXPECT_EQ(GlobalV::MIXING_RESTART,0); + EXPECT_EQ(GlobalV::MIXING_DMR,false); } TEST_F(InputConvTest, ConvRelax) diff --git a/source/module_io/test/input_test_para.cpp b/source/module_io/test/input_test_para.cpp index d005fdfccc..f74544ecb6 100644 --- a/source/module_io/test/input_test_para.cpp +++ b/source/module_io/test/input_test_para.cpp @@ -382,6 +382,7 @@ TEST_F(InputParaTest, Bcast) EXPECT_FALSE(INPUT.mixing_tau); EXPECT_FALSE(INPUT.mixing_dftu); EXPECT_EQ(INPUT.mixing_restart,0); + EXPECT_EQ(INPUT.mixing_dmr,false); EXPECT_EQ(INPUT.out_bandgap, 0); EXPECT_EQ(INPUT.out_mat_t, 0); diff --git a/source/module_io/test/write_input_test.cpp b/source/module_io/test/write_input_test.cpp index 8dccb5627a..0103c0a14b 100644 --- a/source/module_io/test/write_input_test.cpp +++ b/source/module_io/test/write_input_test.cpp @@ -394,6 +394,7 @@ TEST_F(write_input, Mixing7) EXPECT_THAT(output, testing::HasSubstr("mixing_tau 0 #whether to mix tau in mGGA calculation")); EXPECT_THAT(output, testing::HasSubstr("mixing_dftu 0 #whether to mix locale in DFT+U calculation")); EXPECT_THAT(output, testing::HasSubstr("mixing_restart 0 #which step to restart mixing during SCF")); + EXPECT_THAT(output, testing::HasSubstr("mixing_dmr 0 #whether to mix real-space density matrix")); EXPECT_THAT(output, testing::HasSubstr("")); ifs.close(); remove("write_input_test.log"); diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp index e4c2c68464..dd93cbd35b 100644 --- a/source/module_io/write_input.cpp +++ b/source/module_io/write_input.cpp @@ -256,6 +256,7 @@ ModuleBase::GlobalFunc::OUTP(ofs, "out_bandgap", out_bandgap, "if true, print ou ModuleBase::GlobalFunc::OUTP(ofs, "mixing_angle", mixing_angle, "angle mixing parameter for non-colinear calculations"); ModuleBase::GlobalFunc::OUTP(ofs, "mixing_tau", mixing_tau, "whether to mix tau in mGGA calculation"); ModuleBase::GlobalFunc::OUTP(ofs, "mixing_dftu", mixing_dftu, "whether to mix locale in DFT+U calculation"); + ModuleBase::GlobalFunc::OUTP(ofs, "mixing_dmr", mixing_dmr, "whether to mix real-space density matrix"); ofs << "\n#Parameters (8.DOS)" << std::endl; ModuleBase::GlobalFunc::OUTP(ofs, "dos_emin_ev", dos_emin_ev, "minimal range for dos"); From 69be0919291e7cb19606863c88a0f3113becf1ed Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Wed, 24 Jan 2024 16:10:04 +0800 Subject: [PATCH 02/11] delete Kerker_screen_real_test() --- source/module_elecstate/module_charge/charge_mixing.h | 1 - 1 file changed, 1 deletion(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 993c229057..32505177e9 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -56,7 +56,6 @@ class Charge_Mixing * */ void Kerker_screen_real(double* rho); - void Kerker_screen_real_test(double* rho); /** * @brief Inner product of two complex vectors From c846caf13d2336ff630da1bc81ac838d93ea39f8 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Wed, 24 Jan 2024 16:33:10 +0800 Subject: [PATCH 03/11] delete autoset() in charge_mixing.h --- .../module_elecstate/module_charge/charge_mixing.h | 14 +------------- 1 file changed, 1 insertion(+), 13 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 32505177e9..d06313edea 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -86,19 +86,7 @@ class Charge_Mixing const int& mixing_ndim_in, const double& mixing_gg0_in, const bool& mixing_tau_in, - const double& mixing_beta_mag_in); // mohan add mixing_gg0_in 2014-09-27 - - // /** - // * @brief use auto set - // * - // */ - // void need_auto_set(); - - // /** - // * @brief auto set mixing gg0 and mixing_beta - // * - // */ - // void auto_set(const double& bandgap_in, const UnitCell& ucell_); + const double& mixing_beta_mag_in); /** * @brief Get the drho From 1b1c19e4271c0301b480d38eb0cecb3d1b4070c8 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Wed, 24 Jan 2024 16:47:37 +0800 Subject: [PATCH 04/11] declare a new pubilc member dmr_mdata, whose memory can be allocated by allocate_mixing_dmr() --- .../module_charge/charge_mixing.cpp | 20 +++++++++++++++++++ .../module_charge/charge_mixing.h | 7 +++++++ 2 files changed, 27 insertions(+) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 227d9a050a..d76abdb9df 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -114,6 +114,26 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, return; } +void Charge_Mixing::allocate_mixing_dmr(int nnr) +{ + ModuleBase::TITLE("Charge_Mixing", "allocate_mixing_dmr"); + ModuleBase::timer::tick("Charge_Mixing", "allocate_mixing_dmr"); + // + const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1; + // allocate memory for dmr_mdata + if (GlobalV::SCF_THR_TYPE == 1) + { + ModuleBase::WARNING_QUIT("Charge_Mixing", "This Mixing of Density Matrix is not supported for PW basis yet"); + } + else if (GlobalV::SCF_THR_TYPE == 2) + { + this->mixing->init_mixing_data(this->dmr_mdata, nnr * dmr_nspin, sizeof(double)); + } + ModuleBase::timer::tick("Charge_Mixing", "allocate_mixing_dmr"); + + return; +} + void Charge_Mixing::set_rhopw(ModulePW::PW_Basis* rhopw_in, ModulePW::PW_Basis* rhodpw_in) { this->rhopw = rhopw_in; diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index d06313edea..b91b5698ce 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -16,6 +16,7 @@ class Charge_Mixing Base_Mixing::Mixing_Data rho_mdata; ///< Mixing data for charge density Base_Mixing::Mixing_Data tau_mdata; ///< Mixing data for kinetic energy density Base_Mixing::Mixing_Data nhat_mdata; ///< Mixing data for compensation density + Base_Mixing::Mixing_Data dmr_mdata; ///< Mixing data for real space density matrix Base_Mixing::Plain_Mixing* mixing_highf = nullptr; ///< The high_frequency part is mixed by plain mixing method. @@ -88,6 +89,12 @@ class Charge_Mixing const bool& mixing_tau_in, const double& mixing_beta_mag_in); + /** + * @brief allocate memory of dmr_mdata + * + */ + void allocate_mixing_dmr(int nnr); + /** * @brief Get the drho * From 21997c8f9017dbed7f530465f3ba2e0152ec5410 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Thu, 25 Jan 2024 13:18:57 +0800 Subject: [PATCH 05/11] allocate dmr_mdata in ESolver_KS_LCAO::eachiterinit() --- source/module_elecstate/module_charge/charge_mixing.cpp | 3 +++ source/module_esolver/esolver_ks_lcao.cpp | 8 ++++++++ 2 files changed, 11 insertions(+) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index d76abdb9df..f6ec11a4a5 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -116,6 +116,9 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, void Charge_Mixing::allocate_mixing_dmr(int nnr) { + // Note that: we cannot allocate memory for dmr_mdata in set_mixing. + // since the size of dmr_mdata is given by the size of HContainer.nnr, which is calculated in DensityMatrix::init_DMR(). + // and DensityMatrix::init_DMR() is called in beforescf(). While set_mixing() is called in ESolver_KS::Init(). ModuleBase::TITLE("Charge_Mixing", "allocate_mixing_dmr"); ModuleBase::timer::tick("Charge_Mixing", "allocate_mixing_dmr"); // diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 4c4d6c342f..b350014051 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -501,6 +501,14 @@ namespace ModuleESolver GlobalV::MIXING_GG0, GlobalV::MIXING_TAU, GlobalV::MIXING_BETA_MAG); + // allocate memory for dmr_mdata + if (GlobalV::MIXING_DMR) + { + const elecstate::DensityMatrix* dm + = dynamic_cast*>(this->pelec)->get_DM(); + int nnr_tmp = dm->get_DMR_pointer(1)->get_nnr(); + this->p_chgmix->allocate_mixing_dmr(nnr_tmp); + } } this->p_chgmix->mix_reset(); } From 7a30ec5d44f2bb90a8c019d27f052245adf6ae02 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Thu, 25 Jan 2024 14:59:09 +0800 Subject: [PATCH 06/11] define a function save_DMR in DensityMatrix, and save perious DMR in ESolver_KS_LCAO::hamilt2density() --- .../module_dm/density_matrix.cpp | 32 +++++++++++++++++++ .../module_dm/density_matrix.h | 6 ++++ source/module_esolver/esolver_ks_lcao.cpp | 7 ++++ 3 files changed, 45 insertions(+) diff --git a/source/module_elecstate/module_dm/density_matrix.cpp b/source/module_elecstate/module_dm/density_matrix.cpp index c284099135..b7e8198bed 100644 --- a/source/module_elecstate/module_dm/density_matrix.cpp +++ b/source/module_elecstate/module_dm/density_matrix.cpp @@ -382,6 +382,38 @@ int DensityMatrix::get_DMK_ncol() const return this->_paraV->ncol; } +template +void DensityMatrix::save_DMR() +{ + ModuleBase::TITLE("DensityMatrix", "save_DMR"); + ModuleBase::timer::tick("DensityMatrix", "save_DMR"); + + const int nnr = this->_DMR[0]->get_nnr(); + // allocate if _DMR_save is empty + if(_DMR_save.size() == 0) + { + _DMR_save.resize(this->_DMR.size()); + for(int is=0;is_DMR.size();is++) + { + _DMR_save[is].resize(nnr); + } + } + // save _DMR to _DMR_save + for(int is=0;is_DMR.size();is++) + { + TR* DMR_pointer = this->_DMR[is]->get_wrapper(); + TR* DMR_save_pointer = _DMR_save[is].data(); + // set to zero + ModuleBase::GlobalFunc::ZEROS(DMR_save_pointer, nnr); + for(int i=0;i void DensityMatrix::cal_DMR_test() diff --git a/source/module_elecstate/module_dm/density_matrix.h b/source/module_elecstate/module_dm/density_matrix.h index 69c3ce9a7c..b3aaeab186 100644 --- a/source/module_elecstate/module_dm/density_matrix.h +++ b/source/module_elecstate/module_dm/density_matrix.h @@ -190,6 +190,11 @@ namespace elecstate * @param ik k-point index */ void read_DMK(const std::string directory, const int ispin, const int ik); + + /** + * @brief save _DMR into _DMR_save + */ + void save_DMR(); std::vector EDMK; // for TD-DFT @@ -200,6 +205,7 @@ namespace elecstate * vector.size() = 2 for spin-polarization */ std::vector*> _DMR; + std::vector> _DMR_save; /** * @brief HContainer for density matrix in real space for gird parallelization diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index b350014051..24e7aa6367 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -609,6 +609,13 @@ namespace ModuleESolver { // save input rho this->pelec->charge->save_rho_before_sum_band(); + // save density matrix for mixing + if (GlobalV::MIXING_DMR && iter >= GlobalV::MIXING_RESTART) + { + elecstate::DensityMatrix* dm + = dynamic_cast*>(this->pelec)->get_DM(); + dm->save_DMR(); + } // using HSolverLCAO::solve() if (this->phsol != nullptr) From 621c2a72a308f90b308cc2a6a9b1966ce02b710c Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Thu, 25 Jan 2024 19:09:56 +0800 Subject: [PATCH 07/11] define a mix_dmr() in charge_mixing, and call it to mix density matrix in esolver_ks_lcao::eachiterfinish() --- .../module_charge/charge_mixing.cpp | 23 +++++++++++++++++++ .../module_charge/charge_mixing.h | 8 +++++++ .../module_dm/density_matrix.h | 2 ++ source/module_esolver/esolver_ks.cpp | 2 +- source/module_esolver/esolver_ks_lcao.cpp | 9 ++++++++ 5 files changed, 43 insertions(+), 1 deletion(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index f6ec11a4a5..7e8a1783ad 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -760,6 +760,29 @@ void Charge_Mixing::mix_rho_real(Charge* chr) } +void Charge_Mixing::mix_dmr(std::vector*> dmr, std::vector> dmr_save) +{ + // Notice that here I do not pass a DensityMatrix object to this function directly, + // since DensityMatrix is a Template class, which can not be used as a function parameter. + ModuleBase::TITLE("Charge_Mixing", "mix_dmr"); + ModuleBase::timer::tick("Charge_Mixing", "mix_dmr"); + // + const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1; + double* dmr_in; + double* dmr_out; + if (GlobalV::NSPIN == 1) + { + dmr_in = dmr_save[0].data(); + dmr_out = dmr[0]->get_wrapper(); + this->mixing->push_data(this->rho_mdata, dmr_in, dmr_out, nullptr, false); + this->mixing->mix_data(this->rho_mdata, dmr_out); + } + + ModuleBase::timer::tick("Charge_Mixing", "mix_dmr"); + + return; +} + void Charge_Mixing::mix_reset() { this->mixing->reset(); diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index b91b5698ce..b2d2f1fb5a 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -2,6 +2,7 @@ #ifndef CHARGE_MIXING_H #define CHARGE_MIXING_H #include "charge.h" +#include "module_elecstate/module_dm/density_matrix.h" #include "module_base/global_function.h" #include "module_base/global_variable.h" #include "module_base/module_mixing/mixing.h" @@ -32,6 +33,13 @@ class Charge_Mixing */ void mix_rho(Charge* chr); + /** + * @brief density matrix mixing, only for LCAO + * + */ + void mix_dmr(std::vector*> dmr, std::vector> dmr_save); + + /** * @brief charge mixing for reciprocal space * diff --git a/source/module_elecstate/module_dm/density_matrix.h b/source/module_elecstate/module_dm/density_matrix.h index b3aaeab186..92ba46e928 100644 --- a/source/module_elecstate/module_dm/density_matrix.h +++ b/source/module_elecstate/module_dm/density_matrix.h @@ -136,6 +136,8 @@ namespace elecstate */ std::vector*> get_DMR_vector() const; + std::vector> get_DMR_save() const {return _DMR_save;} + /** * @brief get pointer of DMK * @param ik k-point index, which is the index of _DMK diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index b4fcd83375..4af414600a 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -447,7 +447,7 @@ namespace ModuleESolver p_chgmix->mix_rho(pelec->charge); // update chr->rho by mixing } if (GlobalV::SCF_THR_TYPE == 2) pelec->charge->renormalize_rho(); // renormalize rho in R-space would induce a error in K-space - //----------charge mixing done----------- + //----------charge mixing done----------- } } #ifdef __MPI diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 24e7aa6367..476a887118 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -19,6 +19,7 @@ #include "module_hamilt_lcao/module_dftu/dftu.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/print_info.h" +#include "module_elecstate/module_dm/density_matrix.h" #ifdef __EXX // #include "module_rpa/rpa.h" #include "module_ri/RPA_LRI.h" @@ -796,6 +797,14 @@ namespace ModuleESolver template void ESolver_KS_LCAO::eachiterfinish(int iter) { + // mix density matrix + if (GlobalV::MIXING_RESTART > 0 && iter > GlobalV::MIXING_RESTART && GlobalV::MIXING_DMR ) + { + elecstate::DensityMatrix* dm + = dynamic_cast*>(this->pelec)->get_DM(); + this->p_chgmix->mix_dmr(dm->get_DMR_vector(), dm->get_DMR_save()); + } + //----------------------------------- // save charge density //----------------------------------- From baec07499e93213118a2ce47a99178444e1dab4c Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Fri, 26 Jan 2024 11:14:53 +0800 Subject: [PATCH 08/11] fix a bug --- .../module_charge/charge_mixing.cpp | 35 +++++++++++++++++-- .../module_charge/charge_mixing.h | 4 +-- source/module_esolver/esolver_ks_lcao.cpp | 6 ++-- 3 files changed, 37 insertions(+), 8 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 7e8a1783ad..e614ddcd00 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -132,6 +132,8 @@ void Charge_Mixing::allocate_mixing_dmr(int nnr) { this->mixing->init_mixing_data(this->dmr_mdata, nnr * dmr_nspin, sizeof(double)); } + + this->dmr_mdata.reset(); ModuleBase::timer::tick("Charge_Mixing", "allocate_mixing_dmr"); return; @@ -760,13 +762,15 @@ void Charge_Mixing::mix_rho_real(Charge* chr) } -void Charge_Mixing::mix_dmr(std::vector*> dmr, std::vector> dmr_save) +void Charge_Mixing::mix_dmr(elecstate::DensityMatrix* DM) { - // Notice that here I do not pass a DensityMatrix object to this function directly, - // since DensityMatrix is a Template class, which can not be used as a function parameter. + // Notice that DensityMatrix object is a Template class ModuleBase::TITLE("Charge_Mixing", "mix_dmr"); ModuleBase::timer::tick("Charge_Mixing", "mix_dmr"); // + std::vector*> dmr = DM->get_DMR_vector(); + std::vector> dmr_save = DM->get_DMR_save(); + // const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1; double* dmr_in; double* dmr_out; @@ -783,6 +787,31 @@ void Charge_Mixing::mix_dmr(std::vector*> dmr, std::v return; } +void Charge_Mixing::mix_dmr(elecstate::DensityMatrix, double>* DM) +{ + // Notice that DensityMatrix object is a Template class + ModuleBase::TITLE("Charge_Mixing", "mix_dmr"); + ModuleBase::timer::tick("Charge_Mixing", "mix_dmr"); + // + std::vector*> dmr = DM->get_DMR_vector(); + std::vector> dmr_save = DM->get_DMR_save(); + // + const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1; + double* dmr_in; + double* dmr_out; + if (GlobalV::NSPIN == 1) + { + dmr_in = dmr_save[0].data(); + dmr_out = dmr[0]->get_wrapper(); + this->mixing->push_data(this->dmr_mdata, dmr_in, dmr_out, nullptr, false); + this->mixing->mix_data(this->dmr_mdata, dmr_out); + } + + ModuleBase::timer::tick("Charge_Mixing", "mix_dmr"); + + return; +} + void Charge_Mixing::mix_reset() { this->mixing->reset(); diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index b2d2f1fb5a..8bfdeeb922 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -37,8 +37,8 @@ class Charge_Mixing * @brief density matrix mixing, only for LCAO * */ - void mix_dmr(std::vector*> dmr, std::vector> dmr_save); - + void mix_dmr(elecstate::DensityMatrix* DM); + void mix_dmr(elecstate::DensityMatrix, double>* DM); /** * @brief charge mixing for reciprocal space diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 476a887118..5919814a68 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -611,7 +611,7 @@ namespace ModuleESolver // save input rho this->pelec->charge->save_rho_before_sum_band(); // save density matrix for mixing - if (GlobalV::MIXING_DMR && iter >= GlobalV::MIXING_RESTART) + if (GlobalV::MIXING_RESTART > 0 && GlobalV::MIXING_DMR && iter >= GlobalV::MIXING_RESTART - 1) { elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); @@ -798,11 +798,11 @@ namespace ModuleESolver void ESolver_KS_LCAO::eachiterfinish(int iter) { // mix density matrix - if (GlobalV::MIXING_RESTART > 0 && iter > GlobalV::MIXING_RESTART && GlobalV::MIXING_DMR ) + if (GlobalV::MIXING_RESTART > 0 && iter >= GlobalV::MIXING_RESTART && GlobalV::MIXING_DMR ) { elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); - this->p_chgmix->mix_dmr(dm->get_DMR_vector(), dm->get_DMR_save()); + this->p_chgmix->mix_dmr(dm); } //----------------------------------- From a98c53016301d36df9723de096f18117897c5618 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Fri, 26 Jan 2024 16:40:43 +0800 Subject: [PATCH 09/11] implement mix_dmr for npsin=2 --- .../module_charge/charge_mixing.cpp | 160 +++++++++++++++++- source/module_esolver/esolver_ks_lcao.cpp | 3 +- 2 files changed, 155 insertions(+), 8 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index e614ddcd00..4919e9b5c9 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -771,15 +771,89 @@ void Charge_Mixing::mix_dmr(elecstate::DensityMatrix* DM) std::vector*> dmr = DM->get_DMR_vector(); std::vector> dmr_save = DM->get_DMR_save(); // - const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1; + //const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1; double* dmr_in; double* dmr_out; - if (GlobalV::NSPIN == 1) + if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4) { dmr_in = dmr_save[0].data(); dmr_out = dmr[0]->get_wrapper(); - this->mixing->push_data(this->rho_mdata, dmr_in, dmr_out, nullptr, false); - this->mixing->mix_data(this->rho_mdata, dmr_out); + this->mixing->push_data(this->dmr_mdata, dmr_in, dmr_out, nullptr, false); + this->mixing->mix_data(this->dmr_mdata, dmr_out); + } + else if (GlobalV::NSPIN == 2) + { + // magnetic density matrix + double* dmr_mag = nullptr; + double* dmr_mag_save = nullptr; + const int nnr = dmr[0]->get_nnr(); + // allocate dmr_mag[is*nnrx] and dmr_mag_save[is*nnrx] + dmr_mag = new double[nnr * GlobalV::NSPIN]; + dmr_mag_save = new double[nnr * GlobalV::NSPIN]; + ModuleBase::GlobalFunc::ZEROS(dmr_mag, nnr * GlobalV::NSPIN); + ModuleBase::GlobalFunc::ZEROS(dmr_mag_save, nnr * GlobalV::NSPIN); + double* dmr_up; + double* dmr_down; + // tranfer dmr into dmr_mag + dmr_up = dmr[0]->get_wrapper(); + dmr_down = dmr[1]->get_wrapper(); + for (int ir = 0; ir < nnr; ir++) + { + dmr_mag[ir] = dmr_up[ir] + dmr_down[ir]; + dmr_mag[ir + nnr] = dmr_up[ir] - dmr_down[ir]; + } + // tranfer dmr_save into dmr_mag_save + dmr_up = dmr_save[0].data(); + dmr_down = dmr_save[1].data(); + for (int ir = 0; ir < nnr; ir++) + { + dmr_mag_save[ir] = dmr_up[ir] + dmr_down[ir]; + dmr_mag_save[ir + nnr] = dmr_up[ir] - dmr_down[ir]; + } + // + dmr_in = dmr_mag_save; + dmr_out = dmr_mag; + // no kerker in mixing_dmr + //auto screen = std::bind(&Charge_Mixing::Kerker_screen_real, this, std::placeholders::_1); + auto twobeta_mix + = [this, nnr](double* out, const double* in, const double* sres) { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = 0; i < nnr; ++i) + { + out[i] = in[i] + this->mixing_beta * sres[i]; + } + // magnetism +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = nnr; i < 2 * nnr; ++i) + { + out[i] = in[i] + this->mixing_beta_mag * sres[i]; + } + }; + this->mixing->push_data(this->dmr_mdata, dmr_in, dmr_out, nullptr, twobeta_mix, false); + //auto inner_product + // = std::bind(&Charge_Mixing::inner_product_real, this, std::placeholders::_1, std::placeholders::_2); + //this->mixing->cal_coef(this->rho_mdata, inner_product); + this->mixing->mix_data(this->dmr_mdata, dmr_out); + // get new dmr from dmr_mag + dmr_up = dmr[0]->get_wrapper(); + dmr_down = dmr[1]->get_wrapper(); + for (int is = 0; is < GlobalV::NSPIN; is++) + { + ModuleBase::GlobalFunc::ZEROS(dmr_up, nnr); + ModuleBase::GlobalFunc::ZEROS(dmr_down, nnr); + } + for (int ir = 0; ir < nnr; ir++) + { + dmr_up[ir] = 0.5 * (dmr_mag[ir] + dmr_mag[ir+nnr]); + dmr_down[ir] = 0.5 * (dmr_mag[ir] - dmr_mag[ir+nnr]); + } + // delete + delete[] dmr_mag; + delete[] dmr_mag_save; } ModuleBase::timer::tick("Charge_Mixing", "mix_dmr"); @@ -796,16 +870,90 @@ void Charge_Mixing::mix_dmr(elecstate::DensityMatrix, doubl std::vector*> dmr = DM->get_DMR_vector(); std::vector> dmr_save = DM->get_DMR_save(); // - const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1; + //const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1; double* dmr_in; double* dmr_out; - if (GlobalV::NSPIN == 1) + if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4) { dmr_in = dmr_save[0].data(); dmr_out = dmr[0]->get_wrapper(); this->mixing->push_data(this->dmr_mdata, dmr_in, dmr_out, nullptr, false); this->mixing->mix_data(this->dmr_mdata, dmr_out); } + else if (GlobalV::NSPIN == 2) + { + // magnetic density matrix + double* dmr_mag = nullptr; + double* dmr_mag_save = nullptr; + const int nnr = dmr[0]->get_nnr(); + // allocate dmr_mag[is*nnrx] and dmr_mag_save[is*nnrx] + dmr_mag = new double[nnr * GlobalV::NSPIN]; + dmr_mag_save = new double[nnr * GlobalV::NSPIN]; + ModuleBase::GlobalFunc::ZEROS(dmr_mag, nnr * GlobalV::NSPIN); + ModuleBase::GlobalFunc::ZEROS(dmr_mag_save, nnr * GlobalV::NSPIN); + double* dmr_up; + double* dmr_down; + // tranfer dmr into dmr_mag + dmr_up = dmr[0]->get_wrapper(); + dmr_down = dmr[1]->get_wrapper(); + for (int ir = 0; ir < nnr; ir++) + { + dmr_mag[ir] = dmr_up[ir] + dmr_down[ir]; + dmr_mag[ir + nnr] = dmr_up[ir] - dmr_down[ir]; + } + // tranfer dmr_save into dmr_mag_save + dmr_up = dmr_save[0].data(); + dmr_down = dmr_save[1].data(); + for (int ir = 0; ir < nnr; ir++) + { + dmr_mag_save[ir] = dmr_up[ir] + dmr_down[ir]; + dmr_mag_save[ir + nnr] = dmr_up[ir] - dmr_down[ir]; + } + // + dmr_in = dmr_mag_save; + dmr_out = dmr_mag; + // no kerker in mixing_dmr + //auto screen = std::bind(&Charge_Mixing::Kerker_screen_real, this, std::placeholders::_1); + auto twobeta_mix + = [this, nnr](double* out, const double* in, const double* sres) { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = 0; i < nnr; ++i) + { + out[i] = in[i] + this->mixing_beta * sres[i]; + } + // magnetism +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (int i = nnr; i < 2 * nnr; ++i) + { + out[i] = in[i] + this->mixing_beta_mag * sres[i]; + } + }; + this->mixing->push_data(this->dmr_mdata, dmr_in, dmr_out, nullptr, twobeta_mix, false); + //auto inner_product + // = std::bind(&Charge_Mixing::inner_product_real, this, std::placeholders::_1, std::placeholders::_2); + //this->mixing->cal_coef(this->rho_mdata, inner_product); + this->mixing->mix_data(this->dmr_mdata, dmr_out); + // get new dmr from dmr_mag + dmr_up = dmr[0]->get_wrapper(); + dmr_down = dmr[1]->get_wrapper(); + for (int is = 0; is < GlobalV::NSPIN; is++) + { + ModuleBase::GlobalFunc::ZEROS(dmr_up, nnr); + ModuleBase::GlobalFunc::ZEROS(dmr_down, nnr); + } + for (int ir = 0; ir < nnr; ir++) + { + dmr_up[ir] = 0.5 * (dmr_mag[ir] + dmr_mag[ir+nnr]); + dmr_down[ir] = 0.5 * (dmr_mag[ir] - dmr_mag[ir+nnr]); + } + // delete + delete[] dmr_mag; + delete[] dmr_mag_save; + } ModuleBase::timer::tick("Charge_Mixing", "mix_dmr"); diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 5919814a68..89f57864af 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -19,7 +19,6 @@ #include "module_hamilt_lcao/module_dftu/dftu.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/print_info.h" -#include "module_elecstate/module_dm/density_matrix.h" #ifdef __EXX // #include "module_rpa/rpa.h" #include "module_ri/RPA_LRI.h" @@ -611,7 +610,7 @@ namespace ModuleESolver // save input rho this->pelec->charge->save_rho_before_sum_band(); // save density matrix for mixing - if (GlobalV::MIXING_RESTART > 0 && GlobalV::MIXING_DMR && iter >= GlobalV::MIXING_RESTART - 1) + if (GlobalV::MIXING_RESTART > 0 && GlobalV::MIXING_DMR && iter >= GlobalV::MIXING_RESTART) { elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); From cdd8b953a55e91365640aac4d2e5272f3ca8d9e1 Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Fri, 26 Jan 2024 17:55:54 +0800 Subject: [PATCH 10/11] move some get_function from .cpp to .h, namely density_matrix.h, hcontainer.h, and atom_pair.h --- .../module_dm/density_matrix.cpp | 7 ------- .../module_dm/density_matrix.h | 2 +- .../module_hcontainer/atom_pair.cpp | 18 ------------------ .../module_hcontainer/atom_pair.h | 11 +++++++++-- .../module_hcontainer/hcontainer.cpp | 19 ------------------- .../module_hcontainer/hcontainer.h | 12 ++++++++++-- 6 files changed, 20 insertions(+), 49 deletions(-) diff --git a/source/module_elecstate/module_dm/density_matrix.cpp b/source/module_elecstate/module_dm/density_matrix.cpp index b7e8198bed..fd6b3a2af8 100644 --- a/source/module_elecstate/module_dm/density_matrix.cpp +++ b/source/module_elecstate/module_dm/density_matrix.cpp @@ -263,13 +263,6 @@ hamilt::HContainer* DensityMatrix::get_DMR_pointer(const int ispin) return this->_DMR[ispin - 1]; } -// get _DMR pointer vector -template -std::vector*> DensityMatrix::get_DMR_vector() const -{ - return this->_DMR; -} - // get _DMK[ik] pointer template TK* DensityMatrix::get_DMK_pointer(const int ik) const diff --git a/source/module_elecstate/module_dm/density_matrix.h b/source/module_elecstate/module_dm/density_matrix.h index 92ba46e928..519f798f5d 100644 --- a/source/module_elecstate/module_dm/density_matrix.h +++ b/source/module_elecstate/module_dm/density_matrix.h @@ -134,7 +134,7 @@ namespace elecstate * @brief get pointer vector of DMR * @return HContainer* vector of DMR */ - std::vector*> get_DMR_vector() const; + std::vector*> get_DMR_vector() const {return this->_DMR;} std::vector> get_DMR_save() const {return _DMR_save;} diff --git a/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp b/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp index 060f4a7ac2..d308bfd8e9 100644 --- a/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp @@ -317,13 +317,6 @@ void AtomPair::set_size(const int& col_size_in, const int& row_size_in) this->row_size = row_size_in; } -// get size -template -int AtomPair::get_size() const -{ - return this->col_size * this->row_size; -} - // get paraV for check template const Parallel_Orbitals* AtomPair::get_paraV() const @@ -788,17 +781,6 @@ T* AtomPair::get_pointer(int ir) const return this->values[ir].get_pointer(); } -// get_R_size -template -size_t AtomPair::get_R_size() const -{ -#ifdef __DEBUG - assert(this->R_index.size() / 3 == this->values.size()); - assert(this->R_index.size() % 3 == 0); -#endif - return this->R_index.size() / 3; -} - // get_memory_size template size_t AtomPair::get_memory_size() const diff --git a/source/module_hamilt_lcao/module_hcontainer/atom_pair.h b/source/module_hamilt_lcao/module_hcontainer/atom_pair.h index c0dbc905d6..309231d79b 100644 --- a/source/module_hamilt_lcao/module_hcontainer/atom_pair.h +++ b/source/module_hamilt_lcao/module_hcontainer/atom_pair.h @@ -124,7 +124,7 @@ class AtomPair * @brief get size = col_size * row_size * @return int */ - int get_size() const; + int get_size() const {return this->col_size * this->row_size;} /** * @brief get Parallel_Orbitals pointer of this AtomPair for checking 2d-block parallel @@ -256,7 +256,14 @@ class AtomPair AtomPair& operator=(AtomPair&& other) noexcept; // interface for getting the size of this->R_index - size_t get_R_size() const; + size_t get_R_size() const + { +#ifdef __DEBUG + assert(this->R_index.size() / 3 == this->values.size()); + assert(this->R_index.size() % 3 == 0); +#endif + return this->R_index.size() / 3; + } /** * @brief get total memory size of AtomPair diff --git a/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp b/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp index 6867c84b70..4cde27154d 100644 --- a/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp @@ -601,25 +601,6 @@ size_t HContainer::get_memory_size() const return memory; } -// get_nnr -template -size_t HContainer::get_nnr() const -{ - size_t sum = 0; - for(int iap=0;iap < this->atom_pairs.size();++iap) - { - sum += this->atom_pairs[iap].get_R_size() * this->atom_pairs[iap].get_size(); - } - return sum; -} - -// get_wrapper -template -T* HContainer::get_wrapper() const -{ - return this->wrapper_pointer; -} - // synchronize template void HContainer::shape_synchron( const HContainer& other) diff --git a/source/module_hamilt_lcao/module_hcontainer/hcontainer.h b/source/module_hamilt_lcao/module_hcontainer/hcontainer.h index e982b6e1d3..43eddf0f53 100644 --- a/source/module_hamilt_lcao/module_hcontainer/hcontainer.h +++ b/source/module_hamilt_lcao/module_hcontainer/hcontainer.h @@ -354,7 +354,15 @@ class HContainer * named nnr inherited from history * all AtomPairs and BaseMatrixs are counted */ - size_t get_nnr() const; + size_t get_nnr() const + { + size_t sum = 0; + for(int iap=0;iap < this->atom_pairs.size();++iap) + { + sum += this->atom_pairs[iap].get_R_size() * this->atom_pairs[iap].get_size(); + } + return sum; + } /** * @brief get infomation of IJR pairs in HContainer @@ -373,7 +381,7 @@ class HContainer /** * @brief return the wrapper_pointer */ - T* get_wrapper() const; + T* get_wrapper() const {return this->wrapper_pointer;} /** * @brief synchronization of atom-pairs for read-in HContainer From 8c6f167d6255a27d74a97dde68baccf0e27c746d Mon Sep 17 00:00:00 2001 From: weiqingzhou Date: Sun, 28 Jan 2024 16:14:40 +0800 Subject: [PATCH 11/11] add the docs about mixing_dmr --- docs/advanced/input_files/input-main.md | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 8c675cab92..56b8a5aa81 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -73,6 +73,7 @@ - [mixing\_beta\_mag](#mixing_beta_mag) - [mixing\_ndim](#mixing_ndim) - [mixing\_restart](#mixing_restart) + - [mixing\_dmr](#mixing_dmr) - [mixing\_gg0](#mixing_gg0) - [mixing\_gg0\_mag](#mixing_gg0_mag) - [mixing\_gg0\_min](#mixing_gg0_min) @@ -1013,6 +1014,14 @@ We recommend the following options: - **Default**: 0 +### mixing_dmr + +- **Type**: bool +- **Availability**: Only for `mixing_restart>=2` +- **Description**: At `mixing_restart`-th iteration, SCF will start a mixing for real-space density matrix by using the same coefficiences as the mixing of charge density. + +- **Default**: false + ### mixing_gg0 - **Type**: Real