Refactor:remove GlobalC::ucell in module_ri

source/module_esolver/esolver_ks_lcao.cpp

@@ -178,12 +178,12 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
             // initialize 2-center radial tables for EXX-LRI
             if (GlobalC::exx_info.info_ri.real_number)
             {
-                this->exx_lri_double->init(MPI_COMM_WORLD, this->kv, orb_);
+                this->exx_lri_double->init(MPI_COMM_WORLD, ucell,this->kv, orb_);
                 this->exd->exx_before_all_runners(this->kv, ucell, this->pv);
             }
             else
             {
-                this->exx_lri_complex->init(MPI_COMM_WORLD, this->kv, orb_);
+                this->exx_lri_complex->init(MPI_COMM_WORLD, ucell,this->kv, orb_);
                 this->exc->exx_before_all_runners(this->kv, ucell, this->pv);
             }
         }
@@ -623,13 +623,15 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
         if (GlobalC::exx_info.info_ri.real_number)
         {
             this->exd->exx_eachiterinit(istep,
+                                        ucell,
                                         *dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
                                         this->kv,
                                         iter);
         }
         else
         {
             this->exc->exx_eachiterinit(istep,
+                                        ucell,
                                         *dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
                                         this->kv,
                                         iter);
@@ -1052,7 +1054,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
                 occ_number_ks(ik, inb) /= this->kv.wk[ik];
             }
         }
-        this->rdmft_solver.update_elec(occ_number_ks, *(this->psi));
+        this->rdmft_solver.update_elec(ucell,occ_number_ks, *(this->psi));
 
         //! initialize the gradients of Etotal with respect to occupation numbers and wfc, 
         //! and set all elements to 0. 
@@ -1074,10 +1076,11 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
         RPA_LRI<TK, double> rpa_lri_double(GlobalC::exx_info.info_ri);
         rpa_lri_double.cal_postSCF_exx(*dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
                                        MPI_COMM_WORLD,
+                                       ucell,
                                        this->kv,
                                        orb_);
         rpa_lri_double.init(MPI_COMM_WORLD, this->kv, orb_.cutoffs());
-        rpa_lri_double.out_for_RPA(this->pv, *(this->psi), this->pelec);
+        rpa_lri_double.out_for_RPA(ucell,this->pv, *(this->psi), this->pelec);
     }
 #endif
 

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp

@@ -375,25 +375,25 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
         {
             if (GlobalC::exx_info.info_ri.real_number)
             {
-                exx_lri_double.cal_exx_force();
+                exx_lri_double.cal_exx_force(ucell.nat);
                 force_exx = GlobalC::exx_info.info_global.hybrid_alpha * exx_lri_double.force_exx;
             }
             else
             {
-                exx_lri_complex.cal_exx_force();
+                exx_lri_complex.cal_exx_force(ucell.nat);
                 force_exx = GlobalC::exx_info.info_global.hybrid_alpha * exx_lri_complex.force_exx;
             }
         }
         if (isstress)
         {
             if (GlobalC::exx_info.info_ri.real_number)
             {
-                exx_lri_double.cal_exx_stress();
+                exx_lri_double.cal_exx_stress(ucell.omega,ucell.lat0);
                 stress_exx = GlobalC::exx_info.info_global.hybrid_alpha * exx_lri_double.stress_exx;
             }
             else
             {
-                exx_lri_complex.cal_exx_stress();
+                exx_lri_complex.cal_exx_stress(ucell.omega,ucell.lat0);
                 stress_exx = GlobalC::exx_info.info_global.hybrid_alpha * exx_lri_complex.stress_exx;
             }
         }

source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp

@@ -21,7 +21,7 @@
 // Peize Lin add 2022.09.13
 
 template <typename Tdata>
