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 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_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 227d9a050a..4919e9b5c9 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -114,6 +114,31 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, return; } +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"); + // + 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)); + } + + this->dmr_mdata.reset(); + 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; @@ -737,6 +762,204 @@ void Charge_Mixing::mix_rho_real(Charge* chr) } +void Charge_Mixing::mix_dmr(elecstate::DensityMatrix* 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 || 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"); + + 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 || 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"); + + 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 993c229057..8bfdeeb922 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" @@ -16,6 +17,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. @@ -31,6 +33,13 @@ class Charge_Mixing */ void mix_rho(Charge* chr); + /** + * @brief density matrix mixing, only for LCAO + * + */ + void mix_dmr(elecstate::DensityMatrix* DM); + void mix_dmr(elecstate::DensityMatrix, double>* DM); + /** * @brief charge mixing for reciprocal space * @@ -56,7 +65,6 @@ class Charge_Mixing * */ void Kerker_screen_real(double* rho); - void Kerker_screen_real_test(double* rho); /** * @brief Inner product of two complex vectors @@ -87,19 +95,13 @@ 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 allocate memory of dmr_mdata + * + */ + void allocate_mixing_dmr(int nnr); /** * @brief Get the drho diff --git a/source/module_elecstate/module_dm/density_matrix.cpp b/source/module_elecstate/module_dm/density_matrix.cpp index c284099135..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 @@ -382,6 +375,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..519f798f5d 100644 --- a/source/module_elecstate/module_dm/density_matrix.h +++ b/source/module_elecstate/module_dm/density_matrix.h @@ -134,7 +134,9 @@ 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;} /** * @brief get pointer of DMK @@ -190,6 +192,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 +207,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.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 4c4d6c342f..89f57864af 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(); } @@ -601,6 +609,13 @@ 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) + { + elecstate::DensityMatrix* dm + = dynamic_cast*>(this->pelec)->get_DM(); + dm->save_DMR(); + } // using HSolverLCAO::solve() if (this->phsol != nullptr) @@ -781,6 +796,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); + } + //----------------------------------- // save charge density //----------------------------------- 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 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");