add test in esolver_ks_lcao.cpp, function after_scf(): update exx to test the exx_energy get by rdmft is exactly equal to exx_energy by abacus

Author: JGHan7

source/module_esolver/esolver_ks_lcao.cpp

@@ -1248,6 +1248,73 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
     // rdmft, added by jghan, 2024-10-17
     if ( PARAM.inp.ab_initio_type == "rdmft" )
     {
+#ifdef __EXX
+        if( GlobalC::exx_info.info_global.cal_exx )
+        {
+            UnitCell* ucell = this->rdmft_solver.ucell;
+            bool exx_spacegroup_symmetry = (PARAM.inp.nspin < 4 && ModuleSymmetry::Symmetry::symm_flag == 1);
+
+            // get the current DM
+            elecstate::DensityMatrix<TK, double> DM_test(&this->pv, PARAM.inp.nspin, this->kv.kvec_d, this->rdmft_solver.nk_total);
+            elecstate::cal_dm_psi(&this->pv, this->pelec->wg, *(this->psi), DM_test);
+            DM_test.init_DMR(&GlobalC::GridD, &GlobalC::ucell);
+            DM_test.cal_DMR();
+            
+            std::vector< std::vector<TK> > DM_XC_test = DM_test.get_DMK_vector();
+            ModuleSymmetry::Symmetry_rotation symrot_exx;
+            if( exx_spacegroup_symmetry )
+            {
+                const std::array<int, 3> period = RI_Util::get_Born_vonKarmen_period(this->kv);
+                symrot_exx.find_irreducible_sector(ucell->symm, ucell->atoms, ucell->st,
+                        RI_Util::get_Born_von_Karmen_cells(period), period, ucell->lat);
+                symrot_exx.cal_Ms(this->kv, *ucell, this->pv);
+
+                // get DM_XC of all k points
+                DM_XC_test = symrot_exx.restore_dm(*this->rdmft_solver.kv, DM_XC_test, this->pv); // class vector could be auto resize()
+            }
+
+            std::vector< const std::vector<TK>* > DM_pointer(DM_XC_test.size());
+            for(int ik=0; ik<DM_XC_test.size(); ++ik) DM_pointer[ik] = &(DM_XC_test[ik]);
+
+            // convert format and update exx
+            if (GlobalC::exx_info.info_ri.real_number)
+            {
+                std::vector<std::map<int,std::map<std::pair<int,std::array<int,3>>,RI::Tensor<double>>>>
+                    Ds_d = std::is_same<TK, double>::value //gamma_only_local
+                ? RI_2D_Comm::split_m2D_ktoR<double>(this->kv, DM_pointer, this->pv, PARAM.inp.nspin)
+                : RI_2D_Comm::split_m2D_ktoR<double>(this->kv, DM_pointer, this->pv, PARAM.inp.nspin, exx_spacegroup_symmetry);
+
+                if (exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace)
+                {
+                    this->exx_lri_double->cal_exx_elec(Ds_d, this->pv, &symrot_exx);
+                }
+                else
+                {
+                    this->exx_lri_double->cal_exx_elec(Ds_d, this->pv);
+                }
+
+                this->pelec->f_en.exx = this->exx_lri_double->Eexx;
+            }
+            else
+            {
+                std::vector<std::map<int,std::map<std::pair<int,std::array<int,3>>,RI::Tensor<std::complex<double>>>>>
+                    Ds_c = std::is_same<TK, double>::value //gamma_only_local
+                ? RI_2D_Comm::split_m2D_ktoR<std::complex<double>>(this->kv, DM_pointer, this->pv, PARAM.inp.nspin)
+                : RI_2D_Comm::split_m2D_ktoR<std::complex<double>>(this->kv, DM_pointer, this->pv, PARAM.inp.nspin, exx_spacegroup_symmetry);
+
+                if (exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace)
+                {
+                    this->exx_lri_complex->cal_exx_elec(Ds_c, this->pv, &symrot_exx);
+                }
+                else
+                {
+                    this->exx_lri_complex->cal_exx_elec(Ds_c, this->pv);
+                }
+
+                this->pelec->f_en.exx = this->exx_lri_complex->Eexx;
+            }
+        }
+#endif
         ModuleBase::matrix occ_number_ks(this->pelec->wg);
         for(int ik=0; ik < occ_number_ks.nr; ++ik)
         {