Refactor: simplify Mixing DM

source/Makefile.Objects

@@ -677,7 +677,6 @@ OBJS_MODULE_RI=conv_coulomb_pot_k.o\
       Matrix_Orbs22.o\
       RI_2D_Comm.o\
       Mix_DMk_2D.o\
-      Mix_Matrix.o\
       symmetry_rotation.o\
       symmetry_rotation_output.o\
       symmetry_irreducible_sector.o\

source/source_lcao/module_ri/CMakeLists.txt

@@ -8,7 +8,6 @@ if (ENABLE_LIBRI)
         Matrix_Orbs22.cpp
         RI_2D_Comm.cpp
         Mix_DMk_2D.cpp
-        Mix_Matrix.cpp
   )
 
   if(ENABLE_LCAO)

source/source_lcao/module_ri/Exx_LRI_interface.h

@@ -124,7 +124,8 @@ class Exx_LRI_Interface
     std::shared_ptr<Exx_LRI<Tdata>> exx_ptr;
 
 private:
-    Mix_DMk_2D mix_DMk_2D;
+
+    Mix_DMk_2D<T> mix_DMk_2D;
 
     bool exx_spacegroup_symmetry = false;
     ModuleSymmetry::Symmetry_rotation symrot_;

source/source_lcao/module_ri/Exx_LRI_interface.hpp

@@ -35,9 +35,7 @@ void Exx_LRI_Interface<T, Tdata>::cal_exx_ions(const UnitCell& ucell, const bool
     ModuleBase::TITLE("Exx_LRI_Interface","cal_exx_ions");
     if(!this->flag_finish.init)
         { throw std::runtime_error("Exx init unfinished when "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); }
-
     this->exx_ptr->cal_exx_ions(ucell, write_cv);
-
     this->flag_finish.ions = true;
 }
 
