Refactor:Remove GlobalC::ucell in module_lr,module_psi (#5691)

A-006 · web-flow · commit f1d6055ef6b5 · 2024-12-06T22:13:15.000+08:00
* change ucell in module_lr

* change ucell in module_psi/wavefunc.cpp

* change module_psi/wf_atomic.cpp

* change module_pwdft/structure_factor_k.cpp

* fix bug in wavefunc
diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
@@ -219,7 +219,8 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
                                   this->kv.ngk.data(),
                                   this->pw_wfc->npwk_max,
                                   &this->sf,
-                                  &this->ppcell);
+                                  &this->ppcell,
+                                  ucell);
 
     this->kspw_psi = PARAM.inp.device == "gpu" || PARAM.inp.precision == "single"
                          ? new psi::Psi<T, Device>(this->psi[0])
@@ -257,7 +258,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
 
         this->pw_wfc->collect_local_pw(PARAM.inp.erf_ecut, PARAM.inp.erf_height, PARAM.inp.erf_sigma);
 
-        this->p_wf_init->make_table(this->kv.get_nks(), &this->sf, &this->ppcell);
+        this->p_wf_init->make_table(this->kv.get_nks(), &this->sf, &this->ppcell,ucell);
     }
     if (ucell.ionic_position_updated)
     {
@@ -373,6 +374,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
                                         this->kspw_psi,
                                         this->p_hamilt,
                                         this->ppcell,
+                                        ucell,
                                         GlobalV::ofs_running,
                                         this->already_initpsi);
 
diff --git a/source/module_hamilt_pw/hamilt_pwdft/structure_factor.cpp b/source/module_hamilt_pw/hamilt_pwdft/structure_factor.cpp
@@ -61,7 +61,7 @@ void Structure_Factor::setup_structure_factor(const UnitCell* Ucell, const Modul
     ModuleBase::TITLE("PW_Basis","setup_structure_factor");
     ModuleBase::timer::tick("PW_Basis","setup_struc_factor");
     const std::complex<double> ci_tpi = ModuleBase::NEG_IMAG_UNIT * ModuleBase::TWO_PI;
-
+    this->ucell = Ucell;
     this->strucFac.create(Ucell->ntype, rho_basis->npw);
     ModuleBase::Memory::record("SF::strucFac", sizeof(std::complex<double>) * Ucell->ntype*rho_basis->npw);
 
diff --git a/source/module_hamilt_pw/hamilt_pwdft/structure_factor.h b/source/module_hamilt_pw/hamilt_pwdft/structure_factor.h
@@ -53,6 +53,7 @@ class Structure_Factor
                                   ModuleBase::Vector3<double> q);
 
   private:
+    const UnitCell* ucell; 
     std::complex<float> * c_eigts1 = nullptr, * c_eigts2 = nullptr, * c_eigts3 = nullptr;
     std::complex<double> * z_eigts1 = nullptr, * z_eigts2 = nullptr, * z_eigts3 = nullptr;
     const ModulePW::PW_Basis* rho_basis = nullptr;
