Skip to content

Commit 2864d6e

Browse files
WHUweiqingzhouWHUweiqingzhou
authored
Feature: mix real-space density matrix for better convergence of DFT+U (#3549)
* add a new input parameter mixing_dmr and corresponding UnitTests * delete Kerker_screen_real_test() * delete autoset() in charge_mixing.h * declare a new pubilc member dmr_mdata, whose memory can be allocated by allocate_mixing_dmr() * allocate dmr_mdata in ESolver_KS_LCAO::eachiterinit() * define a function save_DMR in DensityMatrix, and save perious DMR in ESolver_KS_LCAO::hamilt2density() * define a mix_dmr() in charge_mixing, and call it to mix density matrix in esolver_ks_lcao::eachiterfinish() * fix a bug * implement mix_dmr for npsin=2 * move some get_function from .cpp to .h, namely density_matrix.h, hcontainer.h, and atom_pair.h * add the docs about mixing_dmr
1 parent 85e0234 commit 2864d6e

File tree

21 files changed

+351
-64
lines changed

21 files changed

+351
-64
lines changed

‎docs/advanced/input_files/input-main.md‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,7 @@
7373
- [mixing\_beta\_mag](#mixing_beta_mag)
7474
- [mixing\_ndim](#mixing_ndim)
7575
- [mixing\_restart](#mixing_restart)
76+
- [mixing\_dmr](#mixing_dmr)
7677
- [mixing\_gg0](#mixing_gg0)
7778
- [mixing\_gg0\_mag](#mixing_gg0_mag)
7879
- [mixing\_gg0\_min](#mixing_gg0_min)
@@ -1013,6 +1014,14 @@ We recommend the following options:
10131014

10141015
- **Default**: 0
10151016

1017+
### mixing_dmr
1018+
1019+
- **Type**: bool
1020+
- **Availability**: Only for `mixing_restart>=2`
1021+
- **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.
1022+
1023+
- **Default**: false
1024+
10161025
### mixing_gg0
10171026

10181027
- **Type**: Real

‎source/module_base/global_variable.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -255,6 +255,7 @@ double MIXING_GG0_MAG = 1.00;
255255
double MIXING_GG0_MIN = 0.1;
256256
double MIXING_ANGLE = 0.0;
257257
bool MIXING_TAU = 0;
258+
bool MIXING_DMR = 0;
258259

259260
//==========================================================
260261
// device flags added by denghui

‎source/module_base/global_variable.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -284,6 +284,7 @@ extern double MIXING_BETA_MAG;
284284
extern double MIXING_GG0_MAG;
285285
extern double MIXING_GG0_MIN;
286286
extern double MIXING_ANGLE;
287+
extern bool MIXING_DMR;
287288

288289
//==========================================================
289290
// device flags added by denghui

‎source/module_elecstate/module_charge/charge_mixing.cpp‎

Lines changed: 223 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -114,6 +114,31 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in,
114114
return;
115115
}
116116

117+
void Charge_Mixing::allocate_mixing_dmr(int nnr)
118+
{
119+
// Note that: we cannot allocate memory for dmr_mdata in set_mixing.
120+
// since the size of dmr_mdata is given by the size of HContainer.nnr, which is calculated in DensityMatrix::init_DMR().
121+
// and DensityMatrix::init_DMR() is called in beforescf(). While set_mixing() is called in ESolver_KS::Init().
122+
ModuleBase::TITLE("Charge_Mixing", "allocate_mixing_dmr");
123+
ModuleBase::timer::tick("Charge_Mixing", "allocate_mixing_dmr");
124+
//
125+
const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1;
126+
// allocate memory for dmr_mdata
127+
if (GlobalV::SCF_THR_TYPE == 1)
128+
{
129+
ModuleBase::WARNING_QUIT("Charge_Mixing", "This Mixing of Density Matrix is not supported for PW basis yet");
130+
}
131+
else if (GlobalV::SCF_THR_TYPE == 2)
132+
{
133+
this->mixing->init_mixing_data(this->dmr_mdata, nnr * dmr_nspin, sizeof(double));
134+
}
135+
136+
this->dmr_mdata.reset();
137+
ModuleBase::timer::tick("Charge_Mixing", "allocate_mixing_dmr");
138+
139+
return;
140+
}
141+
117142
void Charge_Mixing::set_rhopw(ModulePW::PW_Basis* rhopw_in, ModulePW::PW_Basis* rhodpw_in)
118143
{
119144
this->rhopw = rhopw_in;
@@ -737,6 +762,204 @@ void Charge_Mixing::mix_rho_real(Charge* chr)
737762

738763
}
739764

765+
void Charge_Mixing::mix_dmr(elecstate::DensityMatrix<double, double>* DM)
766+
{
767+
// Notice that DensityMatrix object is a Template class
768+
ModuleBase::TITLE("Charge_Mixing", "mix_dmr");
769+
ModuleBase::timer::tick("Charge_Mixing", "mix_dmr");
770+
//
771+
std::vector<hamilt::HContainer<double>*> dmr = DM->get_DMR_vector();
772+
std::vector<std::vector<double>> dmr_save = DM->get_DMR_save();
773+
//
774+
//const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1;
775+
double* dmr_in;
776+
double* dmr_out;
777+
if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
778+
{
779+
dmr_in = dmr_save[0].data();
780+
dmr_out = dmr[0]->get_wrapper();
781+
this->mixing->push_data(this->dmr_mdata, dmr_in, dmr_out, nullptr, false);
782+
this->mixing->mix_data(this->dmr_mdata, dmr_out);
783+
}
784+
else if (GlobalV::NSPIN == 2)
785+
{
786+
// magnetic density matrix
787+
double* dmr_mag = nullptr;
788+
double* dmr_mag_save = nullptr;
789+
const int nnr = dmr[0]->get_nnr();
790+
// allocate dmr_mag[is*nnrx] and dmr_mag_save[is*nnrx]
791+
dmr_mag = new double[nnr * GlobalV::NSPIN];
792+
dmr_mag_save = new double[nnr * GlobalV::NSPIN];
793+
ModuleBase::GlobalFunc::ZEROS(dmr_mag, nnr * GlobalV::NSPIN);
794+
ModuleBase::GlobalFunc::ZEROS(dmr_mag_save, nnr * GlobalV::NSPIN);
795+
double* dmr_up;
796+
double* dmr_down;
797+
// tranfer dmr into dmr_mag
798+
dmr_up = dmr[0]->get_wrapper();
799+
dmr_down = dmr[1]->get_wrapper();
800+
for (int ir = 0; ir < nnr; ir++)
801+
{
802+
dmr_mag[ir] = dmr_up[ir] + dmr_down[ir];
803+
dmr_mag[ir + nnr] = dmr_up[ir] - dmr_down[ir];
804+
}
805+
// tranfer dmr_save into dmr_mag_save
806+
dmr_up = dmr_save[0].data();
807+
dmr_down = dmr_save[1].data();
808+
for (int ir = 0; ir < nnr; ir++)
809+
{
810+
dmr_mag_save[ir] = dmr_up[ir] + dmr_down[ir];
811+
dmr_mag_save[ir + nnr] = dmr_up[ir] - dmr_down[ir];
812+
}
813+
//
814+
dmr_in = dmr_mag_save;
815+
dmr_out = dmr_mag;
816+
// no kerker in mixing_dmr
817+
//auto screen = std::bind(&Charge_Mixing::Kerker_screen_real, this, std::placeholders::_1);
818+
auto twobeta_mix
819+
= [this, nnr](double* out, const double* in, const double* sres) {
820+
#ifdef _OPENMP
821+
#pragma omp parallel for schedule(static, 256)
822+
#endif
823+
for (int i = 0; i < nnr; ++i)
824+
{
825+
out[i] = in[i] + this->mixing_beta * sres[i];
826+
}
827+
// magnetism
828+
#ifdef _OPENMP
829+
#pragma omp parallel for schedule(static, 256)
830+
#endif
831+
for (int i = nnr; i < 2 * nnr; ++i)
832+
{
833+
out[i] = in[i] + this->mixing_beta_mag * sres[i];
834+
}
835+
};
836+
this->mixing->push_data(this->dmr_mdata, dmr_in, dmr_out, nullptr, twobeta_mix, false);
837+
//auto inner_product
838+
// = std::bind(&Charge_Mixing::inner_product_real, this, std::placeholders::_1, std::placeholders::_2);
839+
//this->mixing->cal_coef(this->rho_mdata, inner_product);
840+
this->mixing->mix_data(this->dmr_mdata, dmr_out);
841+
// get new dmr from dmr_mag
842+
dmr_up = dmr[0]->get_wrapper();
843+
dmr_down = dmr[1]->get_wrapper();
844+
for (int is = 0; is < GlobalV::NSPIN; is++)
845+
{
846+
ModuleBase::GlobalFunc::ZEROS(dmr_up, nnr);
847+
ModuleBase::GlobalFunc::ZEROS(dmr_down, nnr);
848+
}
849+
for (int ir = 0; ir < nnr; ir++)
850+
{
851+
dmr_up[ir] = 0.5 * (dmr_mag[ir] + dmr_mag[ir+nnr]);
852+
dmr_down[ir] = 0.5 * (dmr_mag[ir] - dmr_mag[ir+nnr]);
853+
}
854+
// delete
855+
delete[] dmr_mag;
856+
delete[] dmr_mag_save;
857+
}
858+
859+
ModuleBase::timer::tick("Charge_Mixing", "mix_dmr");
860+
861+
return;
862+
}
863+
864+
void Charge_Mixing::mix_dmr(elecstate::DensityMatrix<std::complex<double>, double>* DM)
865+
{
866+
// Notice that DensityMatrix object is a Template class
867+
ModuleBase::TITLE("Charge_Mixing", "mix_dmr");
868+
ModuleBase::timer::tick("Charge_Mixing", "mix_dmr");
869+
//
870+
std::vector<hamilt::HContainer<double>*> dmr = DM->get_DMR_vector();
871+
std::vector<std::vector<double>> dmr_save = DM->get_DMR_save();
872+
//
873+
//const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1;
874+
double* dmr_in;
875+
double* dmr_out;
876+
if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
877+
{
878+
dmr_in = dmr_save[0].data();
879+
dmr_out = dmr[0]->get_wrapper();
880+
this->mixing->push_data(this->dmr_mdata, dmr_in, dmr_out, nullptr, false);
881+
this->mixing->mix_data(this->dmr_mdata, dmr_out);
882+
}
883+
else if (GlobalV::NSPIN == 2)
884+
{
885+
// magnetic density matrix
886+
double* dmr_mag = nullptr;
887+
double* dmr_mag_save = nullptr;
888+
const int nnr = dmr[0]->get_nnr();
889+
// allocate dmr_mag[is*nnrx] and dmr_mag_save[is*nnrx]
890+
dmr_mag = new double[nnr * GlobalV::NSPIN];
891+
dmr_mag_save = new double[nnr * GlobalV::NSPIN];
892+
ModuleBase::GlobalFunc::ZEROS(dmr_mag, nnr * GlobalV::NSPIN);
893+
ModuleBase::GlobalFunc::ZEROS(dmr_mag_save, nnr * GlobalV::NSPIN);
894+
double* dmr_up;
895+
double* dmr_down;
896+
// tranfer dmr into dmr_mag
897+
dmr_up = dmr[0]->get_wrapper();
898+
dmr_down = dmr[1]->get_wrapper();
899+
for (int ir = 0; ir < nnr; ir++)
900+
{
901+
dmr_mag[ir] = dmr_up[ir] + dmr_down[ir];
902+
dmr_mag[ir + nnr] = dmr_up[ir] - dmr_down[ir];
903+
}
904+
// tranfer dmr_save into dmr_mag_save
905+
dmr_up = dmr_save[0].data();
906+
dmr_down = dmr_save[1].data();
907+
for (int ir = 0; ir < nnr; ir++)
908+
{
909+
dmr_mag_save[ir] = dmr_up[ir] + dmr_down[ir];
910+
dmr_mag_save[ir + nnr] = dmr_up[ir] - dmr_down[ir];
911+
}
912+
//
913+
dmr_in = dmr_mag_save;
914+
dmr_out = dmr_mag;
915+
// no kerker in mixing_dmr
916+
//auto screen = std::bind(&Charge_Mixing::Kerker_screen_real, this, std::placeholders::_1);
917+
auto twobeta_mix
918+
= [this, nnr](double* out, const double* in, const double* sres) {
919+
#ifdef _OPENMP
920+
#pragma omp parallel for schedule(static, 256)
921+
#endif
922+
for (int i = 0; i < nnr; ++i)
923+
{
924+
out[i] = in[i] + this->mixing_beta * sres[i];
925+
}
926+
// magnetism
927+
#ifdef _OPENMP
928+
#pragma omp parallel for schedule(static, 256)
929+
#endif
930+
for (int i = nnr; i < 2 * nnr; ++i)
931+
{
932+
out[i] = in[i] + this->mixing_beta_mag * sres[i];
933+
}
934+
};
935+
this->mixing->push_data(this->dmr_mdata, dmr_in, dmr_out, nullptr, twobeta_mix, false);
936+
//auto inner_product
937+
// = std::bind(&Charge_Mixing::inner_product_real, this, std::placeholders::_1, std::placeholders::_2);
938+
//this->mixing->cal_coef(this->rho_mdata, inner_product);
939+
this->mixing->mix_data(this->dmr_mdata, dmr_out);
940+
// get new dmr from dmr_mag
941+
dmr_up = dmr[0]->get_wrapper();
942+
dmr_down = dmr[1]->get_wrapper();
943+
for (int is = 0; is < GlobalV::NSPIN; is++)
944+
{
945+
ModuleBase::GlobalFunc::ZEROS(dmr_up, nnr);
946+
ModuleBase::GlobalFunc::ZEROS(dmr_down, nnr);
947+
}
948+
for (int ir = 0; ir < nnr; ir++)
949+
{
950+
dmr_up[ir] = 0.5 * (dmr_mag[ir] + dmr_mag[ir+nnr]);
951+
dmr_down[ir] = 0.5 * (dmr_mag[ir] - dmr_mag[ir+nnr]);
952+
}
953+
// delete
954+
delete[] dmr_mag;
955+
delete[] dmr_mag_save;
956+
}
957+
958+
ModuleBase::timer::tick("Charge_Mixing", "mix_dmr");
959+
960+
return;
961+
}
962+
740963
void Charge_Mixing::mix_reset()
741964
{
742965
this->mixing->reset();

‎source/module_elecstate/module_charge/charge_mixing.h‎

Lines changed: 16 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
#ifndef CHARGE_MIXING_H
33
#define CHARGE_MIXING_H
44
#include "charge.h"
5+
#include "module_elecstate/module_dm/density_matrix.h"
56
#include "module_base/global_function.h"
67
#include "module_base/global_variable.h"
78
#include "module_base/module_mixing/mixing.h"
@@ -16,6 +17,7 @@ class Charge_Mixing
1617
Base_Mixing::Mixing_Data rho_mdata; ///< Mixing data for charge density
1718
Base_Mixing::Mixing_Data tau_mdata; ///< Mixing data for kinetic energy density
1819
Base_Mixing::Mixing_Data nhat_mdata; ///< Mixing data for compensation density
20+
Base_Mixing::Mixing_Data dmr_mdata; ///< Mixing data for real space density matrix
1921

2022
Base_Mixing::Plain_Mixing* mixing_highf = nullptr; ///< The high_frequency part is mixed by plain mixing method.
2123

@@ -31,6 +33,13 @@ class Charge_Mixing
3133
*/
3234
void mix_rho(Charge* chr);
3335

36+
/**
37+
* @brief density matrix mixing, only for LCAO
38+
*
39+
*/
40+
void mix_dmr(elecstate::DensityMatrix<double, double>* DM);
41+
void mix_dmr(elecstate::DensityMatrix<std::complex<double>, double>* DM);
42+
3443
/**
3544
* @brief charge mixing for reciprocal space
3645
*
@@ -56,7 +65,6 @@ class Charge_Mixing
5665
*
5766
*/
5867
void Kerker_screen_real(double* rho);
59-
void Kerker_screen_real_test(double* rho);
6068

6169
/**
6270
* @brief Inner product of two complex vectors
@@ -87,19 +95,13 @@ class Charge_Mixing
8795
const int& mixing_ndim_in,
8896
const double& mixing_gg0_in,
8997
const bool& mixing_tau_in,
90-
const double& mixing_beta_mag_in); // mohan add mixing_gg0_in 2014-09-27
91-
92-
// /**
93-
// * @brief use auto set
94-
// *
95-
// */
96-
// void need_auto_set();
97-
98-
// /**
99-
// * @brief auto set mixing gg0 and mixing_beta
100-
// *
101-
// */
102-
// void auto_set(const double& bandgap_in, const UnitCell& ucell_);
98+
const double& mixing_beta_mag_in);
99+
100+
/**
101+
* @brief allocate memory of dmr_mdata
102+
*
103+
*/
104+
void allocate_mixing_dmr(int nnr);
103105

104106
/**
105107
* @brief Get the drho

‎source/module_elecstate/module_dm/density_matrix.cpp‎

Lines changed: 32 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -263,13 +263,6 @@ hamilt::HContainer<TR>* DensityMatrix<TK, TR>::get_DMR_pointer(const int ispin)
263263
return this->_DMR[ispin - 1];
264264
}
265265

266-
// get _DMR pointer vector
267-
template <typename TK, typename TR>
268-
std::vector<hamilt::HContainer<TR>*> DensityMatrix<TK, TR>::get_DMR_vector() const
269-
{
270-
return this->_DMR;
271-
}
272-
273266
// get _DMK[ik] pointer
274267
template <typename TK, typename TR>
275268
TK* DensityMatrix<TK, TR>::get_DMK_pointer(const int ik) const
@@ -382,6 +375,38 @@ int DensityMatrix<TK, TR>::get_DMK_ncol() const
382375
return this->_paraV->ncol;
383376
}
384377

378+
template <typename TK, typename TR>
379+
void DensityMatrix<TK, TR>::save_DMR()
380+
{
381+
ModuleBase::TITLE("DensityMatrix", "save_DMR");
382+
ModuleBase::timer::tick("DensityMatrix", "save_DMR");
383+
384+
const int nnr = this->_DMR[0]->get_nnr();
385+
// allocate if _DMR_save is empty
386+
if(_DMR_save.size() == 0)
387+
{
388+
_DMR_save.resize(this->_DMR.size());
389+
for(int is=0;is<this->_DMR.size();is++)
390+
{
391+
_DMR_save[is].resize(nnr);
392+
}
393+
}
394+
// save _DMR to _DMR_save
395+
for(int is=0;is<this->_DMR.size();is++)
396+
{
397+
TR* DMR_pointer = this->_DMR[is]->get_wrapper();
398+
TR* DMR_save_pointer = _DMR_save[is].data();
399+
// set to zero
400+
ModuleBase::GlobalFunc::ZEROS(DMR_save_pointer, nnr);
401+
for(int i=0;i<nnr;i++)
402+
{
403+
DMR_save_pointer[i] = DMR_pointer[i];
404+
}
405+
}
406+
407+
ModuleBase::timer::tick("DensityMatrix", "save_DMR");
408+
}
409+
385410
// calculate DMR from DMK using add_element
386411
template <typename TK, typename TR>
387412
void DensityMatrix<TK,TR>::cal_DMR_test()

0 commit comments

Comments
 (0)