prepare for PR

Author: JGHan7

source/Makefile.Objects

@@ -736,4 +736,3 @@ OBJS_TENSOR=tensor.o\
 
 OBJS_RDMFT=rdmft.o\
         rdmft_tools.o\
-        rdmft_test.o\

source/module_esolver/esolver_ks_lcao.cpp

@@ -57,6 +57,8 @@
 
 // test RDMFT
 #include "module_rdmft/rdmft.h"
+#include "module_elecstate/module_dm/cal_dm_psi.h"  //temp ,delete
+
 #include <iostream>
 
 namespace ModuleESolver
@@ -953,15 +955,15 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
     }
 
     // 2.5) determine whether rdmft needs to get the initial value, added by jghan, 2024-10-25
-    bool one_step_exx = false;
+    int one_step_exx = false;
     bool get_init_value_rdmft = false;
     if( iter == 1 ) get_init_value_rdmft = true; // the case without hybrid functionals
 #ifdef __EXX
     if( GlobalC::exx_info.info_global.cal_exx )
     {
         if( this->conv_esolver ) one_step_exx = true;
-        // the case with hybrid functionals
-        if( one_step_exx && iter==1 ) get_init_value_rdmft = true;
+        // the case with hybrid functionals, calculate rdmft just after updateExx, cal once in  one inner loop
+        if( one_step_exx ) get_init_value_rdmft = true;
         else get_init_value_rdmft = false;
     }
 #endif
@@ -1132,6 +1134,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
             }
 
             this->rdmft_solver.update_elec(occ_number_ks, *(this->psi));
+            std::cout << "\n******\n" << "update elec in rdmft successfully" << "\n******\n" << std::endl;
 
             //initialize the gradients of Etotal on occupation numbers and wfc, and set all elements to 0. 
             ModuleBase::matrix dE_dOccNum(this->pelec->wg.nr, this->pelec->wg.nc, true);
@@ -1248,8 +1251,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
         ModuleBase::matrix occ_number_ks(this->pelec->wg);
         for(int ik=0; ik < occ_number_ks.nr; ++ik)
         {
-            for(int inb=0; inb < occ_number_ks.nc; ++inb) { occ_number_ks(ik, inb) /= this->kv.wk[ik];
-}
+            for(int inb=0; inb < occ_number_ks.nc; ++inb) occ_number_ks(ik, inb) /= this->kv.wk[ik];
         } 
         this->rdmft_solver.update_elec(occ_number_ks, *(this->psi));
 

source/module_rdmft/CMakeLists.txt

@@ -4,7 +4,6 @@ if(ENABLE_LCAO)
         OBJECT
         rdmft.cpp
         rdmft_tools.cpp
-        rdmft_test.cpp
     )
 endif()
 

source/module_rdmft/rdmft.cpp