@@ -49,13 +47,9 @@ void Exx_LRI_Interface<T, Tdata>::cal_exx_elec(const std::vector<std::map<TA, st
 {
     ModuleBase::TITLE("Exx_LRI_Interface","cal_exx_elec");
     if(!this->flag_finish.init || !this->flag_finish.ions)
-    { 
-        throw std::runtime_error("Exx init unfinished when "
-        +std::string(__FILE__)+" line "+std::to_string(__LINE__)); 
-    }
+        { throw std::runtime_error("Exx init unfinished when "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); }
 
     this->exx_ptr->cal_exx_elec(Ds, ucell, pv, p_symrot);
-
     this->flag_finish.elec = true;
 }
 
@@ -64,17 +58,11 @@ void Exx_LRI_Interface<T, Tdata>::cal_exx_force(const int& nat)
 {
     ModuleBase::TITLE("Exx_LRI_Interface","cal_exx_force");
     if(!this->flag_finish.init || !this->flag_finish.ions)
-    { 
-        throw std::runtime_error("Exx init unfinished when "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); 
-    }
+        { throw std::runtime_error("Exx init unfinished when "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); }
     if(!this->flag_finish.elec)
-    { 
-        throw std::runtime_error("Exx Hamiltonian unfinished when "+std::string(__FILE__)
-        +" line "+std::to_string(__LINE__)); 
-    }
+        { throw std::runtime_error("Exx Hamiltonian unfinished when "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); }
 
     this->exx_ptr->cal_exx_force(nat);
-
     this->flag_finish.force = true;
 }
 
@@ -83,24 +71,17 @@ void Exx_LRI_Interface<T, Tdata>::cal_exx_stress(const double& omega, const doub
 {
     ModuleBase::TITLE("Exx_LRI_Interface","cal_exx_stress");
     if(!this->flag_finish.init || !this->flag_finish.ions)
-    { 
-        throw std::runtime_error("Exx init unfinished when "
-                +std::string(__FILE__)+" line "+std::to_string(__LINE__)); 
-    }
+        { throw std::runtime_error("Exx init unfinished when "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); }
     if(!this->flag_finish.elec)
-    { 
-        throw std::runtime_error("Exx Hamiltonian unfinished when "
-                +std::string(__FILE__)+" line "+std::to_string(__LINE__)); 
-    }
+        { throw std::runtime_error("Exx Hamiltonian unfinished when "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); }
 
     this->exx_ptr->cal_exx_stress(omega, lat0);
-
     this->flag_finish.stress = true;
 }
 
 template<typename T, typename Tdata>
 void Exx_LRI_Interface<T, Tdata>::exx_before_all_runners(
-    const K_Vectors& kv, 
+    const K_Vectors& kv,
     const UnitCell& ucell,
     const Parallel_2D& pv)
 {
@@ -147,14 +128,15 @@ void Exx_LRI_Interface<T, Tdata>::exx_beforescf(const int istep,
     if(GlobalC::exx_info.info_global.cal_exx)
     {
         if (this->exx_spacegroup_symmetry)
-            {this->mix_DMk_2D.set_nks(kv.get_nkstot_full() * (PARAM.inp.nspin == 2 ? 2 : 1), PARAM.globalv.gamma_only_local);}
+            { this->mix_DMk_2D.set_nks(kv.get_nkstot_full() * (PARAM.inp.nspin == 2 ? 2 : 1)); }
         else
-            {this->mix_DMk_2D.set_nks(kv.get_nks(), PARAM.globalv.gamma_only_local);}
+            { this->mix_DMk_2D.set_nks(kv.get_nks()); }
 
-        if(GlobalC::exx_info.info_global.separate_loop)
-            { this->mix_DMk_2D.set_mixing(nullptr); }
+        if (GlobalC::exx_info.info_global.separate_loop)
+            { this->mix_DMk_2D.set_mixing_plain(GlobalC::exx_info.info_global.mixing_beta_for_loop1); }
         else
             { this->mix_DMk_2D.set_mixing(chgmix.get_mixing()); }
+
         // for exx two_level scf
         this->two_level_step = 0;
     }
@@ -171,13 +153,13 @@ void Exx_LRI_Interface<T, Tdata>::exx_eachiterinit(const int istep,
     ModuleBase::TITLE("Exx_LRI_Interface","exx_eachiterinit");
     if (GlobalC::exx_info.info_global.cal_exx)
     {
-        if (!GlobalC::exx_info.info_global.separate_loop 
-            && (this->two_level_step 
-                || istep > 0 
+        if (!GlobalC::exx_info.info_global.separate_loop
+            && (this->two_level_step
+                || istep > 0
                 || PARAM.inp.init_wfc == "file") // non separate loop case
-            || (GlobalC::exx_info.info_global.separate_loop 
-                && PARAM.inp.init_wfc == "file" 
-                && this->two_level_step == 0 
+            || (GlobalC::exx_info.info_global.separate_loop
+                && PARAM.inp.init_wfc == "file"
+                && this->two_level_step == 0
                 && iter == 1)
            )  // the first iter in separate loop case
         {
@@ -186,33 +168,25 @@ void Exx_LRI_Interface<T, Tdata>::exx_eachiterinit(const int istep,
             auto cal = [this, &ucell,&kv, &flag_restart](const elecstate::DensityMatrix<T, double>& dm_in)
             {
                 if (this->exx_spacegroup_symmetry)
-                    { this->mix_DMk_2D.mix(symrot_.restore_dm(kv,dm_in.get_DMK_vector(), *dm_in.get_paraV_pointer()), flag_restart); }
+                    { this->mix_DMk_2D.mix(symrot_.restore_dm(kv, dm_in.get_DMK_vector(), *dm_in.get_paraV_pointer()), flag_restart); }
                 else
                     { this->mix_DMk_2D.mix(dm_in.get_DMK_vector(), flag_restart); }
-                const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>>
-                    Ds = PARAM.globalv.gamma_only_local
-                        ? RI_2D_Comm::split_m2D_ktoR<Tdata>(
-                            ucell,
-                            *this->exx_ptr->p_kv,
-                            this->mix_DMk_2D.get_DMk_gamma_out(),
-                            *dm_in.get_paraV_pointer(),
-                            PARAM.inp.nspin)
-                        : RI_2D_Comm::split_m2D_ktoR<Tdata>(
-                            ucell,
-                            *this->exx_ptr->p_kv,
-                            this->mix_DMk_2D.get_DMk_k_out(),
-                            *dm_in.get_paraV_pointer(),
-                            PARAM.inp.nspin,
-                            this->exx_spacegroup_symmetry);
-
-                if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_ri.exx_symmetry_realspace)
+                const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>> Ds =
+                    RI_2D_Comm::split_m2D_ktoR<Tdata>(
+                        ucell,
+                        *this->exx_ptr->p_kv,
+                        this->mix_DMk_2D.get_DMk_out(),
+                        *dm_in.get_paraV_pointer(),
+                        PARAM.inp.nspin,
+                        this->exx_spacegroup_symmetry);
+                if(this->exx_spacegroup_symmetry && GlobalC::exx_info.info_ri.exx_symmetry_realspace)
                     { this->cal_exx_elec(Ds, ucell,*dm_in.get_paraV_pointer(), &this->symrot_); }
                 else
                     { this->cal_exx_elec(Ds, ucell,*dm_in.get_paraV_pointer()); }
             };
 
             if(istep > 0 && flag_restart)
-                { cal(*dm_last_step); }
+                { cal(*this->dm_last_step); }
             else
                 { cal(dm); }
         }
@@ -387,21 +361,23 @@ bool Exx_LRI_Interface<T, Tdata>::exx_after_converge(
             // if init_wfc == "file", DM is calculated in the 1st iter of the 1st two-level step, so we mix it here
             const bool flag_restart = (this->two_level_step == 0 && PARAM.inp.init_wfc != "file") ? true : false;
 
-            if (this->exx_spacegroup_symmetry)
-                {this->mix_DMk_2D.mix(symrot_.restore_dm(kv, dm.get_DMK_vector(), *dm.get_paraV_pointer()), flag_restart);}
+            if(this->exx_spacegroup_symmetry)
+                { this->mix_DMk_2D.mix(symrot_.restore_dm(kv, dm.get_DMK_vector(), *dm.get_paraV_pointer()), flag_restart); }
             else
-                {this->mix_DMk_2D.mix(dm.get_DMK_vector(), flag_restart);}
-
-            // GlobalC::exx_lcao.cal_exx_elec(p_esolver->LOC, p_esolver->LOWF.wfc_k_grid);
-            const std::vector<std::map<int, std::map<std::pair<int, std::array<int, 3>>, RI::Tensor<Tdata>>>>
-                Ds = std::is_same<T, double>::value //gamma_only_local
-                ? RI_2D_Comm::split_m2D_ktoR<Tdata>(ucell,*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_gamma_out(), *dm.get_paraV_pointer(), nspin)
-                : RI_2D_Comm::split_m2D_ktoR<Tdata>(ucell,*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_k_out(), *dm.get_paraV_pointer(), nspin, this->exx_spacegroup_symmetry);
-
-            if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_ri.exx_symmetry_realspace)
+                { this->mix_DMk_2D.mix(dm.get_DMK_vector(), flag_restart); }
+            const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>> Ds =
+                RI_2D_Comm::split_m2D_ktoR<Tdata>(
+                    ucell,
+                    *this->exx_ptr->p_kv,
+                    this->mix_DMk_2D.get_DMk_out(),
+                    *dm.get_paraV_pointer(),
+                    nspin,
+                    this->exx_spacegroup_symmetry);
+            if(this->exx_spacegroup_symmetry && GlobalC::exx_info.info_ri.exx_symmetry_realspace)
                 { this->cal_exx_elec(Ds, ucell, *dm.get_paraV_pointer(), &this->symrot_); }
             else
                 { this->cal_exx_elec(Ds, ucell, *dm.get_paraV_pointer()); }    // restore DM but not Hexx
+
             iter = 0;
             this->two_level_step++;
 

source/source_lcao/module_ri/Mix_DMk_2D.cpp

@@ -4,69 +4,85 @@
 //=======================
 
 #include "Mix_DMk_2D.h"
+#include "source_base/module_mixing/plain_mixing.h"
 #include "source_base/tool_title.h"
 
-Mix_DMk_2D &Mix_DMk_2D::set_nks(const int nks, const bool gamma_only_in)
+#include <cassert>
+
+template <typename Tdata>
+Mix_DMk_2D<Tdata>::~Mix_DMk_2D<Tdata>()
+{
+    if(this->flag_del_mixing)
+        delete this->mixing;
+}
+
+template <typename Tdata>
+void Mix_DMk_2D<Tdata>::set_nks(const int nks)
 {
-    ModuleBase::TITLE("Mix_DMk_2D", "set_nks");
-    this->gamma_only = gamma_only_in;
-    if (this->gamma_only)
-		this->mix_DMk_gamma.resize(nks);
-	else
-		this->mix_DMk_k.resize(nks);
-	return *this;
+    this->mix_DMk.clear();
+    this->mix_DMk.resize(nks);
 }
 
-Mix_DMk_2D &Mix_DMk_2D::set_mixing(Base_Mixing::Mixing* mixing_in)
+template <typename Tdata>
+void Mix_DMk_2D<Tdata>::set_mixing(Base_Mixing::Mixing* mixing_in)
 {
-	ModuleBase::TITLE("Mix_DMk_2D","set_mixing");
-	if(this->gamma_only)
-        for (Mix_Matrix<std::vector<double>>& mix_one : this->mix_DMk_gamma)
-			mix_one.init(mixing_in);
-	else
-        for (Mix_Matrix<std::vector<std::complex<double>>>& mix_one : this->mix_DMk_k)
-			mix_one.init(mixing_in);
-	return *this;
+    if(this->flag_del_mixing)
+        delete this->mixing;
+    this->mixing = mixing_in;
+    this->flag_del_mixing = false;
 }
 
-Mix_DMk_2D &Mix_DMk_2D::set_mixing_beta(const double mixing_beta)
+template <typename Tdata>
+void Mix_DMk_2D<Tdata>::set_mixing_plain(const double& mixing_beta)
 {
-	ModuleBase::TITLE("Mix_DMk_2D","set_mixing_beta");
-	if(this->gamma_only)
-        for (Mix_Matrix<std::vector<double>>& mix_one : this->mix_DMk_gamma)
-			mix_one.mixing_beta = mixing_beta;
-	else
-        for (Mix_Matrix<std::vector<std::complex<double>>>& mix_one : this->mix_DMk_k)
-			mix_one.mixing_beta = mixing_beta;
-	return *this;
+    if(this->flag_del_mixing)
+        delete this->mixing;
+    this->mixing = new Base_Mixing::Plain_Mixing(mixing_beta);
+    this->flag_del_mixing = true;
 }
 
-void Mix_DMk_2D::mix(const std::vector<std::vector<double>>& dm, const bool flag_restart)
+template <typename Tdata>
+void Mix_DMk_2D<Tdata>::mix(const std::vector<std::vector<Tdata>>& dm, const bool flag_restart)
 {
-	ModuleBase::TITLE("Mix_DMk_2D","mix");
-	assert(this->mix_DMk_gamma.size() == dm.size());
-	for(int ik=0; ik<dm.size(); ++ik)
-		this->mix_DMk_gamma[ik].mix(dm[ik], flag_restart);
+    ModuleBase::TITLE("Mix_DMk_2D", "mix");
+    if (flag_restart)
+        { this->restart_all(dm); }
+    else
+        { this->mix_all(dm); }
 }
-void Mix_DMk_2D::mix(const std::vector<std::vector<std::complex<double>>>& dm, const bool flag_restart)
+
+template <typename Tdata>
+std::vector<const std::vector<Tdata>*> Mix_DMk_2D<Tdata>::get_DMk_out() const
 {
-	ModuleBase::TITLE("Mix_DMk_2D","mix");
-	assert(this->mix_DMk_k.size() == dm.size());
-	for(int ik=0; ik<dm.size(); ++ik)
-		this->mix_DMk_k[ik].mix(dm[ik], flag_restart);
+    std::vector<const std::vector<Tdata>*> DMk_out(this->mix_DMk.size());
+    for (int ik = 0; ik < this->mix_DMk.size(); ++ik)
+        { DMk_out[ik] = &this->mix_DMk[ik].data_out; }
+    return DMk_out;
 }
 
-std::vector<const std::vector<double>*> Mix_DMk_2D::get_DMk_gamma_out() const
+template <typename Tdata>
+void Mix_DMk_2D<Tdata>::restart_all(const std::vector<std::vector<Tdata>>& data_in)
 {
-    std::vector<const std::vector<double>*> DMk_out(this->mix_DMk_gamma.size());
-	for(int ik=0; ik<this->mix_DMk_gamma.size(); ++ik)
-		DMk_out[ik] = &this->mix_DMk_gamma[ik].get_data_out();
-	return DMk_out;
+    assert(this->mix_DMk.size() == data_in.size());
+    assert(this->mixing != nullptr);
+    for (int ik = 0; ik < data_in.size(); ++ik)
+    {
+        this->mix_DMk[ik].data_out = data_in[ik];
+        this->mixing->init_mixing_data(this->mix_DMk[ik].mixing_data, data_in[ik].size(), sizeof(Tdata));
+    }
 }
-std::vector<const std::vector<std::complex<double>>*> Mix_DMk_2D::get_DMk_k_out() const
+
+template <typename Tdata>
+void Mix_DMk_2D<Tdata>::mix_all(const std::vector<std::vector<Tdata>>& data_in)
 {
-    std::vector<const std::vector<std::complex<double>>*> DMk_out(this->mix_DMk_k.size());
-	for(int ik=0; ik<this->mix_DMk_k.size(); ++ik)
-		DMk_out[ik] = &this->mix_DMk_k[ik].get_data_out();
-	return DMk_out;
-}
+    assert(this->mix_DMk.size() == data_in.size());
+    assert(this->mixing != nullptr);
+    for (int ik = 0; ik < data_in.size(); ++ik)
+    {
+        this->mixing->push_data(this->mix_DMk[ik].mixing_data, this->mix_DMk[ik].data_out.data(), data_in[ik].data(), nullptr, false);
+        this->mixing->mix_data(this->mix_DMk[ik].mixing_data, this->mix_DMk[ik].data_out.data());
+    }
+}
+
+template class Mix_DMk_2D<double>;
+template class Mix_DMk_2D<std::complex<double>>;