-void sparse_format::cal_HR_exx(
+void sparse_format::cal_HR_exx(const UnitCell& ucell,
     const Parallel_Orbitals& pv,
     LCAO_HS_Arrays& HS_Arrays,
     const int& current_spin,
@@ -36,15 +36,15 @@ void sparse_format::cal_HR_exx(
     const Tdata frac = GlobalC::exx_info.info_global.hybrid_alpha;
 
     std::map<int, std::array<double, 3>> atoms_pos;
-    for (int iat = 0; iat < GlobalC::ucell.nat; ++iat) {
+    for (int iat = 0; iat < ucell.nat; ++iat) {
         atoms_pos[iat] = RI_Util::Vector3_to_array3(
-            GlobalC::ucell.atoms[GlobalC::ucell.iat2it[iat]]
-                .tau[GlobalC::ucell.iat2ia[iat]]);
+            ucell.atoms[ucell.iat2it[iat]]
+                .tau[ucell.iat2ia[iat]]);
     }
     const std::array<std::array<double, 3>, 3> latvec
-        = {RI_Util::Vector3_to_array3(GlobalC::ucell.a1), // too bad to use GlobalC here, 
-           RI_Util::Vector3_to_array3(GlobalC::ucell.a2),
-           RI_Util::Vector3_to_array3(GlobalC::ucell.a3)};
+        = {RI_Util::Vector3_to_array3(ucell.a1), // too bad to use GlobalC here, 
+           RI_Util::Vector3_to_array3(ucell.a2),
+           RI_Util::Vector3_to_array3(ucell.a3)};
 
     const std::array<int, 3> Rs_period = {nmp[0], nmp[1], nmp[2]};
 
@@ -84,7 +84,7 @@ void sparse_format::cal_HR_exx(
 
                 for (size_t iw0 = 0; iw0 < Hexx.shape[0]; ++iw0) 
                 {
-                    const int iwt0 = RI_2D_Comm::get_iwt(iat0, iw0, is0_b);
+                    const int iwt0 = RI_2D_Comm::get_iwt(ucell,iat0, iw0, is0_b);
                     const int iwt0_local = pv.global2local_row(iwt0);
 
                     if (iwt0_local < 0) 
@@ -94,7 +94,7 @@ void sparse_format::cal_HR_exx(
 
                     for (size_t iw1 = 0; iw1 < Hexx.shape[1]; ++iw1) 
                     {
-                        const int iwt1 = RI_2D_Comm::get_iwt(iat1, iw1, is1_b);
+                        const int iwt1 = RI_2D_Comm::get_iwt(ucell,iat1, iw1, is1_b);
                         const int iwt1_local = pv.global2local_col(iwt1);
 
                         if (iwt1_local < 0) 

source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp

@@ -403,6 +403,7 @@ HamiltLCAO<TK, TR>::HamiltLCAO(Gint_Gamma* GG_in,
         // and calculate Cs, Vs
         Operator<TK>* exx = new OperatorEXX<OperatorLCAO<TK, TR>>(this->hsk,
                                                                   this->hR,
+                                                                  ucell,
                                                                   *this->kv,
                                                                   Hexxd,
                                                                   Hexxc,

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h

@@ -28,6 +28,7 @@ class OperatorEXX<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
 public:
     OperatorEXX<OperatorLCAO<TK, TR>>(HS_Matrix_K<TK>* hsk_in,
         hamilt::HContainer<TR>* hR_in,
+        const UnitCell& ucell,
         const K_Vectors& kv_in,
         std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd_in = nullptr,
         std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc_in = nullptr,
@@ -58,7 +59,9 @@ class OperatorEXX<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
       const int istep = 0; // the ion step
 
       void add_loaded_Hexx(const int ik);
-
+
+      const UnitCell& ucell;
+
       const K_Vectors& kv;
 
       // if k points has no shift, use cell_nearest to reduce the memory cost

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp

@@ -49,9 +49,9 @@ namespace hamilt
         auto* pv = hR->get_paraV();
         auto Rs = RI_Util::get_Born_von_Karmen_cells(Rs_period);
         bool need_allocate = false;
-        for (int iat0 = 0;iat0 < GlobalC::ucell.nat;++iat0)
+        for (int iat0 = 0;iat0 < nat;++iat0)
         {
-            for (int iat1 = 0;iat1 < GlobalC::ucell.nat;++iat1)
+            for (int iat1 = 0;iat1 < nat;++iat1)
             {
                 // complete the atom pairs that has orbitals in this processor but not in hR due to the adj_list 
                 // but adj_list is not enought for EXX, which is more nonlocal than Nonlocal 
@@ -81,6 +81,7 @@ namespace hamilt
 template <typename TK, typename TR>
 OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
     HContainer<TR>*hR_in,
+    const UnitCell& ucell_in,
 	const K_Vectors& kv_in,
 	std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd_in,
 	std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc_in,
@@ -89,6 +90,7 @@ OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
     int* two_level_step_in,
 	const bool restart_in)
     : OperatorLCAO<TK, TR>(hsk_in, kv_in.kvec_d, hR_in),
+    ucell(ucell_in),
     kv(kv_in),
     Hexxd(Hexxd_in),
     Hexxc(Hexxc_in),
@@ -115,12 +117,12 @@ OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
             // Read HexxR in CSR format
             if (GlobalC::exx_info.info_ri.real_number)
             {
-                ModuleIO::read_Hexxs_csr(file_name_exx, GlobalC::ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd);
+                ModuleIO::read_Hexxs_csr(file_name_exx, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd);
                 if (this->add_hexx_type == Add_Hexx_Type::R) { reallocate_hcontainer(*Hexxd, this->hR); }
             }
             else
             {
-                ModuleIO::read_Hexxs_csr(file_name_exx, GlobalC::ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc);
+                ModuleIO::read_Hexxs_csr(file_name_exx, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc);
                 if (this->add_hexx_type == Add_Hexx_Type::R) { reallocate_hcontainer(*Hexxc, this->hR); }
             }
         }
@@ -154,19 +156,19 @@ OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
             {
                 // set cell_nearest
                 std::map<int, std::array<double, 3>> atoms_pos;
-                for (int iat = 0; iat < GlobalC::ucell.nat; ++iat) {
+                for (int iat = 0; iat < ucell.nat; ++iat) {
                     atoms_pos[iat] = RI_Util::Vector3_to_array3(
-                        GlobalC::ucell.atoms[GlobalC::ucell.iat2it[iat]]
-                        .tau[GlobalC::ucell.iat2ia[iat]]);
+                        ucell.atoms[ucell.iat2it[iat]]
+                        .tau[ucell.iat2ia[iat]]);
                 }
                 const std::array<std::array<double, 3>, 3> latvec
-                    = { RI_Util::Vector3_to_array3(GlobalC::ucell.a1),
-                       RI_Util::Vector3_to_array3(GlobalC::ucell.a2),
-                       RI_Util::Vector3_to_array3(GlobalC::ucell.a3) };
+                    = { RI_Util::Vector3_to_array3(ucell.a1),
+                       RI_Util::Vector3_to_array3(ucell.a2),
+                       RI_Util::Vector3_to_array3(ucell.a3) };
                 this->cell_nearest.init(atoms_pos, latvec, Rs_period);
-                reallocate_hcontainer(GlobalC::ucell.nat, this->hR, Rs_period, &this->cell_nearest);
+                reallocate_hcontainer(ucell.nat, this->hR, Rs_period, &this->cell_nearest);
             }
-            else { reallocate_hcontainer(GlobalC::ucell.nat, this->hR, Rs_period); }
+            else { reallocate_hcontainer(ucell.nat, this->hR, Rs_period); }
         }
 
         if (this->restart)
@@ -216,10 +218,10 @@ OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
                 {
                     // Read HexxR in CSR format
                     if (GlobalC::exx_info.info_ri.real_number) {
-                        ModuleIO::read_Hexxs_csr(restart_HR_path, GlobalC::ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd);
+                        ModuleIO::read_Hexxs_csr(restart_HR_path, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd);
                     }
                     else {
-                        ModuleIO::read_Hexxs_csr(restart_HR_path, GlobalC::ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc);
+                        ModuleIO::read_Hexxs_csr(restart_HR_path, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc);
                     }
                 }
                 else
@@ -317,6 +319,7 @@ void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHk(int ik)
 
         if (GlobalC::exx_info.info_ri.real_number) {
             RI_2D_Comm::add_Hexx(
+                ucell,
                 this->kv,
                 ik,
                 GlobalC::exx_info.info_global.hybrid_alpha,
@@ -325,6 +328,7 @@ void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHk(int ik)
                 this->hsk->get_hk());
         } else {
             RI_2D_Comm::add_Hexx(
+                ucell,
                 this->kv,
                 ik,
                 GlobalC::exx_info.info_global.hybrid_alpha,

source/module_hamilt_lcao/hamilt_lcaodft/spar_exx.h

@@ -11,12 +11,12 @@
 
 #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_HS_arrays.hpp"
 #include "module_basis/module_ao/parallel_orbitals.h"
-
+#include "module_cell/unitcell.h"
 namespace sparse_format
 {
 
 template <typename Tdata>
-void cal_HR_exx(
+void cal_HR_exx(const UnitCell& ucell,
     const Parallel_Orbitals& pv,
     LCAO_HS_Arrays& HS_Arrays,
     const int& current_spin,

source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp

@@ -7,7 +7,8 @@
 #include "spar_exx.h"
 #include "spar_u.h"
 
-void sparse_format::cal_HSR(const Parallel_Orbitals& pv,
+void sparse_format::cal_HSR(const UnitCell& ucell,
+    const Parallel_Orbitals& pv,
     LCAO_HS_Arrays& HS_Arrays,
     Grid_Driver& grid,
     const int& current_spin,
@@ -99,7 +100,8 @@ void sparse_format::cal_HSR(const Parallel_Orbitals& pv,
     if (GlobalC::exx_info.info_global.cal_exx) {
         if (Hexxd && GlobalC::exx_info.info_ri.real_number)
         {
-            sparse_format::cal_HR_exx(pv,
+            sparse_format::cal_HR_exx(ucell,
+                pv,
                 HS_Arrays,
                 current_spin,
                 sparse_thr,
@@ -108,7 +110,8 @@ void sparse_format::cal_HSR(const Parallel_Orbitals& pv,
         }
         else if (Hexxc && !GlobalC::exx_info.info_ri.real_number)
         {
-            sparse_format::cal_HR_exx(pv,
+            sparse_format::cal_HR_exx(ucell,
+                pv,
                 HS_Arrays,
                 current_spin,
                 sparse_thr,

source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h

@@ -5,7 +5,8 @@
 
 namespace sparse_format {
     using TAC = std::pair<int, std::array<int, 3>>;
-    void cal_HSR(const Parallel_Orbitals& pv,
+    void cal_HSR(const UnitCell& ucell,
+        const Parallel_Orbitals& pv,
         LCAO_HS_Arrays& HS_Arrays,
         Grid_Driver& grid,
         const int& current_spin,

source/module_io/output_mat_sparse.cpp

@@ -45,7 +45,7 @@ void output_mat_sparse(const bool& out_mat_hsR,
     //! generate a file containing the Hamiltonian and S(overlap) matrices
     if (out_mat_hsR)
     {
-        output_HSR(istep, v_eff, pv, HS_Arrays, grid, kv, p_ham);
+        output_HSR(ucell,istep, v_eff, pv, HS_Arrays, grid, kv, p_ham);
     }
 
     //! generate a file containing the kinetic energy matrix

source/module_io/write_HS_R.cpp

@@ -12,7 +12,8 @@
 // The 'sparse_thr' is the accuracy of the sparse matrix.
 // If the absolute value of the matrix element is less than or equal to the
 // 'sparse_thr', it will be ignored.
-void ModuleIO::output_HSR(const int& istep,
+void ModuleIO::output_HSR(const UnitCell& ucell,
+    const int& istep,
     const ModuleBase::matrix& v_eff,
     const Parallel_Orbitals& pv,
     LCAO_HS_Arrays& HS_Arrays,
@@ -37,7 +38,7 @@ void ModuleIO::output_HSR(const int& istep,
     if (nspin == 1 || nspin == 4) {
         const int spin_now = 0;
         // jingan add 2021-6-4, modify 2021-12-2
-        sparse_format::cal_HSR(pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham
+        sparse_format::cal_HSR(ucell,pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham
 #ifdef __EXX
             , Hexxd, Hexxc
 #endif
@@ -47,7 +48,7 @@ void ModuleIO::output_HSR(const int& istep,
         int spin_now = 1;
 
         // save HR of spin down first (the current spin always be down)
-        sparse_format::cal_HSR(pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham
+        sparse_format::cal_HSR(ucell,pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham
 #ifdef __EXX
             , Hexxd, Hexxc
 #endif
@@ -61,7 +62,7 @@ void ModuleIO::output_HSR(const int& istep,
             spin_now = 0;
         }
 
-        sparse_format::cal_HSR(pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham
+        sparse_format::cal_HSR(ucell,pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham
 #ifdef __EXX
             , Hexxd, Hexxc
 #endif

source/module_io/write_HS_R.h

@@ -11,7 +11,8 @@
 namespace ModuleIO
 {
     using TAC = std::pair<int, std::array<int, 3>>;
-void output_HSR(const int& istep,
+void output_HSR(const UnitCell& ucell,
+    const int& istep,
     const ModuleBase::matrix& v_eff,
     const Parallel_Orbitals& pv,
     LCAO_HS_Arrays& HS_Arrays,

source/module_io/write_vxc.hpp

@@ -246,10 +246,10 @@ void write_Vxc(const int nspin,
     std::vector<std::vector<double>> e_orb_tot;   // orbital energy (total)
 #ifdef __EXX
     hamilt::OperatorEXX<hamilt::OperatorLCAO<TK, TR>> vexx_op_ao(&vxc_k_ao,
-        &vxcs_R_ao[0] /*for paraV*/, kv, Hexxd, Hexxc, hamilt::Add_Hexx_Type::k);
+        &vxcs_R_ao[0],ucell,/*for paraV*/ kv, Hexxd, Hexxc, hamilt::Add_Hexx_Type::k);
     hamilt::HS_Matrix_K<TK> vexxonly_k_ao(pv, 1); // only hk is needed, sk is skipped
     hamilt::OperatorEXX<hamilt::OperatorLCAO<TK, TR>> vexxonly_op_ao(&vexxonly_k_ao,
-        &vxcs_R_ao[0]/*for paraV*/, kv, Hexxd, Hexxc, hamilt::Add_Hexx_Type::k);
+        &vxcs_R_ao[0],ucell,/*for paraV*/ kv, Hexxd, Hexxc, hamilt::Add_Hexx_Type::k);
     std::vector<std::vector<double>> e_orb_exx; // orbital energy (EXX)
 #endif
     hamilt::OperatorDFTU<hamilt::OperatorLCAO<TK, TR>> vdftu_op_ao(&vxc_k_ao, kv.kvec_d, nullptr, kv.isk);