@@ -193,6 +193,15 @@ void RDMFT<TK, TR>::init(Gint_Gamma& GG_in, Gint_k& GK_in, Parallel_Orbitals& Pa
 #ifdef __EXX
     if( GlobalC::exx_info.info_global.cal_exx )
     {
+        exx_spacegroup_symmetry = (PARAM.inp.nspin < 4 && ModuleSymmetry::Symmetry::symm_flag == 1);
+        if (exx_spacegroup_symmetry)
+        {
+            const std::array<int, 3>& period = RI_Util::get_Born_vonKarmen_period(*kv);
+            this->symrot_exx.find_irreducible_sector(ucell->symm, ucell->atoms, ucell->st,
+                    RI_Util::get_Born_von_Karmen_cells(period), period, ucell->lat);
+            this->symrot_exx.cal_Ms(*kv, *ucell, *ParaV);
+        }
+
         if (GlobalC::exx_info.info_ri.real_number)
         {
             Vxc_fromRI_d = new Exx_LRI<double>(GlobalC::exx_info.info_ri);
@@ -571,16 +580,6 @@ void RDMFT<TK, TR>::cal_V_hartree()
 template <typename TK, typename TR>
 void RDMFT<TK, TR>::cal_V_XC()
 {
-    HR_dft_XC->set_zero();
-    HR_exx_XC->set_zero();
-
-    std::vector< std::vector<TK> > DM_XC(nk_total, std::vector<TK>(ParaV->nloc));
-    std::vector< const std::vector<TK>* > DM_XC_pointer(nk_total);
-    for(int ik=0; ik<nk_total; ++ik) { DM_XC_pointer[ik] = &DM_XC[ik];
-}
-    
-    get_DM_XC(DM_XC);
-
     // // //test
     // DM_XC_pass = DM_XC;
 
@@ -609,6 +608,7 @@ void RDMFT<TK, TR>::cal_V_XC()
 
     if( !only_exx_type )
     {
+        HR_dft_XC->set_zero();
         if( PARAM.inp.gamma_only )
         {
             // this can be optimized, use potXC.update_from_charge()
@@ -661,17 +661,35 @@ void RDMFT<TK, TR>::cal_V_XC()
 #ifdef __EXX
     if(GlobalC::exx_info.info_global.cal_exx)
     {
+        HR_exx_XC->set_zero();
+
+        std::vector< std::vector<TK> > DM_XC(nk_total, std::vector<TK>(ParaV->nloc));
+        get_DM_XC(DM_XC);
+        // get DM_XC of all k points
+        if( exx_spacegroup_symmetry )
+        {
+            DM_XC = symrot_exx.restore_dm(*this->kv, DM_XC, *ParaV); // class vector could be auto resize()
+        }
+        std::vector< const std::vector<TK>* > DM_XC_pointer(DM_XC.size());
+        for(int ik=0; ik<DM_XC.size(); ++ik) DM_XC_pointer[ik] = &DM_XC[ik];
+
         if (GlobalC::exx_info.info_ri.real_number)
         {
             // transfer the DM_XC to appropriate format
-            std::vector<std::map<int,std::map<std::pair<int,std::array<int,3>>,RI::Tensor<double>>>> Ds_XC_d = 
-                RI_2D_Comm::split_m2D_ktoR<double>(*kv, DM_XC_pointer, *ParaV, nspin);
+            std::vector<std::map<int,std::map<std::pair<int,std::array<int,3>>,RI::Tensor<double>>>> 
+                Ds_XC_d = std::is_same<TK, double>::value //gamma_only_local
+                ? RI_2D_Comm::split_m2D_ktoR<double>(*kv, DM_XC_pointer, *ParaV, nspin)
+                : RI_2D_Comm::split_m2D_ktoR<double>(*kv, DM_XC_pointer, *ParaV, nspin, this->exx_spacegroup_symmetry);
 
             // provide the Ds_XC to Vxc_fromRI(V_exx_XC)
-            // Vxc_fromRI_d = new Exx_LRI<double>(GlobalC::exx_info.info_ri);
-            // Vxc_fromRI_d->init(MPI_COMM_WORLD, *kv);
-            // Vxc_fromRI_d->cal_exx_ions();
-            Vxc_fromRI_d->cal_exx_elec(Ds_XC_d, *ParaV);
+            if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace)
+            {
+                Vxc_fromRI_d->cal_exx_elec(Ds_XC_d, *ParaV, &this->symrot_exx);
+            }
+            else
+            {
+                Vxc_fromRI_d->cal_exx_elec(Ds_XC_d, *ParaV);
+            }
 
             // when we doing V_exx_XC.contributeHk(ik), we get HK_XC constructed by the special DM_XC
             V_exx_XC = new hamilt::OperatorEXX<hamilt::OperatorLCAO<TK, TR>>(
@@ -682,24 +700,24 @@ void RDMFT<TK, TR>::cal_V_XC()
                 nullptr,
                 hamilt::Add_Hexx_Type::k
             );
-            // V_exx_XC = new hamilt::OperatorEXX<hamilt::OperatorLCAO<TK, TR>>(
-            //     hsk_exx_XC,
-            //     HR_exx_XC,
-            //     *kv
-            // );
-
         }
         else
         {
             // transfer the DM_XC to appropriate format
-            std::vector<std::map<int,std::map<std::pair<int,std::array<int,3>>,RI::Tensor<std::complex<double>>>>> Ds_XC_c = 
-                RI_2D_Comm::split_m2D_ktoR<std::complex<double>>(*kv, DM_XC_pointer, *ParaV, nspin);
-
-            // provide the Ds_XC to Vxc_fromRI(V_exx_XC)
-            // Vxc_fromRI_c = new Exx_LRI<std::complex<double>>(GlobalC::exx_info.info_ri);
-            // Vxc_fromRI_c->init(MPI_COMM_WORLD, *kv);
-            // Vxc_fromRI_c->cal_exx_ions();
-            Vxc_fromRI_c->cal_exx_elec(Ds_XC_c, *ParaV);
+            std::vector<std::map<int,std::map<std::pair<int,std::array<int,3>>,RI::Tensor<std::complex<double>>>>> 
+                Ds_XC_c = std::is_same<TK, double>::value //gamma_only_local
+                ? RI_2D_Comm::split_m2D_ktoR<std::complex<double>>(*kv, DM_XC_pointer, *ParaV, nspin)
+                : RI_2D_Comm::split_m2D_ktoR<std::complex<double>>(*kv, DM_XC_pointer, *ParaV, nspin, this->exx_spacegroup_symmetry);
+
+            // // provide the Ds_XC to Vxc_fromRI(V_exx_XC)
+            if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace)
+            {
+                Vxc_fromRI_c->cal_exx_elec(Ds_XC_c, *ParaV, &this->symrot_exx);
+            }
+            else
+            {
+                Vxc_fromRI_c->cal_exx_elec(Ds_XC_c, *ParaV);
+            }
 
             // when we doing V_exx_XC.contributeHk(ik), we get HK_XC constructed by the special DM_XC
             V_exx_XC = new hamilt::OperatorEXX<hamilt::OperatorLCAO<TK, TR>>(
@@ -840,6 +858,8 @@ void RDMFT<TK, TR>::cal_Energy(const int cal_type)
     double E_deband_KS = pelec->f_en.deband;
     double E_deband_harris_KS = pelec->f_en.deband_harris;
 
+    double E_exxType_rdmft = 0.0; // delete in the future
+
     if( cal_type == 1 )
     {
         // for E_TV
@@ -862,6 +882,7 @@ void RDMFT<TK, TR>::cal_Energy(const int cal_type)
             occNum_Mul_wfcHwfc(wk_fun_occNum, wfcHwfc_exx_XC, Exc_n_k, 1);
             E_RDMFT[2] = getEnergy(Exc_n_k);
             Parallel_Reduce::reduce_all(E_RDMFT[2]);
+            E_exxType_rdmft = E_RDMFT[2];
         }
 #endif
         E_RDMFT[2] += etxc;
@@ -934,7 +955,8 @@ void RDMFT<TK, TR>::cal_Energy(const int cal_type)
                 << "\nE_descf:         " << E_descf 
                 << "\n\nEtotal_RDMFT:    " << Etotal 
                 << "\n\nExc_ksdft:       " << E_xc_KS 
-                << "\nE_exx_ksdft:     " << E_exx_KS 
+                << "\nE_exx_ksdft:     " << E_exx_KS
+                << "\nE_exxType_rdmft: " << E_exxType_rdmft
                 <<"\n******\n" << std::endl;
     }
     std::cout << std::defaultfloat;

source/module_rdmft/rdmft.h

@@ -41,6 +41,7 @@
 #include "module_ri/Exx_LRI.h"
 #include "module_ri/RI_2D_Comm.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h"
+#include "module_ri/module_exx_symmetry/symmetry_rotation.h"
 #endif
 
 // there are some operator reload to print data in different formats
@@ -179,6 +180,8 @@ class RDMFT
 #ifdef __EXX
     Exx_LRI<double>* Vxc_fromRI_d = nullptr;
     Exx_LRI<std::complex<double>>* Vxc_fromRI_c = nullptr;
+    ModuleSymmetry::Symmetry_rotation symrot_exx;
+    bool exx_spacegroup_symmetry = false;
 #endif
 
     double Etotal = 0.0;