Feature: mix real-space density matrix for better convergence of DFT+U

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

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

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

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<double, 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<hamilt::HContainer<double>*> dmr = DM->get_DMR_vector();
+    std::vector<std::vector<double>> 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<std::complex<double>, 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<hamilt::HContainer<double>*> dmr = DM->get_DMR_vector();
+    std::vector<std::vector<double>> 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();

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<double, double>* DM);
+    void mix_dmr(elecstate::DensityMatrix<std::complex<double>, 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

source/module_elecstate/module_dm/density_matrix.cpp

@@ -263,13 +263,6 @@ hamilt::HContainer<TR>* DensityMatrix<TK, TR>::get_DMR_pointer(const int ispin)
     return this->_DMR[ispin - 1];
 }
 
-// get _DMR pointer vector
-template <typename TK, typename TR>
-std::vector<hamilt::HContainer<TR>*> DensityMatrix<TK, TR>::get_DMR_vector() const
-{
-    return this->_DMR;
-}
-
 // get _DMK[ik] pointer
 template <typename TK, typename TR>
 TK* DensityMatrix<TK, TR>::get_DMK_pointer(const int ik) const
@@ -382,6 +375,38 @@ int DensityMatrix<TK, TR>::get_DMK_ncol() const
     return this->_paraV->ncol;
 }
 
+template <typename TK, typename TR>
+void DensityMatrix<TK, TR>::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<this->_DMR.size();is++)
+        {
+            _DMR_save[is].resize(nnr);
+        }
+    }
+    // save _DMR to _DMR_save
+    for(int is=0;is<this->_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<nnr;i++)
+        {
+            DMR_save_pointer[i] = DMR_pointer[i];
+        }
+    }
+
+    ModuleBase::timer::tick("DensityMatrix", "save_DMR");
+}
+
 // calculate DMR from DMK using add_element
 template <typename TK, typename TR>
 void DensityMatrix<TK,TR>::cal_DMR_test()