diff --git a/source/module_hamilt_pw/hamilt_pwdft/structure_factor_k.cpp b/source/module_hamilt_pw/hamilt_pwdft/structure_factor_k.cpp
@@ -10,7 +10,7 @@ std::complex<double>* Structure_Factor::get_sk(const int ik,
                                                const ModulePW::PW_Basis_K* wfc_basis) const
 {
     ModuleBase::timer::tick("Structure_Factor", "get_sk");
-    const double arg = (wfc_basis->kvec_c[ik] * GlobalC::ucell.atoms[it].tau[ia]) * ModuleBase::TWO_PI;
+    const double arg = (wfc_basis->kvec_c[ik] * ucell->atoms[it].tau[ia]) * ModuleBase::TWO_PI;
     const std::complex<double> kphase = std::complex<double>(cos(arg), -sin(arg));
     const int npw = wfc_basis->npwk[ik];
     std::complex<double> *sk = new std::complex<double>[npw];
@@ -26,19 +26,22 @@ std::complex<double>* Structure_Factor::get_sk(const int ik,
         const int ixy = wfc_basis->is2fftixy[is];
         int ix = ixy / wfc_basis->fftny;
         int iy = ixy % wfc_basis->fftny;
-        if (ix >= int(nx / 2) + 1) {
+        if (ix >= int(nx / 2) + 1) 
+        {
             ix -= nx;
-}
-        if (iy >= int(ny / 2) + 1) {
+        }
+        if (iy >= int(ny / 2) + 1) 
+        {
             iy -= ny;
-}
-        if (iz >= int(nz / 2) + 1) {
+        }
+        if (iz >= int(nz / 2) + 1) 
+        {
             iz -= nz;
-}
+        }
         ix += this->rho_basis->nx;
         iy += this->rho_basis->ny;
         iz += this->rho_basis->nz;
-        const int iat = GlobalC::ucell.itia2iat(it, ia);
+        const int iat = ucell->itia2iat(it, ia);
         sk[igl] = kphase * this->eigts1(iat, ix) * this->eigts2(iat, iy) * this->eigts3(iat, iz);
     }
     ModuleBase::timer::tick("Structure_Factor", "get_sk");
@@ -66,33 +69,33 @@ void Structure_Factor::get_sk(Device* ctx,
 
     int iat = 0, _npw = wfc_basis->npwk[ik], eigts1_nc = this->eigts1.nc, eigts2_nc = this->eigts2.nc,
             eigts3_nc = this->eigts3.nc;
-    int *igl2isz = nullptr, *is2fftixy = nullptr, *atom_na = nullptr, *h_atom_na = new int[GlobalC::ucell.ntype];
-    FPTYPE *atom_tau = nullptr, *h_atom_tau = new FPTYPE[GlobalC::ucell.nat * 3], *kvec = wfc_basis->get_kvec_c_data<FPTYPE>();
+    int *igl2isz = nullptr, *is2fftixy = nullptr, *atom_na = nullptr, *h_atom_na = new int[ucell->ntype];
+    FPTYPE *atom_tau = nullptr, *h_atom_tau = new FPTYPE[ucell->nat * 3], *kvec = wfc_basis->get_kvec_c_data<FPTYPE>();
     std::complex<FPTYPE> *eigts1 = this->get_eigts1_data<FPTYPE>(), *eigts2 = this->get_eigts2_data<FPTYPE>(),
             *eigts3 = this->get_eigts3_data<FPTYPE>();
-    for (int it = 0; it < GlobalC::ucell.ntype; it++)
+    for (int it = 0; it < ucell->ntype; it++)
     {
-        h_atom_na[it] = GlobalC::ucell.atoms[it].na;
+        h_atom_na[it] = ucell->atoms[it].na;
     }
 #ifdef _OPENMP
 #pragma omp parallel for
 #endif
-    for (int iat = 0; iat < GlobalC::ucell.nat; iat++)
+    for (int iat = 0; iat < ucell->nat; iat++)
     {
-        int it = GlobalC::ucell.iat2it[iat];
-        int ia = GlobalC::ucell.iat2ia[iat];
-        auto *tau = reinterpret_cast<double *>(GlobalC::ucell.atoms[it].tau.data());
+        int it = ucell->iat2it[iat];
+        int ia = ucell->iat2ia[iat];
+        auto *tau = reinterpret_cast<double *>(ucell->atoms[it].tau.data());
         h_atom_tau[iat * 3 + 0] = static_cast<FPTYPE>(tau[ia * 3 + 0]);
         h_atom_tau[iat * 3 + 1] = static_cast<FPTYPE>(tau[ia * 3 + 1]);
         h_atom_tau[iat * 3 + 2] = static_cast<FPTYPE>(tau[ia * 3 + 2]);
     }
     if (device == base_device::GpuDevice)
     {
-        resmem_int_op()(ctx, atom_na, GlobalC::ucell.ntype);
-        syncmem_int_op()(ctx, cpu_ctx, atom_na, h_atom_na, GlobalC::ucell.ntype);
+        resmem_int_op()(ctx, atom_na, ucell->ntype);
+        syncmem_int_op()(ctx, cpu_ctx, atom_na, h_atom_na, ucell->ntype);
 
-        resmem_var_op()(ctx, atom_tau, GlobalC::ucell.nat * 3);
-        syncmem_var_op()(ctx, cpu_ctx, atom_tau, h_atom_tau, GlobalC::ucell.nat * 3);
+        resmem_var_op()(ctx, atom_tau, ucell->nat * 3);
+        syncmem_var_op()(ctx, cpu_ctx, atom_tau, h_atom_tau, ucell->nat * 3);
 
         igl2isz = wfc_basis->d_igl2isz_k;
         is2fftixy = wfc_basis->d_is2fftixy;
@@ -107,7 +110,7 @@ void Structure_Factor::get_sk(Device* ctx,
 
     cal_sk_op()(ctx,
                 ik,
-                GlobalC::ucell.ntype,
+                ucell->ntype,
                 wfc_basis->nx,
                 wfc_basis->ny,
                 wfc_basis->nz,
@@ -152,7 +155,7 @@ std::complex<double>* Structure_Factor::get_skq(int ik,
     for (int ig = 0; ig < npw; ig++)
     {
         ModuleBase::Vector3<double> qkq = wfc_basis->getgpluskcar(ik, ig) + q;
-        double arg = (qkq * GlobalC::ucell.atoms[it].tau[ia]) * ModuleBase::TWO_PI;
+        double arg = (qkq * ucell->atoms[it].tau[ia]) * ModuleBase::TWO_PI;
         skq[ig] = std::complex<double>(cos(arg), -sin(arg));
     }
 
diff --git a/source/module_hsolver/test/hsolver_pw_sup.h b/source/module_hsolver/test/hsolver_pw_sup.h
@@ -192,6 +192,7 @@ void diago_PAO_in_pw_k2(
     wavefunc* p_wf,
     const ModuleBase::realArray& tab_at,
     const int& lmaxkb,
+    const UnitCell& ucell,
     hamilt::Hamilt<std::complex<float>, base_device::DEVICE_CPU>* phm_in) {
     for (int i = 0; i < wvf.size(); i++) {
         wvf.get_pointer()[i] = std::complex<float>((float)i + 1, 0);
@@ -207,6 +208,7 @@ void diago_PAO_in_pw_k2(
     wavefunc* p_wf,
     const ModuleBase::realArray& tab_at,
     const int& lmaxkb,
+    const UnitCell& ucell,
     hamilt::Hamilt<std::complex<double>, base_device::DEVICE_CPU>* phm_in) {
     for (int i = 0; i < wvf.size(); i++) {
         wvf.get_pointer()[i] = std::complex<double>((double)i + 1, 0);
diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp
@@ -255,7 +255,7 @@ LR::ESolver_LR<T, TR>::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu
     // necessary steps in ESolver_KS::before_all_runners : symmetry and k-points
     if (ModuleSymmetry::Symmetry::symm_flag == 1)
     {
-        GlobalC::ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
+        ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
     }
     this->kv.set(ucell,ucell.symm, PARAM.inp.kpoint_file, PARAM.inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
@@ -318,12 +318,12 @@ LR::ESolver_LR<T, TR>::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu
     this->init_pot(chg_gs);
 
     // search adjacent atoms and init Gint
-    std::cout << "ucell.infoNL.get_rcutmax_Beta(): " << GlobalC::ucell.infoNL.get_rcutmax_Beta() << std::endl;
+    std::cout << "ucell.infoNL.get_rcutmax_Beta(): " << ucell.infoNL.get_rcutmax_Beta() << std::endl;
     double search_radius = -1.0;
     search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running,
         PARAM.inp.out_level,
         orb.get_rcutmax_Phi(),
-        GlobalC::ucell.infoNL.get_rcutmax_Beta(),
+        ucell.infoNL.get_rcutmax_Beta(),
         PARAM.globalv.gamma_only_local);
     atom_arrange::search(PARAM.inp.search_pbc,
         GlobalV::ofs_running,
@@ -341,7 +341,7 @@ LR::ESolver_LR<T, TR>::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu
     std::vector<std::vector<double>> dpsi_u;
     std::vector<std::vector<double>> d2psi_u;
 
-    Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, orb, psi_u, dpsi_u, d2psi_u);
+    Gint_Tools::init_orb(dr_uniform, rcuts, ucell, orb, psi_u, dpsi_u, d2psi_u);
     this->gt_.set_pbc_grid(this->pw_rho->nx,
         this->pw_rho->ny,
         this->pw_rho->nz,
@@ -357,7 +357,7 @@ LR::ESolver_LR<T, TR>::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu
         this->pw_rho->ny,
         this->pw_rho->nplane,
         this->pw_rho->startz_current,
-        GlobalC::ucell,
+        ucell,
         GlobalC::GridD,
         dr_uniform,
         rcuts,
diff --git a/source/module_lr/operator_casida/operator_lr_exx.cpp b/source/module_lr/operator_casida/operator_lr_exx.cpp
@@ -65,7 +65,7 @@ namespace LR
         for (auto cell : this->BvK_cells)
         {
             std::complex<double> frac = RI::Global_Func::convert<std::complex<double>>(std::exp(
-                -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * (this->kv.kvec_c.at(ik) * (RI_Util::array3_to_Vector3(cell) * GlobalC::ucell.latvec))));
+                -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * (this->kv.kvec_c.at(ik) * (RI_Util::array3_to_Vector3(cell) * ucell.latvec))));
             for (int it1 = 0;it1 < ucell.ntype;++it1)
                 for (int ia1 = 0; ia1 < ucell.atoms[it1].na; ++ia1)
                     for (int it2 = 0;it2 < ucell.ntype;++it2)
diff --git a/source/module_lr/operator_casida/operator_lr_hxc.cpp b/source/module_lr/operator_casida/operator_lr_hxc.cpp
@@ -68,15 +68,15 @@ namespace LR
 
         // 3. v_hxc = f_hxc * rho_trans
         ModuleBase::matrix vr_hxc(1, nrxx);   //grid
-        this->pot.lock()->cal_v_eff(rho_trans, GlobalC::ucell, vr_hxc, ispin_ks);
+        this->pot.lock()->cal_v_eff(rho_trans, ucell, vr_hxc, ispin_ks);
         LR_Util::_deallocate_2order_nested_ptr(rho_trans, 1);
 
         // 4. V^{Hxc}_{\mu,\nu}=\int{dr} \phi_\mu(r) v_{Hxc}(r) \phi_\mu(r)
         Gint_inout inout_vlocal(vr_hxc.c, 0, Gint_Tools::job_type::vlocal);
         this->gint->get_hRGint()->set_zero();
         this->gint->cal_gint(&inout_vlocal);
         this->hR->set_zero();   // clear hR for each bands
-        this->gint->transfer_pvpR(&*this->hR, &GlobalC::ucell);    //grid to 2d block
+        this->gint->transfer_pvpR(&*this->hR, &ucell);    //grid to 2d block
         ModuleBase::timer::tick("OperatorLRHxc", "grid_calculation");
     }
 
@@ -88,7 +88,7 @@ namespace LR
 
         elecstate::DensityMatrix<std::complex<double>, double> DM_trans_real_imag(&pmat, 1, kv.kvec_d, kv.get_nks() / nspin);
         DM_trans_real_imag.init_DMR(*this->hR);
-        hamilt::HContainer<double> HR_real_imag(GlobalC::ucell, &this->pmat);
+        hamilt::HContainer<double> HR_real_imag(ucell, &this->pmat);
         LR_Util::initialize_HR<std::complex<double>, double>(HR_real_imag, ucell, gd, orb_cutoff_);
 
         auto dmR_to_hR = [&, this](const char& type) -> void
@@ -111,7 +111,7 @@ namespace LR
 
                 // 3. v_hxc = f_hxc * rho_trans
                 ModuleBase::matrix vr_hxc(1, nrxx);   //grid
-                this->pot.lock()->cal_v_eff(rho_trans, GlobalC::ucell, vr_hxc, ispin_ks);
+                this->pot.lock()->cal_v_eff(rho_trans, ucell, vr_hxc, ispin_ks);
                 // print_grid_nonzero(vr_hxc.c, this->poticab->nrxx, 10, "vr_hxc");
 
                 LR_Util::_deallocate_2order_nested_ptr(rho_trans, 1);
@@ -123,9 +123,9 @@ namespace LR
 
                 // LR_Util::print_HR(*this->gint->get_hRGint(), this->ucell.nat, "VR(grid)");
                 HR_real_imag.set_zero();
-                this->gint->transfer_pvpR(&HR_real_imag, &GlobalC::ucell, &GlobalC::GridD);
+                this->gint->transfer_pvpR(&HR_real_imag, &ucell, &GlobalC::GridD);
                 // LR_Util::print_HR(HR_real_imag, this->ucell.nat, "VR(real, 2d)");
-                LR_Util::set_HR_real_imag_part(HR_real_imag, *this->hR, GlobalC::ucell.nat, type);
+                LR_Util::set_HR_real_imag_part(HR_real_imag, *this->hR, ucell.nat, type);
             };
         this->hR->set_zero();
         dmR_to_hR('R');   //real
diff --git a/source/module_psi/psi_init.cpp b/source/module_psi/psi_init.cpp
@@ -93,7 +93,8 @@ void PSIInit<T, Device>::allocate_psi(Psi<std::complex<double>>*& psi,
                                       const int* ngk,
                                       const int npwx,
                                       Structure_Factor* p_sf,
-                                      pseudopot_cell_vnl* p_ppcell)
+                                      pseudopot_cell_vnl* p_ppcell,
+                                      const UnitCell& ucell)
 {
     // allocate memory for std::complex<double> datatype psi
     // New psi initializer in ABACUS, Developer's note:
@@ -126,7 +127,7 @@ void PSIInit<T, Device>::allocate_psi(Psi<std::complex<double>>*& psi,
         // however, init_at_1 does not actually initialize the psi, instead, it is a
         // function to calculate a interpolate table saving overlap intergral or say
         // Spherical Bessel Transform of atomic orbitals.
-        this->wf_old.init_at_1(p_sf, &p_ppcell->tab_at);
+        this->wf_old.init_at_1(ucell,p_sf, &p_ppcell->tab_at);
         // similarly, wfcinit not really initialize any wavefunction, instead, it initialize
         // the mapping from ixy, the 1d flattened index of point on fft grid (x, y) plane,
         // to the index of "stick", composed of grid points.
@@ -135,15 +136,18 @@ void PSIInit<T, Device>::allocate_psi(Psi<std::complex<double>>*& psi,
 }
 
 template <typename T, typename Device>
-void PSIInit<T, Device>::make_table(const int nks, Structure_Factor* p_sf, pseudopot_cell_vnl* p_ppcell)
+void PSIInit<T, Device>::make_table(const int nks, 
+                                    Structure_Factor* p_sf, 
+                                    pseudopot_cell_vnl* p_ppcell,
+                                    const UnitCell& ucell)
 {
     if (this->use_psiinitializer)
     {
     }    // do not need to do anything because the interpolate table is unchanged
     else // old initialization method, used in EXX calculation
     {
         this->wf_old.init_after_vc(nks); // reallocate wanf2, the planewave expansion of lcao
-        this->wf_old.init_at_1(p_sf, &p_ppcell->tab_at);    // re-calculate tab_at, the overlap matrix between atomic pswfc and jlq
+        this->wf_old.init_at_1(ucell,p_sf, &p_ppcell->tab_at);    // re-calculate tab_at, the overlap matrix between atomic pswfc and jlq
     }
 }
 
@@ -152,6 +156,7 @@ void PSIInit<T, Device>::initialize_psi(Psi<std::complex<double>>* psi,
                                         psi::Psi<T, Device>* kspw_psi,
                                         hamilt::Hamilt<T, Device>* p_hamilt,
                                         const pseudopot_cell_vnl& nlpp,
+                                        const UnitCell& ucell,
                                         std::ofstream& ofs_running,
                                         const bool is_already_initpsi)
 {
@@ -278,6 +283,7 @@ void PSIInit<T, Device>::initialize_psi(Psi<std::complex<double>>* psi,
                                                    &this->wf_old,
                                                    nlpp.tab_at,
                                                    nlpp.lmaxkb,
+                                                   ucell,
                                                    p_hamilt);
                     }
                 }
@@ -294,6 +300,7 @@ void PSIInit<T, Device>::initialize_psi(Psi<std::complex<double>>* psi,
                                                &this->wf_old,
                                                nlpp.tab_at,
                                                nlpp.lmaxkb,
+                                               ucell,
                                                p_hamilt);
                 }
             }
diff --git a/source/module_psi/psi_init.h b/source/module_psi/psi_init.h
@@ -36,10 +36,14 @@ class PSIInit
                       const int* ngk,                  //< number of G-vectors in the current pool
                       const int npwx,                  //< max number of plane waves of all pools
                       Structure_Factor* p_sf,          //< structure factor
-                      pseudopot_cell_vnl* p_ppcell);   //< nonlocal pseudopotential
+                      pseudopot_cell_vnl* p_ppcell,    //< nonlocal pseudopotential
+                      const UnitCell& ucell);          //< unit cell
 
     // make interpolate table
-    void make_table(const int nks, Structure_Factor* p_sf, pseudopot_cell_vnl* p_ppcell);
+    void make_table(const int nks, 
+                    Structure_Factor* p_sf, 
+                    pseudopot_cell_vnl* p_ppcell,
+                    const UnitCell& ucell);
 
     //------------------------ only for psi_initializer --------------------
     /**
@@ -54,6 +58,7 @@ class PSIInit
                         psi::Psi<T, Device>* kspw_psi,
                         hamilt::Hamilt<T, Device>* p_hamilt,
                         const pseudopot_cell_vnl& nlpp,
+                        const UnitCell& ucell,
                         std::ofstream& ofs_running,
                         const bool is_already_initpsi);
 
diff --git a/source/module_psi/wavefunc.cpp b/source/module_psi/wavefunc.cpp
diff --git a/source/module_psi/wavefunc.h b/source/module_psi/wavefunc.h
diff --git a/source/module_psi/wf_atomic.cpp b/source/module_psi/wf_atomic.cpp
diff --git a/source/module_psi/wf_atomic.h b/source/module_psi/wf_atomic.h

Original file line number	Diff line number	Diff line change
`@@ -65,7 +65,7 @@ namespace LR`
`65`	`65`	`for (auto cell : this->BvK_cells)`
`66`	`66`	`{`
`67`	`67`	`std::complex<double> frac = RI::Global_Func::convert<std::complex<double>>(std::exp(`
`68`		`- -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * (this->kv.kvec_c.at(ik) * (RI_Util::array3_to_Vector3(cell) * GlobalC::ucell.latvec))));`
	`68`	`+ -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * (this->kv.kvec_c.at(ik) * (RI_Util::array3_to_Vector3(cell) * ucell.latvec))));`
`69`	`69`	`for (int it1 = 0;it1 < ucell.ntype;++it1)`
`70`	`70`	`for (int ia1 = 0; ia1 < ucell.atoms[it1].na; ++ia1)`
`71`	`71`	`for (int it2 = 0;it2 < ucell.ntype;++it2)`
Original file line number	Diff line number	Diff line change
`@@ -93,7 +93,8 @@ void PSIInit<T, Device>::allocate_psi(Psi<std::complex<double>>*& psi,`
`93`	`93`	`const int* ngk,`
`94`	`94`	`const int npwx,`
`95`	`95`	`Structure_Factor* p_sf,`
`96`		`- pseudopot_cell_vnl* p_ppcell)`
	`96`	`+ pseudopot_cell_vnl* p_ppcell,`
	`97`	`+ const UnitCell& ucell)`
`97`	`98`	`{`
`98`	`99`	`// allocate memory for std::complex<double> datatype psi`
`99`	`100`	`// New psi initializer in ABACUS, Developer's note:`
`@@ -126,7 +127,7 @@ void PSIInit<T, Device>::allocate_psi(Psi<std::complex<double>>*& psi,`
`126`	`127`	`// however, init_at_1 does not actually initialize the psi, instead, it is a`
`127`	`128`	`// function to calculate a interpolate table saving overlap intergral or say`
`128`	`129`	`// Spherical Bessel Transform of atomic orbitals.`
`129`		`- this->wf_old.init_at_1(p_sf, &p_ppcell->tab_at);`
	`130`	`+ this->wf_old.init_at_1(ucell,p_sf, &p_ppcell->tab_at);`
`130`	`131`	`// similarly, wfcinit not really initialize any wavefunction, instead, it initialize`
`131`	`132`	`// the mapping from ixy, the 1d flattened index of point on fft grid (x, y) plane,`
`132`	`133`	`// to the index of "stick", composed of grid points.`
`@@ -135,15 +136,18 @@ void PSIInit<T, Device>::allocate_psi(Psi<std::complex<double>>*& psi,`
`135`	`136`	`}`
`136`	`137`
`137`	`138`	`template <typename T, typename Device>`
`138`		`-void PSIInit<T, Device>::make_table(const int nks, Structure_Factor* p_sf, pseudopot_cell_vnl* p_ppcell)`
	`139`	`+void PSIInit<T, Device>::make_table(const int nks,`
	`140`	`+ Structure_Factor* p_sf,`
	`141`	`+ pseudopot_cell_vnl* p_ppcell,`
	`142`	`+ const UnitCell& ucell)`
`139`	`143`	`{`
`140`	`144`	`if (this->use_psiinitializer)`
`141`	`145`	`{`
`142`	`146`	`} // do not need to do anything because the interpolate table is unchanged`
`143`	`147`	`else // old initialization method, used in EXX calculation`
`144`	`148`	`{`
`145`	`149`	`this->wf_old.init_after_vc(nks); // reallocate wanf2, the planewave expansion of lcao`
`146`		`- this->wf_old.init_at_1(p_sf, &p_ppcell->tab_at); // re-calculate tab_at, the overlap matrix between atomic pswfc and jlq`
	`150`	`+ this->wf_old.init_at_1(ucell,p_sf, &p_ppcell->tab_at); // re-calculate tab_at, the overlap matrix between atomic pswfc and jlq`
`147`	`151`	`}`
`148`	`152`	`}`
`149`	`153`
`@@ -152,6 +156,7 @@ void PSIInit<T, Device>::initialize_psi(Psi<std::complex<double>>* psi,`
`152`	`156`	`psi::Psi<T, Device>* kspw_psi,`
`153`	`157`	`hamilt::Hamilt<T, Device>* p_hamilt,`
`154`	`158`	`const pseudopot_cell_vnl& nlpp,`
	`159`	`+ const UnitCell& ucell,`
`155`	`160`	`std::ofstream& ofs_running,`
`156`	`161`	`const bool is_already_initpsi)`
`157`	`162`	`{`
`@@ -278,6 +283,7 @@ void PSIInit<T, Device>::initialize_psi(Psi<std::complex<double>>* psi,`
`278`	`283`	`&this->wf_old,`
`279`	`284`	`nlpp.tab_at,`
`280`	`285`	`nlpp.lmaxkb,`
	`286`	`+ ucell,`
`281`	`287`	`p_hamilt);`
`282`	`288`	`}`
`283`	`289`	`}`
`@@ -294,6 +300,7 @@ void PSIInit<T, Device>::initialize_psi(Psi<std::complex<double>>* psi,`
`294`	`300`	`&this->wf_old,`
`295`	`301`	`nlpp.tab_at,`
`296`	`302`	`nlpp.lmaxkb,`
	`303`	`+ ucell,`
`297`	`304`	`p_hamilt);`
`298`	`305`	`}`
`299`	`306`	`}`