diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt
index 7273c6b8d72..bbc6a9025f1 100644
--- a/source/CMakeLists.txt
+++ b/source/CMakeLists.txt
@@ -615,7 +615,7 @@ target_link_libraries(
     cell
     parameter
     psi_overall_init
-    psi_initializer
+    psi_init
     psi
     dftu
     deltaspin
diff --git a/source/Makefile.Objects b/source/Makefile.Objects
index d6a33659185..6fb86642274 100644
--- a/source/Makefile.Objects
+++ b/source/Makefile.Objects
@@ -443,7 +443,7 @@ OBJS_ORBITAL=ORB_atomic.o\
 
 OBJS_PSI=psi.o\
 
-OBJS_PSI_INITIALIZER=psi_initializer.o\
+OBJS_PSI_INITIALIZER=psi_base.o\
                      psi_init_random.o\
                      psi_init_file.o\
                      psi_init_atomic.o\
diff --git a/source/source_cell/test/klist_test.cpp b/source/source_cell/test/klist_test.cpp
index 79bea7c2433..c05744bb785 100644
--- a/source/source_cell/test/klist_test.cpp
+++ b/source/source_cell/test/klist_test.cpp
@@ -70,7 +70,7 @@ pseudopot_cell_vnl::pseudopot_cell_vnl()
 pseudopot_cell_vnl::~pseudopot_cell_vnl()
 {
 }
-Soc::~Soc()
+SocBase::~SocBase()
 {
 }
 Fcoef::~Fcoef()
diff --git a/source/source_cell/test/klist_test_para.cpp b/source/source_cell/test/klist_test_para.cpp
index d5417cc5492..3415cadc7e7 100644
--- a/source/source_cell/test/klist_test_para.cpp
+++ b/source/source_cell/test/klist_test_para.cpp
@@ -73,7 +73,7 @@ pseudopot_cell_vnl::pseudopot_cell_vnl()
 pseudopot_cell_vnl::~pseudopot_cell_vnl()
 {
 }
-Soc::~Soc()
+SocBase::~SocBase()
 {
 }
 Fcoef::~Fcoef()
diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp
index f08d2aa99f3..743ac4c9c62 100644
--- a/source/source_esolver/esolver_ks_pw.cpp
+++ b/source/source_esolver/esolver_ks_pw.cpp
@@ -93,7 +93,7 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
       this->locpp, this->ppcell, this->vsep_cell, this->pw_wfc, this->pw_rho,
       this->pw_rhod, this->pw_big, this->solvent, inp);
 
-    this->stp.before_runner(ucell, this->kv, this->sf, *this->pw_wfc, this->ppcell, PARAM.inp);
+    this->stp.before_runner(ucell, this->kv, this->sf, *this->pw_wfc, this->ppcell.lmaxkb, PARAM.inp);
 
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS");
 
@@ -143,7 +143,18 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
 
     if (ucell.cell_parameter_updated)
     {
-        this->stp.p_psi_init->prepare_init(PARAM.inp.pw_seed);
+        this->stp.p_psi_init->prepare_init(PARAM.inp.pw_seed,
+                                            PARAM.globalv.nbands_l,
+                                            PARAM.inp.nspin,
+                                            PARAM.globalv.nqx,
+                                            PARAM.globalv.dq,
+                                            PARAM.globalv.domag,
+                                            PARAM.globalv.domag_z,
+                                            PARAM.inp.pseudo_mesh,
+                                            PARAM.inp.orbital_dir,
+                                            PARAM.globalv.global_readin_dir,
+                                            PARAM.globalv.ks_run,
+                                            PARAM.globalv.npol);
     }
 
     //! Init Hamiltonian (cell changed)
@@ -168,7 +179,16 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
 	      this->pw_wfc, this->pw_rhod, PARAM.globalv.global_out_dir, PARAM.inp);
 
     // setup psi (electronic wave functions)
-    this->stp.init(this->p_hamilt);
+    this->stp.init(this->p_hamilt,
+                     PARAM.inp.device,
+                     PARAM.inp.precision,
+                     PARAM.inp.ks_solver,
+                     PARAM.inp.bndpar,
+                     PARAM.inp.use_k_continuity,
+                     PARAM.inp.calculation,
+                     PARAM.inp.mem_saver,
+                     PARAM.globalv.npol,
+                     PARAM.globalv.ks_run);
 
     //! Setup EXX helper for Hamiltonian and psi
     exx_helper->before_scf(this->p_hamilt, this->stp.template get_psi_t<T, Device>(), PARAM.inp);
diff --git a/source/source_estate/test/elecstate_pw_test.cpp b/source/source_estate/test/elecstate_pw_test.cpp
index b1e9e354c30..837ffe37ce9 100644
--- a/source/source_estate/test/elecstate_pw_test.cpp
+++ b/source/source_estate/test/elecstate_pw_test.cpp
@@ -8,7 +8,7 @@
 #include "source_hamilt/module_xc/xc_functional.h"
 #include "source_pw/module_pwdft/vl_pw.h"
 #include "source_pw/module_pwdft/vnl_pw.h"
-#include "source_pw/module_pwdft/soc.h"
+#include "source_psi/soc_base.h"
 #include "source_io/module_parameter/parameter.h"
 // mock functions for testing
 int XC_Functional::func_type = 1;
@@ -115,7 +115,7 @@ void pseudopot_cell_vnl::getvnl<double, base_device::DEVICE_CPU>(base_device::DE
                                                                  std::complex<double>*) const
 {
 }
-Soc::~Soc()
+SocBase::~SocBase()
 {
 }
 Fcoef::~Fcoef()
diff --git a/source/source_io/module_ctrl/ctrl_output_pw.h b/source/source_io/module_ctrl/ctrl_output_pw.h
index 00b7509990c..ffed81dfe33 100644
--- a/source/source_io/module_ctrl/ctrl_output_pw.h
+++ b/source/source_io/module_ctrl/ctrl_output_pw.h
@@ -4,7 +4,9 @@
 #include "source_base/module_device/device.h" // use Device
 #include "source_psi/psi.h"                   // define psi
 #include "source_estate/elecstate_lcao.h"     // use pelec
-#include "source_psi/setup_psi_pw.h" // use Setup_Psi class
+#include "source_psi/setup_psi_pw.h"          // use Setup_Psi class
+
+class pseudopot_cell_vnl;
 
 namespace ModuleIO
 {
diff --git a/source/source_io/module_parameter/read_input_item_postprocess.cpp b/source/source_io/module_parameter/read_input_item_postprocess.cpp
index aad6f47a85a..6955d840410 100644
--- a/source/source_io/module_parameter/read_input_item_postprocess.cpp
+++ b/source/source_io/module_parameter/read_input_item_postprocess.cpp
@@ -403,7 +403,7 @@ void ReadInput::item_postprocess()
 
                        In the future lcao_in_pw will have its own ESolver.
 
-                       2023/12/22 use new psi_initializer to expand numerical
+                       2023/12/22 use new psi_base to expand numerical
                       atomic orbitals, ykhuang
                    */
             if (para.input.towannier90 && para.input.basis_type == "lcao_in_pw")
diff --git a/source/source_io/module_wannier/to_wannier90_lcao_in_pw.cpp b/source/source_io/module_wannier/to_wannier90_lcao_in_pw.cpp
index 48d3e7ae65c..3826308f235 100644
--- a/source/source_io/module_wannier/to_wannier90_lcao_in_pw.cpp
+++ b/source/source_io/module_wannier/to_wannier90_lcao_in_pw.cpp
@@ -41,11 +41,12 @@ void toWannier90_LCAO_IN_PW::calculate(
 {
     this->ParaV = pv;
 
-    Structure_Factor* sf_ptr = const_cast<Structure_Factor*>(&sf);
     ModulePW::PW_Basis_K* wfcpw_ptr = const_cast<ModulePW::PW_Basis_K*>(wfcpw);
     delete this->psi_initer_;
-    this->psi_initer_ = new psi_init_nao<std::complex<double>>();
-    this->psi_initer_->initialize(sf_ptr, wfcpw_ptr, &ucell, &kv, 1, nullptr, GlobalV::MY_RANK);
+    psi_init_nao<std::complex<double>>* nao_initer = new psi_init_nao<std::complex<double>>();
+    nao_initer->prepare_params(PARAM.globalv.nqx, PARAM.globalv.dq, PARAM.inp.nspin, PARAM.inp.orbital_dir, sf);
+    this->psi_initer_ = nao_initer;
+    this->psi_initer_->initialize(wfcpw_ptr, &ucell, kv.ik2iktot, kv.get_nkstot(), 1, 0, GlobalV::MY_RANK, PARAM.globalv.npol, PARAM.inp.nbands);
     this->psi_initer_->tabulate();
     delete this->psi;
     const int nks_psi = (PARAM.inp.calculation == "nscf" && PARAM.inp.mem_saver == 1)? 1 : wfcpw->nks;
diff --git a/source/source_io/module_wannier/to_wannier90_lcao_in_pw.h b/source/source_io/module_wannier/to_wannier90_lcao_in_pw.h
index a6f204858f1..bb0100c1b1b 100644
--- a/source/source_io/module_wannier/to_wannier90_lcao_in_pw.h
+++ b/source/source_io/module_wannier/to_wannier90_lcao_in_pw.h
@@ -21,7 +21,7 @@
 
 #ifdef __LCAO
 #include "source_basis/module_ao/parallel_orbitals.h"
-#include "source_psi/psi_initializer.h"
+#include "source_psi/psi_base.h"
 
 class toWannier90_LCAO_IN_PW : public toWannier90_PW
 {
@@ -59,7 +59,7 @@ class toWannier90_LCAO_IN_PW : public toWannier90_PW
   protected:
     const Parallel_Orbitals* ParaV = nullptr;
     /// @brief psi initializer for expanding nao in planewave basis
-    psi_initializer<std::complex<double>>* psi_initer_ = nullptr;
+    psi_base<std::complex<double>>* psi_initer_ = nullptr;
 
     psi::Psi<std::complex<double>, base_device::DEVICE_CPU>* psi = nullptr;
 
diff --git a/source/source_io/test/for_testing_input_conv.h b/source/source_io/test/for_testing_input_conv.h
index 893f94f76c1..917a36a2f93 100644
--- a/source/source_io/test/for_testing_input_conv.h
+++ b/source/source_io/test/for_testing_input_conv.h
@@ -107,7 +107,7 @@ pseudopot_cell_vnl::pseudopot_cell_vnl()
 pseudopot_cell_vnl::~pseudopot_cell_vnl()
 {
 }
-Soc::~Soc()
+SocBase::~SocBase()
 {
 }
 Fcoef::~Fcoef()
diff --git a/source/source_io/test/for_testing_klist.h b/source/source_io/test/for_testing_klist.h
index 0a0e2a14c39..48850444a72 100644
--- a/source/source_io/test/for_testing_klist.h
+++ b/source/source_io/test/for_testing_klist.h
@@ -30,7 +30,7 @@ pseudopot_cell_vl::pseudopot_cell_vl(){}
 pseudopot_cell_vl::~pseudopot_cell_vl(){}
 pseudopot_cell_vnl::pseudopot_cell_vnl(){}
 pseudopot_cell_vnl::~pseudopot_cell_vnl(){}
-Soc::~Soc()
+SocBase::~SocBase()
 {
 }
 Fcoef::~Fcoef()
diff --git a/source/source_lcao/module_deepks/test/CMakeLists.txt b/source/source_lcao/module_deepks/test/CMakeLists.txt
index 99389de9bdf..3e70dd80fb1 100644
--- a/source/source_lcao/module_deepks/test/CMakeLists.txt
+++ b/source/source_lcao/module_deepks/test/CMakeLists.txt
@@ -37,7 +37,6 @@ set(DEEPKS_UNIT_COMMON_SOURCES
     ../../../source_cell/read_pp_blps.cpp
     ../../../source_cell/sep.cpp
     ../../../source_cell/sep_cell.cpp
-    ../../../source_pw/module_pwdft/soc.cpp
     ../../../source_io/module_output/sparse_matrix.cpp
     ../../../source_cell/read_pseudo.cpp
     ../../../source_estate/param_update.cpp
diff --git a/source/source_lcao/module_rdmft/CMakeLists.txt b/source/source_lcao/module_rdmft/CMakeLists.txt
index c936787ba62..9223c2d6169 100644
--- a/source/source_lcao/module_rdmft/CMakeLists.txt
+++ b/source/source_lcao/module_rdmft/CMakeLists.txt
@@ -11,7 +11,7 @@ endif()
 
 # if(ENABLE_COVERAGE)
 #   add_coverage(psi)
-#   add_coverage(psi_initializer)
+#   add_coverage(psi_init)
 # endif()
 
 # if (BUILD_TESTING)
diff --git a/source/source_lcao/module_ri/exx_lip.hpp b/source/source_lcao/module_ri/exx_lip.hpp
index 5a86681c21b..120aec65164 100644
--- a/source/source_lcao/module_ri/exx_lip.hpp
+++ b/source/source_lcao/module_ri/exx_lip.hpp
@@ -19,7 +19,7 @@
 #include "source_estate/elecstate.h"
 #include "source_basis/module_pw/pw_basis_k.h"
 #include "source_cell/module_symmetry/symmetry.h"
-#include "source_psi/psi_initializer.h"
+#include "source_psi/psi_base.h"
 #include "source_pw/module_pwdft/structure_factor.h"
 #include "source_base/tool_title.h"
 #include "source_base/timer.h"
diff --git a/source/source_lcao/module_rt/test/CMakeLists.txt b/source/source_lcao/module_rt/test/CMakeLists.txt
index 1cc5455b003..985793690a6 100644
--- a/source/source_lcao/module_rt/test/CMakeLists.txt
+++ b/source/source_lcao/module_rt/test/CMakeLists.txt
@@ -52,5 +52,5 @@ AddTest(
           ../../../source_io/module_hs/cal_r_overlap_R.cpp
           ../../../source_io/module_hs/single_R_io.cpp
           ../../../source_io/module_hs/rr_sparse_writer.cpp
-          ../../../source_pw/module_pwdft/soc.cpp
+          ../../../source_psi/soc_base.cpp
 )
diff --git a/source/source_lcao/setup_nonlocal.cpp b/source/source_lcao/setup_nonlocal.cpp
index 7a93b65308e..0b1b9f91b3f 100644
--- a/source/source_lcao/setup_nonlocal.cpp
+++ b/source/source_lcao/setup_nonlocal.cpp
@@ -4,7 +4,7 @@
 #include "source_io/module_parameter/parameter.h"
 
 #ifdef __LCAO
-#include "source_pw/module_pwdft/soc.h"
+#include "source_psi/soc_base.h"
 // mohan add 2013-08-02
 // In order to get rid of the read in file .NONLOCAL.
 
@@ -53,7 +53,7 @@ void InfoNonlocal::Set_NonLocal(const int& it,
     {
         lmaxkb = std::max(lmaxkb, atom->ncpp.lll[ibeta]);
     }
-    Soc soc;
+    SocBase soc;
     if (atom->ncpp.has_so)
     {
         soc.rot_ylm(lmaxkb);
diff --git a/source/source_psi/CMakeLists.txt b/source/source_psi/CMakeLists.txt
index 8be3a2ba396..8bdcf5354c9 100644
--- a/source/source_psi/CMakeLists.txt
+++ b/source/source_psi/CMakeLists.txt
@@ -2,6 +2,7 @@ add_library(
     psi
     OBJECT
     psi.cpp
+    soc_base.cpp
 )
 
 add_library(
@@ -13,9 +14,9 @@ add_library(
 )
 
 add_library(
-  psi_initializer
+  psi_init
   OBJECT
-  psi_initializer.cpp
+  psi_base.cpp
   psi_init_random.cpp
   psi_init_file.cpp
   psi_init_atomic.cpp
@@ -27,7 +28,7 @@ add_library(
 
 if(ENABLE_COVERAGE)
   add_coverage(psi)
-  add_coverage(psi_initializer)
+  add_coverage(psi_init)
 endif()
 
 if (BUILD_TESTING)
diff --git a/source/source_psi/psi_initializer.cpp b/source/source_psi/psi_base.cpp
similarity index 86%
rename from source/source_psi/psi_initializer.cpp
rename to source/source_psi/psi_base.cpp
index bed67ccd1c3..f5bc4d44a36 100644
--- a/source/source_psi/psi_initializer.cpp
+++ b/source/source_psi/psi_base.cpp
@@ -1,5 +1,11 @@
-#include "psi_initializer.h"
+#include "psi_base.h"
 
+#include <vector>
+#include <cassert>
+#include <cstdlib>
+
+#include "source_cell/unitcell.h"
+#include "source_basis/module_pw/pw_basis_k.h"
 #include "source_base/parallel_global.h"
 // basic functions support
 #include "source_base/timer.h"
@@ -13,38 +19,42 @@
 #endif
 
 template <typename T>
-void psi_initializer<T>::initialize(const Structure_Factor* sf,
-                                    const ModulePW::PW_Basis_K* pw_wfc,
+void psi_base<T>::initialize(const ModulePW::PW_Basis_K* pw_wfc,
                                     const UnitCell* p_ucell,
-                                    const K_Vectors* p_kv_in,
+                                    const std::vector<int>& ik2iktot,
+                                    const int& nkstot,
                                     const int& random_seed,
-                                    const pseudopot_cell_vnl* p_pspot_nl,
-                                    const int& rank)
+                                    const int& lmaxkb,
+                                    const int& rank,
+                                    const int& npol,
+                                    const int& nbands)
 {
-    this->sf_ = sf;
     this->pw_wfc_ = pw_wfc;
     this->p_ucell_ = p_ucell;
-    this->p_kv = p_kv_in;
+    this->ik2iktot_ = ik2iktot;
+    this->nkstot_ = nkstot;
     this->random_seed_ = random_seed;
-    this->p_pspot_nl_ = p_pspot_nl;
+    this->lmaxkb_ = lmaxkb;
+    this->npol_ = npol;
+    this->nbands_ = nbands;
 }
 
 template <typename T>
-void psi_initializer<T>::random_t(T* psi, const int iw_start, const int iw_end, const int ik, const int mode)
+void psi_base<T>::random_t(T* psi, const int iw_start, const int iw_end, const int ik, const int mode)
 {
     ModuleBase::timer::start("psi_init", "random_t");
     assert(mode <= 1);
     assert(iw_start >= 0);
     const int ng = this->pw_wfc_->npwk[ik];
     const int npwk_max = this->pw_wfc_->npwk_max;
-    const int npol = PARAM.globalv.npol;
+    const int npol = this->npol_;
 
     // If random seed is specified, then generate random wavefunction satisfying that
     // it can generate the same results using different number of processors.
     if (this->random_seed_ > 0) // qianrui add 2021-8-13
     {
 #ifdef __MPI
-        srand(unsigned(this->random_seed_ + this->p_kv->ik2iktot[ik]));
+        srand(unsigned(this->random_seed_ + this->ik2iktot_[ik]));
 #else
         srand(unsigned(this->random_seed_ + ik));
 #endif
@@ -184,7 +194,7 @@ void psi_initializer<T>::random_t(T* psi, const int iw_start, const int iw_end,
 
 #ifdef __MPI
 template <typename T>
-void psi_initializer<T>::stick_to_pool(Real* stick, const int& ir, Real* out) const
+void psi_base<T>::stick_to_pool(Real* stick, const int& ir, Real* out) const
 {
     ModuleBase::timer::start("psi_init", "stick_to_pool");
     MPI_Status ierror;
@@ -211,7 +221,7 @@ void psi_initializer<T>::stick_to_pool(Real* stick, const int& ir, Real* out) co
         }
         else
         {
-            ModuleBase::WARNING_QUIT("psi_initializer", "stick_to_pool: Real type not supported");
+            ModuleBase::WARNING_QUIT("psi_base", "stick_to_pool: Real type not supported");
         }
         for (int iz = 0; iz < nz; iz++)
         {
@@ -230,7 +240,7 @@ void psi_initializer<T>::stick_to_pool(Real* stick, const int& ir, Real* out) co
         }
         else
         {
-            ModuleBase::WARNING_QUIT("psi_initializer", "stick_to_pool: Real type not supported");
+            ModuleBase::WARNING_QUIT("psi_base", "stick_to_pool: Real type not supported");
         }
     }
 
@@ -240,8 +250,8 @@ void psi_initializer<T>::stick_to_pool(Real* stick, const int& ir, Real* out) co
 #endif
 
 // explicit instantiation
-template class psi_initializer<std::complex<double>>;
-template class psi_initializer<std::complex<float>>;
+template class psi_base<std::complex<double>>;
+template class psi_base<std::complex<float>>;
 // gamma point calculation
-template class psi_initializer<double>;
-template class psi_initializer<float>;
+template class psi_base<double>;
+template class psi_base<float>;
diff --git a/source/source_psi/psi_initializer.h b/source/source_psi/psi_base.h
similarity index 68%
rename from source/source_psi/psi_initializer.h
rename to source/source_psi/psi_base.h
index 589355d2710..e77d529c08f 100644
--- a/source/source_psi/psi_initializer.h
+++ b/source/source_psi/psi_base.h
@@ -1,22 +1,21 @@
-#ifndef PSI_INITIALIZER_H
-#define PSI_INITIALIZER_H
-// data structure support
-#include "source_basis/module_pw/pw_basis_k.h" // for kpoint related data structure
-#include "source_pw/module_pwdft/vnl_pw.h"
-#include "source_pw/module_pwdft/structure_factor.h"
-#include "source_psi/psi.h" // for psi data structure
-// smart pointer for auto-memory management
+#ifndef PSI_BASE_H
+#define PSI_BASE_H
+#include "source_basis/module_pw/pw_basis_k.h"
+#include "source_psi/psi.h" 
 #include <memory>
-// numerical algorithm support
 #ifdef __MPI
 #include <mpi.h>
 #endif
 #include "source_base/macros.h"
-#include "source_cell/klist.h"
+#include "source_cell/unitcell.h"
 
 #include <type_traits>
+#include <vector>
+
+using namespace std;
+
 /*
-Psi (planewave based wavefunction) initializer
+Psi (planewave based wavefunction) base class
 Auther: Kirk0830
 Institute: AI for Science Institute, BEIJING
 
@@ -24,14 +23,14 @@ This class is used to allocate memory and give initial guess for psi
 therefore only double datatype is needed to be supported.
 Following methods are available:
     1. file: use wavefunction file to initialize psi
-             implemented in psi_initializer_file.h
+             implemented in psi_init_file.h
     2. random: use random number to initialize psi
-               implemented in psi_initializer_random.h
+               implemented in psi_init_random.h
     3. atomic: use pseudo-wavefunction in pseudopotential file to initialize psi
-               implemented in psi_initializer_atomic.h
+               implemented in psi_init_atomic.h
     4. atomic+random: mix 'atomic' with some random numbers to initialize psi
     5. nao: use numerical orbitals to initialize psi
-            implemented in psi_initializer_nao.h
+            implemented in psi_init_nao.h
     6. nao+random: mix 'nao' with some random numbers to initialize psi
 
 To use:
@@ -39,30 +38,32 @@ To use:
 A practical example would be in ESolver_KS_PW, because polymorphism is achieved by
 pointer, while a raw pointer is risky, therefore std::unique_ptr is a better
 choice.
-1. new a std::unique_ptr<psi_initializer<T> with specific derived class
-2. initialize() to link psi_initializer with external data and methods
+1. new a std::unique_ptr<psi_base<T> with specific derived class
+2. initialize() to link psi_base with external data and methods
 3. tabulate() to calculate the interpolate table
 4. init_psig() to calculate projection of atomic radial function onto planewave basis
 In summary:
 new->initialize->tabulate->init_psig
 */
 template <typename T>
-class psi_initializer
+class psi_base
 {
   private:
     using Real = typename GetTypeReal<T>::type;
 
   public:
-    psi_initializer(){};
-    virtual ~psi_initializer(){};
-    /// @brief initialize the psi_initializer with external data and methods
-    virtual void initialize(const Structure_Factor*,             //< structure factor
-                            const ModulePW::PW_Basis_K*,         //< planewave basis
-                            const UnitCell*,                     //< unit cell
-                            const K_Vectors* = nullptr,          //< parallel kpoints
-                            const int& = 1,                      //< random seed
-                            const pseudopot_cell_vnl* = nullptr, //< nonlocal pseudopotential
-                            const int& = 0);                     //< rank
+    psi_base(){};
+    virtual ~psi_base(){};
+    /// @brief initialize the psi_base with external data and methods
+    virtual void initialize(const ModulePW::PW_Basis_K* pw_wfc,         //< planewave basis
+                            const UnitCell* p_ucell,                     //< unit cell
+                            const std::vector<int>& ik2iktot,             //< ik2iktot: local->global k-point mapping
+                            const int& nkstot,                          //< nkstot: total number of k-points
+                            const int& random_seed,                      //< random seed
+                            const int& lmaxkb,                          //< lmaxkb: max angular momentum for non-local projectors
+                            const int& rank,                            //< rank
+                            const int& npol,                            //< npol
+                            const int& nbands);                         //< nbands
 
     /// @brief CENTRAL FUNCTION: calculate the interpolate table if needed
     virtual void tabulate()
@@ -112,6 +113,7 @@ class psi_initializer
     }
 
   protected:
+
 #ifdef __MPI // MPI additional implementation
     /// @brief mapping from (ix, iy) to is
     void stick_to_pool(Real* stick,      //< stick
@@ -123,17 +125,33 @@ class psi_initializer
                   const int iw_end,                  ///< iw_end, ending band index
                   const int ik,                      ///< ik, kpoint index
                   const int mode = 1);               ///< mode, 0 for rr*exp(i*arg), 1 for rr/(1+gk2)*exp(i*arg)
-    const Structure_Factor* sf_ = nullptr;           ///< Structure_Factor
+
     const ModulePW::PW_Basis_K* pw_wfc_ = nullptr;   ///< use |k+G>, |G>, getgpluskcar and so on in PW_Basis_K
+
     const UnitCell* p_ucell_ = nullptr;              ///< UnitCell
-    const K_Vectors* p_kv = nullptr;                 ///< Parallel_Kpoints
-    const pseudopot_cell_vnl* p_pspot_nl_ = nullptr; ///< pseudopot_cell_vnl
+
+    int lmaxkb_ = 0;                                 ///< max angular momentum for non-local projectors
+
+    std::vector<int> ik2iktot_;                      ///< local->global k-point mapping
+
+    int nkstot_ = 0;                                 ///< total number of k-points
+
     int random_seed_ = 1;                            ///< random seed, shared by random, atomic+random, nao+random
+
     std::vector<int> ixy2is_;                        ///< used by stick_to_pool function
+
     int mem_saver_ = 0;                              ///< if save memory, only for nscf
+
     std::string method_ = "none";                    ///< method name
+
     int nbands_complem_ = 0; ///< complement number of bands, which is nbands_start_ - ucell.natomwfc
+
     double mixing_coef_ = 0; ///< mixing coefficient for atomic+random and nao+random
+
     int nbands_start_ = 0;   ///< starting nbands, which is no less than PARAM.inp.nbands
+
+    int npol_ = 1;           ///< number of polarizations
+
+    int nbands_ = 1;         ///< number of bands
 };
-#endif
\ No newline at end of file
+#endif
diff --git a/source/source_psi/psi_init_atomic.cpp b/source/source_psi/psi_init_atomic.cpp
index fa110f970a6..446f8a9312a 100644
--- a/source/source_psi/psi_init_atomic.cpp
+++ b/source/source_psi/psi_init_atomic.cpp
@@ -1,35 +1,18 @@
 #include "psi_init_atomic.h"
-#include "source_pw/module_pwdft/soc.h"
-// numerical algorithm support
+#include "soc_base.h"
 #include "source_base/math_integral.h" // for numerical integration
 #include "source_base/math_polyint.h" // for polynomial interpolation
 #include "source_base/math_ylmreal.h" // for real spherical harmonics
 #include "source_base/math_sphbes.h" // for spherical bessel functions
-// basic functions support
 #include "source_base/tool_quit.h"
 #include "source_base/timer.h"
-// global variables definition
 #include "source_base/global_variable.h"
-#include "source_io/module_parameter/parameter.h"
-// io support
 #include "source_io/module_output/write_pao.h"
 
-// free function, compared with common radial function normalization, it does not multiply r to function
-// due to pswfc is already multiplied by r
-// template <typename T>
-// void normalize(int n_rgrid, std::vector<T>& pswfcr, double* rab)
-// {
-//     std::vector<T> pswfc2r2(pswfcr.size());
-//     std::transform(pswfcr.begin(), pswfcr.end(), pswfc2r2.begin(), [](T pswfc) { return pswfc * pswfc; });
-//     T norm = ModuleBase::Integral::simpson(n_rgrid, pswfc2r2.data(), rab);
-//     norm = sqrt(norm);
-//     std::transform(pswfcr.begin(), pswfcr.end(), pswfcr.begin(), [norm](T pswfc) { return pswfc / norm; });
-// }
-
 template <typename T>
 void psi_init_atomic<T>::allocate_ps_table()
 {
-   // find correct dimension for ovlp_flzjlq
+    // find correct dimension for ovlp_flzjlq
     int dim1 = this->p_ucell_->ntype;
     int dim2 = 0; // dim2 should be the maximum number of pseudo atomic orbitals
     for (int it = 0; it < this->p_ucell_->ntype; it++)
@@ -38,36 +21,66 @@ void psi_init_atomic<T>::allocate_ps_table()
     }
     if (dim2 == 0)
     {
-        ModuleBase::WARNING_QUIT("psi_init_atomic<T>::allocate_table", "there is not ANY pseudo atomic orbital read in present system, recommand other methods, quit.");
+        ModuleBase::WARNING_QUIT("psi_init_atomic<T>::allocate_table", 
+            "there is not ANY pseudo atomic orbital read in present system, recommand other methods, quit.");
+    }
+    if (this->nqx_ <= 0)
+    {
+        ModuleBase::WARNING_QUIT("psi_init_atomic<T>::allocate_ps_table",
+            "nqx_ must be greater than 0. Did you forget to call prepare_params() before initialize()?");
     }
-    int dim3 = PARAM.globalv.nqx;
+    int dim3 = this->nqx_;
     // allocate memory for ovlp_flzjlq
     this->ovlp_pswfcjlq_.create(dim1, dim2, dim3);
     this->ovlp_pswfcjlq_.zero_out();
 }
 
 template <typename T>
-void psi_init_atomic<T>::initialize(const Structure_Factor* sf,         //< structure factor
-                                           const ModulePW::PW_Basis_K* pw_wfc, //< planewave basis
-                                           const UnitCell* p_ucell,            //< unit cell
-                                           const K_Vectors* p_kv_in,
-                                           const int& random_seed,       //< random seed
-                                           const pseudopot_cell_vnl* p_pspot_nl,
-                                           const int& rank)
+void psi_init_atomic<T>::prepare_params(const int& nqx,
+                                        const double& dq,
+                                        const int& nspin,
+                                        const bool& domag,
+                                        const bool& domag_z,
+                                        const bool& pseudo_mesh,
+                                        const Structure_Factor& sf)
+{
+    this->nqx_ = nqx;
+    this->dq_ = dq;
+    this->nspin_ = nspin;
+    this->domag_ = domag;
+    this->domag_z_ = domag_z;
+    this->pseudo_mesh_ = pseudo_mesh;
+    this->sf_ = &sf;
+    this->params_prepared_ = true;
+}
+
+template <typename T>
+void psi_init_atomic<T>::initialize(const ModulePW::PW_Basis_K* pw_wfc,
+    const UnitCell* p_ucell,
+    const std::vector<int>& ik2iktot,
+    const int& nkstot,
+    const int& random_seed,
+    const int& lmaxkb,
+    const int& rank,
+    const int& npol,
+    const int& nbands)
 {
     ModuleBase::timer::start("psi_init_atomic", "initialize");
 
-    if(p_pspot_nl == nullptr)
+    if (!this->params_prepared_)
     {
-        ModuleBase::WARNING_QUIT("psi_init_atomic<T>::initialize", 
-                                 "pseudopot_cell_vnl object cannot be nullptr for atomic, quit.");
+        ModuleBase::WARNING_QUIT("psi_init_atomic<T>::initialize",
+            "prepare_params() must be called before initialize()");
     }
-    // import
-    psi_initializer<T>::initialize(sf, pw_wfc, p_ucell, p_kv_in, random_seed, p_pspot_nl, rank);
-    this->nbands_start_ = std::max(this->p_ucell_->natomwfc, PARAM.inp.nbands);
+
+    psi_base<T>::initialize(pw_wfc, p_ucell, ik2iktot, nkstot, random_seed, lmaxkb, rank, npol, nbands);
+
+    this->nbands_start_ = std::max(this->p_ucell_->natomwfc, nbands);
     this->nbands_complem_ = this->nbands_start_ - this->p_ucell_->natomwfc;
+
     // allocate
     this->allocate_ps_table();
+
     // then for generate random number to fill in the wavefunction
     this->ixy2is_.clear();
     this->ixy2is_.resize(this->pw_wfc_->fftnxy);
@@ -90,7 +103,7 @@ void psi_init_atomic<T>::tabulate()
     {
         max_msh = (this->p_ucell_->atoms[it].ncpp.msh > max_msh) ? this->p_ucell_->atoms[it].ncpp.msh : max_msh;
     }
-	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"max mesh points in Pseudopotential",max_msh);
+    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"max mesh points in Pseudopotential",max_msh);
     
     this->ovlp_pswfcjlq_.zero_out();
     const int startq = 0;
@@ -98,17 +111,17 @@ void psi_init_atomic<T>::tabulate()
     std::vector<double> aux(max_msh);
     std::vector<double> vchi(max_msh);
 
-	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"dq(describe PAO in reciprocal space)",PARAM.globalv.dq);
-	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"max q",PARAM.globalv.nqx);
+    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"dq(describe PAO in reciprocal space)",this->dq_);
+    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"max q",this->nqx_);
 
     for (int it=0; it<this->p_ucell_->ntype; it++)
     {
-		Atom* atom = &this->p_ucell_->atoms[it];
+	Atom* atom = &this->p_ucell_->atoms[it];
 
-		GlobalV::ofs_running<<"\n number of pseudo atomic orbitals for "<<atom->label<<" is "<< atom->ncpp.nchi << std::endl;
+	GlobalV::ofs_running<<"\n number of pseudo atomic orbitals for "<<atom->label<<" is "<< atom->ncpp.nchi << std::endl;
 
         // QE uses atom->ncpp.mesh
-        const int n_rgrid = (PARAM.inp.pseudo_mesh) ? atom->ncpp.mesh : atom->ncpp.msh;
+        const int n_rgrid = (this->pseudo_mesh_) ? atom->ncpp.mesh : atom->ncpp.msh;
         std::vector<double> chi2(n_rgrid);
 
         for (int ic = 0; ic < atom->ncpp.nchi ;ic++)
@@ -201,9 +214,9 @@ void psi_init_atomic<T>::tabulate()
             }
 
             const int l = atom->ncpp.lchi[ic];
-            for (int iq = startq; iq < PARAM.globalv.nqx; iq++)
+            for (int iq = startq; iq < this->nqx_; iq++)
             {
-                const double q = PARAM.globalv.dq * iq;
+                const double q = this->dq_ * iq;
                 ModuleBase::Sphbes::Spherical_Bessel(atom->ncpp.msh, atom->ncpp.r.data(), q, l, aux.data());
                 for (int ir = 0; ir < atom->ncpp.msh; ir++)
                 {
@@ -232,12 +245,13 @@ template <typename T>
 void psi_init_atomic<T>::init_psig(T* psig,  const int& ik)
 {
     ModuleBase::timer::start("psi_init_atomic", "init_psig");
+
     const int npw = this->pw_wfc_->npwk[ik];
     const int npwk_max = this->pw_wfc_->npwk_max;
     int lmax = this->p_ucell_->lmax_ppwf;
     const int total_lm = (lmax + 1) * (lmax + 1);
     ModuleBase::matrix ylm(total_lm, npw);
-    ModuleBase::GlobalFunc::ZEROS(psig, PARAM.globalv.npol * this->nbands_start_ * npwk_max);
+    ModuleBase::GlobalFunc::ZEROS(psig, this->npol_ * this->nbands_start_ * npwk_max);
 
     std::vector<std::complex<double>> aux(npw);
     std::vector<double> chiaux(npw);
@@ -273,17 +287,18 @@ void psi_init_atomic<T>::init_psig(T* psig,  const int& ik)
                     {
                         ovlp_pswfcjlg[ig] = ModuleBase::PolyInt::Polynomial_Interpolation(
                             this->ovlp_pswfcjlq_, it, ipswfc, 
-                            PARAM.globalv.nqx, PARAM.globalv.dq, gk[ig].norm() * this->p_ucell_->tpiba );
+                            this->nqx_, this->dq_, gk[ig].norm() * this->p_ucell_->tpiba );
                     }
 /* NSPIN == 4 */
-                    if(PARAM.inp.nspin == 4)
+                    if(this->nspin_ == 4)
                     {
                         if(this->p_ucell_->atoms[it].ncpp.has_so)
                         {
-                            Soc soc; soc.rot_ylm(l + 1);
+                            SocBase soc; 
+			    soc.rot_ylm(l + 1);
                             const double j = this->p_ucell_->atoms[it].ncpp.jchi[ipswfc];
     /* NOT NONCOLINEAR CASE, rotation matrix become identity */
-                            if (!(PARAM.globalv.domag||PARAM.globalv.domag_z))
+                            if (!(this->domag_||this->domag_z_))
                             {
                                 double cg_coeffs[2];
                                 for(int m = -l-1; m < l+1; m++)
@@ -297,7 +312,7 @@ void psi_init_atomic<T>::init_psig(T* psig,  const int& ik)
                                             if(fabs(cg_coeffs[is]) > 1e-8)
                                             {
         /* GET COMPLEX SPHERICAL HARMONIC FUNCTION */
-                                                const int ind = this->p_pspot_nl_->lmaxkb + soc.sph_ind(l,j,m,is); // ind can be l+m, l+m+1, l+m-1
+                                                const int ind = this->lmaxkb_ + soc.sph_ind(l,j,m,is); // ind can be l+m, l+m+1, l+m-1
                                                 std::fill(aux.begin(), aux.end(), std::complex<double>(0.0, 0.0));
                                                 for(int n1 = 0; n1 < 2*l+1; n1++)
                                                 {
@@ -366,7 +381,7 @@ void psi_init_atomic<T>::init_psig(T* psig,  const int& ik)
                                         chiaux[ig] =  l *
                                             ModuleBase::PolyInt::Polynomial_Interpolation(
                                                 this->ovlp_pswfcjlq_, it, ipswfc_noncolin_soc, 
-                                                PARAM.globalv.nqx, PARAM.globalv.dq, gk[ig].norm() * this->p_ucell_->tpiba);
+                                                this->nqx_, this->dq_, gk[ig].norm() * this->p_ucell_->tpiba);
                                         chiaux[ig] += ovlp_pswfcjlg[ig] * (l + 1.0) ;
                                         chiaux[ig] *= 1/(2.0*l+1.0);
                                     }
@@ -487,4 +502,4 @@ template class psi_init_atomic<std::complex<double>>;
 template class psi_init_atomic<std::complex<float>>;
 // gamma point calculation
 template class psi_init_atomic<double>;
-template class psi_init_atomic<float>;
\ No newline at end of file
+template class psi_init_atomic<float>;
diff --git a/source/source_psi/psi_init_atomic.h b/source/source_psi/psi_init_atomic.h
index 4cdfabdc6cd..53a4a0e4635 100644
--- a/source/source_psi/psi_init_atomic.h
+++ b/source/source_psi/psi_init_atomic.h
@@ -1,16 +1,30 @@
 #ifndef PSI_INIT_ATOMIC_H
 #define PSI_INIT_ATOMIC_H
+#include <vector>
+#include <string>
+#include <complex>
 #include "source_base/realarray.h"
-#include "psi_initializer.h"
+#include "psi_base.h"
+#include "source_basis/module_pw/pw_basis_k.h"
+#include "source_cell/unitcell.h"
+#include "source_pw/module_pwdft/structure_factor.h"
 
 /*
 Psi (planewave based wavefunction) initializer: atomic
 */
 template <typename T>
-class psi_init_atomic : public psi_initializer<T>
+class psi_init_atomic : public psi_base<T>
 {
   private:
     using Real = typename GetTypeReal<T>::type;
+    int nqx_ = 0;
+    double dq_ = 0.0;
+    int nspin_ = 1;
+    bool domag_ = false;
+    bool domag_z_ = false;
+    bool pseudo_mesh_ = false;
+    bool params_prepared_ = false;
+    const Structure_Factor* sf_ = nullptr;
 
   public:
     psi_init_atomic()
@@ -19,21 +33,66 @@ class psi_init_atomic : public psi_initializer<T>
     }
     ~psi_init_atomic(){};
 
-    /// @brief initialize the psi_init with external data and methods
-    virtual void initialize(const Structure_Factor*,             //< structure factor
-                            const ModulePW::PW_Basis_K*,         //< planewave basis
-                            const UnitCell*,                     //< unit cell
-                            const K_Vectors*,                    //< kpoints
-                            const int& = 1,                      //< random seed
-                            const pseudopot_cell_vnl* = nullptr, //< nonlocal pseudopotential
-                            const int& = 0) override;            //< MPI rank
+    /**
+     * @brief Prepare parameters before initialization.
+     * 
+     * This method must be called before initialize(). It sets up the necessary
+     * parameters for the psi initialization process.
+     * 
+     * @param nqx Number of q-points for interpolation
+     * @param dq Spacing between q-points
+     * @param nspin Number of spin components
+     * @param domag Whether to use non-collinear magnetism
+     * @param domag_z Whether to use z-axis only non-collinear magnetism
+     * @param pseudo_mesh Whether to use pseudo mesh for radial grid
+     * 
+     * @see initialize()
+     */
+    void prepare_params(const int& nqx,
+                        const double& dq,
+                        const int& nspin,
+                        const bool& domag,
+                        const bool& domag_z,
+                        const bool& pseudo_mesh,
+                        const Structure_Factor& sf);
+
+    /**
+     * @brief Initialize the psi_init with external data and methods.
+     * 
+     * This method must be called after prepare_params(). It initializes the
+     * psi initializer with the provided planewave basis and unit cell information.
+     * 
+     * @param pw_wfc Planewave basis
+     * @param p_ucell Unit cell
+     * @param ik2iktot Local->global k-point mapping
+     * @param nkstot Total number of k-points
+     * @param random_seed Random seed
+     * @param lmaxkb Max angular momentum for non-local projectors
+     * @param rank MPI rank
+     * @param npol Number of polarization components
+     * @param nbands Number of bands
+     * 
+     * @see prepare_params()
+     */
+    virtual void initialize(const ModulePW::PW_Basis_K* pw_wfc,         //< planewave basis
+                            const UnitCell* p_ucell,                     //< unit cell
+                            const std::vector<int>& ik2iktot,             //< ik2iktot: local->global k-point mapping
+                            const int& nkstot,                          //< nkstot: total number of k-points
+                            const int& random_seed,                      //< random seed
+                            const int& lmaxkb,                          //< lmaxkb: max angular momentum for non-local projectors
+                            const int& rank,                            //< MPI rank
+                            const int& npol,                            //< npol
+                            const int& nbands) override;                //< nbands
     virtual void tabulate() override;
     virtual void init_psig(T* psig, const int& ik) override;
 
   protected:
+
     // allocate memory for overlap table
     void allocate_ps_table();
+
     std::vector<std::string> pseudopot_files_;
+
     ModuleBase::realArray ovlp_pswfcjlq_;
 };
-#endif
\ No newline at end of file
+#endif
diff --git a/source/source_psi/psi_init_atomic_random.cpp b/source/source_psi/psi_init_atomic_random.cpp
index 91fd4c10bc8..d7b51c8d3c2 100644
--- a/source/source_psi/psi_init_atomic_random.cpp
+++ b/source/source_psi/psi_init_atomic_random.cpp
@@ -3,15 +3,17 @@
 #include "source_io/module_parameter/parameter.h"
 
 template <typename T>
-void psi_init_atomic_random<T>::initialize(const Structure_Factor* sf,         //< structure factor
-                                                  const ModulePW::PW_Basis_K* pw_wfc, //< planewave basis
+void psi_init_atomic_random<T>::initialize(const ModulePW::PW_Basis_K* pw_wfc, //< planewave basis
                                                   const UnitCell* p_ucell,            //< unit cell
-                                                  const K_Vectors* p_kv_in,
+                                                  const std::vector<int>& ik2iktot,
+                                                  const int& nkstot,
                                                   const int& random_seed, //< random seed
-                                                  const pseudopot_cell_vnl* p_pspot_nl,
-                                                  const int& rank)
+                                                  const int& lmaxkb,
+                                                  const int& rank,
+                                                  const int& npol,
+                                                  const int& nbands)
 {
-    psi_init_atomic<T>::initialize(sf, pw_wfc, p_ucell, p_kv_in, random_seed, p_pspot_nl, rank);
+    psi_init_atomic<T>::initialize(pw_wfc, p_ucell, ik2iktot, nkstot, random_seed, lmaxkb, rank, npol, nbands);
 }
 
 template <typename T>
@@ -19,7 +21,7 @@ void psi_init_atomic_random<T>::init_psig(T* psig, const int& ik)
 {
     double rm = this->mixing_coef_;
     psi_init_atomic<T>::init_psig(psig, ik);
-    const int npol = PARAM.globalv.npol;
+    const int npol = this->npol_;
     const int nbasis = this->pw_wfc_->npwk_max * npol;
     psi::Psi<T> psi_random(1, this->nbands_start_, nbasis, nbasis, true);
     psi_random.fix_k(0);
diff --git a/source/source_psi/psi_init_atomic_random.h b/source/source_psi/psi_init_atomic_random.h
index 2c8e49fc8d6..b43c633d5e5 100644
--- a/source/source_psi/psi_init_atomic_random.h
+++ b/source/source_psi/psi_init_atomic_random.h
@@ -1,6 +1,5 @@
 #ifndef PSI_INIT_ATOMIC_RANDOM_H
 #define PSI_INIT_ATOMIC_RANDOM_H
-#include "source_pw/module_pwdft/vnl_pw.h"
 #include "psi_init_atomic.h"
 
 /*
@@ -20,17 +19,36 @@ class psi_init_atomic_random : public psi_init_atomic<T>
     }
     ~psi_init_atomic_random(){};
 
-    /// @brief initialize the psi_initializer with external data and methods
-    virtual void initialize(const Structure_Factor*,             //< structure factor
-                            const ModulePW::PW_Basis_K*,         //< planewave basis
-                            const UnitCell*,                     //< unit cell
-                            const K_Vectors*,                    //< kpoints
-                            const int& = 1,                      //< random seed
-                            const pseudopot_cell_vnl* = nullptr, //< nonlocal pseudopotential
-                            const int& = 0) override;            //< MPI rank
+    /**
+     * @brief Initialize the psi_init with external data and methods.
+     * 
+     * This method must be called after prepare_params(). It initializes the
+     * psi initializer with the provided planewave basis and unit cell information.
+     * 
+     * @param pw_wfc Planewave basis
+     * @param p_ucell Unit cell
+     * @param ik2iktot Local->global k-point mapping
+     * @param nkstot Total number of k-points
+     * @param random_seed Random seed
+     * @param lmaxkb Max angular momentum for non-local projectors
+     * @param rank MPI rank
+     * @param npol Number of polarization components
+     * @param nbands Number of bands
+     * 
+     * @see prepare_params()
+     */
+    virtual void initialize(const ModulePW::PW_Basis_K* pw_wfc,         //< planewave basis
+                            const UnitCell* p_ucell,                     //< unit cell
+                            const std::vector<int>& ik2iktot,             //< ik2iktot: local->global k-point mapping
+                            const int& nkstot,                          //< nkstot: total number of k-points
+                            const int& random_seed,                      //< random seed
+                            const int& lmaxkb,                          //< lmaxkb: max angular momentum for non-local projectors
+                            const int& rank,                            //< MPI rank
+                            const int& npol,                            //< npol
+                            const int& nbands) override;                //< nbands
 
     virtual void init_psig(T* psig, const int& ik) override;
 
   private:
 };
-#endif
\ No newline at end of file
+#endif
diff --git a/source/source_psi/psi_init_file.cpp b/source/source_psi/psi_init_file.cpp
index 2e633573bb4..607bf800c03 100644
--- a/source/source_psi/psi_init_file.cpp
+++ b/source/source_psi/psi_init_file.cpp
@@ -1,48 +1,65 @@
 #include "psi_init_file.h"
 
+#include <vector>
+#include <cassert>
+#include <string>
+#include <complex>
+
 #include "source_base/timer.h"
-#include "source_cell/klist.h"
 #include "source_io/module_wf/read_wfc_pw.h"
 #include "source_io/module_output/filename.h"
-#include "source_io/module_parameter/parameter.h"
 
 template <typename T>
-void psi_init_file<T>::initialize(const Structure_Factor* sf,
-                                               const ModulePW::PW_Basis_K* pw_wfc,
+void psi_init_file<T>::initialize(const ModulePW::PW_Basis_K* pw_wfc,
                                                const UnitCell* p_ucell,
-                                               const K_Vectors* p_kv_in,
+                                               const std::vector<int>& ik2iktot,
+                                               const int& nkstot,
                                                const int& random_seed,
-                                               const pseudopot_cell_vnl* p_pspot_nl,
-                                               const int& rank)
+                                               const int& lmaxkb,
+                                               const int& rank,
+                                               const int& npol,
+                                               const int& nbands)
 {
-    psi_initializer<T>::initialize(sf, pw_wfc, p_ucell, p_kv_in, random_seed, p_pspot_nl, rank);
-    this->nbands_start_ = PARAM.inp.nbands;
+    psi_base<T>::initialize(pw_wfc, p_ucell, ik2iktot, nkstot, random_seed, lmaxkb, rank, npol, nbands);
+    this->nbands_start_ = nbands;
     this->nbands_complem_ = 0;
 }
 
+template <typename T>
+void psi_init_file<T>::prepare_params(const int& nspin,
+                                      const std::string& global_readin_dir,
+                                      const int& rank_in_pool,
+                                      const int& nproc_in_pool)
+{
+    this->nspin_ = nspin;
+    this->global_readin_dir_ = global_readin_dir;
+    this->rank_in_pool_ = rank_in_pool;
+    this->nproc_in_pool_ = nproc_in_pool;
+}
+
 template <typename T>
 void psi_init_file<T>::init_psig(T* psig, const int& ik)
 {
     ModuleBase::timer::start("psi_init_file", "init_psig");
-    const int npol = PARAM.globalv.npol;
+    const int npol = this->npol_;
     const int nbasis = this->pw_wfc_->npwk_max * npol;
-    const int nkstot = this->p_kv->get_nkstot();
+    const int nkstot = this->nkstot_;
     ModuleBase::ComplexMatrix wfcatom(this->nbands_start_, nbasis);
-    int ik_tot = this->p_kv->ik2iktot[ik];
+    int ik_tot = this->ik2iktot_[ik];
 
     // mohan update, this is for plane wave, 2025-05-17
-	const int out_type = 2;
-	const bool out_app_flag = false;
-	const bool gamma_only = false;
-	const int istep = -1;
+    const int out_type = 2;
+    const bool out_app_flag = false;
+    const bool gamma_only = false;
+    const int istep = -1;
 
-	std::string fn = ModuleIO::filename_output(PARAM.globalv.global_readin_dir,"wf","pw",
-			ik,this->p_kv->ik2iktot,PARAM.inp.nspin,nkstot,
+    std::string fn = ModuleIO::filename_output(this->global_readin_dir_,"wf","pw",
+			ik,this->ik2iktot_,this->nspin_,nkstot,
 			out_type,out_app_flag,gamma_only,istep);
 
-	ModuleIO::read_wfc_pw(fn, this->pw_wfc_, 
-			GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL,
-			PARAM.inp.nbands, PARAM.globalv.npol,
+    ModuleIO::read_wfc_pw(fn, this->pw_wfc_, 
+			this->rank_in_pool_, this->nproc_in_pool_,
+			this->nbands_start_, this->npol_,
 			ik, ik_tot, nkstot, wfcatom);
 
     assert(this->nbands_start_ <= wfcatom.nr);
diff --git a/source/source_psi/psi_init_file.h b/source/source_psi/psi_init_file.h
index 72fb18ed1eb..f8b8ebd796b 100644
--- a/source/source_psi/psi_init_file.h
+++ b/source/source_psi/psi_init_file.h
@@ -1,17 +1,22 @@
 #ifndef PSI_INIT_FILE_H
 #define PSI_INIT_FILE_H
 
-#include "source_pw/module_pwdft/vnl_pw.h"
-#include "psi_initializer.h"
+#include <vector>
+#include <string>
+#include "psi_base.h"
 
 /*
 Psi (planewave based wavefunction) initializer: random method
 */
 template <typename T>
-class psi_init_file : public psi_initializer<T>
+class psi_init_file : public psi_base<T>
 {
   private:
     using Real = typename GetTypeReal<T>::type;
+    int nspin_ = 1;
+    std::string global_readin_dir_;
+    int rank_in_pool_ = 0;
+    int nproc_in_pool_ = 1;
 
   public:
     psi_init_file()
@@ -20,18 +25,25 @@ class psi_init_file : public psi_initializer<T>
     };
     ~psi_init_file(){};
 
-    /// @brief initialize the psi_initializer with external data and methods
-    virtual void initialize(const Structure_Factor*,             //< structure factor
-                            const ModulePW::PW_Basis_K*,         //< planewave basis
-                            const UnitCell*,                     //< unit cell
-                            const K_Vectors*,                    //< kpoints
-                            const int& = 1,                      //< random seed
-                            const pseudopot_cell_vnl* = nullptr, //< nonlocal pseudopotential
-                            const int& = 0) override;            //< MPI rank
+    /// @brief initialize the psi_base with external data and methods
+    virtual void initialize(const ModulePW::PW_Basis_K* pw_wfc,         //< planewave basis
+                            const UnitCell* p_ucell,                     //< unit cell
+                            const std::vector<int>& ik2iktot,             //< ik2iktot: local->global k-point mapping
+                            const int& nkstot,                          //< nkstot: total number of k-points
+                            const int& random_seed,                      //< random seed
+                            const int& lmaxkb,                          //< lmaxkb: max angular momentum for non-local projectors
+                            const int& rank,                            //< MPI rank
+                            const int& npol,                            //< npol
+                            const int& nbands) override;                //< nbands
 
     /// @brief calculate and output planewave wavefunction
     /// @param ik kpoint index
     /// @return initialized planewave wavefunction (psi::Psi<std::complex<double>>*)
     virtual void init_psig(T* psig, const int& ik) override;
+
+    void prepare_params(const int& nspin,
+                        const std::string& global_readin_dir,
+                        const int& rank_in_pool,
+                        const int& nproc_in_pool);
 };
 #endif
diff --git a/source/source_psi/psi_init_nao.cpp b/source/source_psi/psi_init_nao.cpp
index 916d87295a9..57008cc5265 100644
--- a/source/source_psi/psi_init_nao.cpp
+++ b/source/source_psi/psi_init_nao.cpp
@@ -17,9 +17,7 @@
 #include "source_base/parallel_reduce.h"
 #endif
 #include "source_io/module_output/orb_io.h"
-#include "source_io/module_parameter/parameter.h"
-// GlobalV::NQX and GlobalV::DQ are here
-#include "source_io/module_parameter/parameter.h"
+
 
 #include <algorithm>
 #include <numeric>
@@ -42,6 +40,21 @@ void normalize(const std::vector<double>& r, std::vector<double>& flz)
     std::transform(flz.begin(), flz.end(), flz.begin(), [norm](double flz) { return flz / norm; });
 }
 
+template <typename T>
+void psi_init_nao<T>::prepare_params(const int& nqx,
+                                     const double& dq,
+                                     const int& nspin,
+                                     const std::string& orbital_dir,
+                                     const Structure_Factor& sf)
+{
+    this->nqx_ = nqx;
+    this->dq_ = dq;
+    this->nspin_ = nspin;
+    this->orbital_dir_ = orbital_dir;
+    this->sf_ = &sf;
+    this->params_prepared_ = true;
+}
+
 template <typename T>
 void psi_init_nao<T>::read_external_orbs(const std::string* orbital_files, const int& rank)
 {
@@ -67,7 +80,7 @@ void psi_init_nao<T>::read_external_orbs(const std::string* orbital_files, const
         bool is_open = false;
         if (rank == 0)
         {
-            ifs_it.open(PARAM.inp.orbital_dir + this->orbital_files_[it]);
+            ifs_it.open(this->orbital_dir_ + this->orbital_files_[it]);
             is_open = ifs_it.is_open();
         }
 #ifdef __MPI
@@ -147,18 +160,26 @@ void psi_init_nao<T>::allocate_ao_table()
 }
 
 template <typename T>
-void psi_init_nao<T>::initialize(const Structure_Factor* sf,
-                                        const ModulePW::PW_Basis_K* pw_wfc,
+void psi_init_nao<T>::initialize(const ModulePW::PW_Basis_K* pw_wfc,
                                         const UnitCell* p_ucell,
-                                        const K_Vectors* p_kv_in,
+                                        const std::vector<int>& ik2iktot,
+                                        const int& nkstot,
                                         const int& random_seed,
-                                        const pseudopot_cell_vnl* p_pspot_nl,
-                                        const int& rank)
+                                        const int& lmaxkb,
+                                        const int& rank,
+                                        const int& npol,
+                                        const int& nbands)
 {
     ModuleBase::timer::start("psi_init_nao", "initialize");
 
+    if (!this->params_prepared_)
+    {
+        ModuleBase::WARNING_QUIT("psi_init_nao<T>::initialize",
+            "prepare_params() must be called before initialize()");
+    }
+
     // import
-    psi_initializer<T>::initialize(sf, pw_wfc, p_ucell, p_kv_in, random_seed, p_pspot_nl, rank);
+    psi_base<T>::initialize(pw_wfc, p_ucell, ik2iktot, nkstot, random_seed, lmaxkb, rank, npol, nbands);
 
     // allocate
     this->allocate_ao_table();
@@ -177,18 +198,10 @@ void psi_init_nao<T>::initialize(const Structure_Factor* sf,
             /* EVERY ZETA FOR (2l+1) ORBS */
             const int nchi = this->p_ucell_->atoms[it].l_nchi[l];
             const int degen_l = (l == 0) ? 1 : 2 * l + 1;
-            nbands_local += nchi * degen_l * PARAM.globalv.npol * this->p_ucell_->atoms[it].na;
-            /*
-                non-rotate basis, nbands_local*=2 for PARAM.globalv.npol = 2 is enough
-            */
-            // nbands_local += this->p_ucell_->atoms[it].l_nchi[l]*(2*l+1) * PARAM.globalv.npol;
-            /*
-                rotate basis, nbands_local*=4 for p, d, f,... orbitals, and nbands_local*=2 for s orbitals
-                risky when NSPIN = 4, problematic psi value, needed to be checked
-            */
+            nbands_local += nchi * degen_l * npol * this->p_ucell_->atoms[it].na;
         }
     }
-    this->nbands_start_ = std::max(nbands_local, PARAM.inp.nbands);
+    this->nbands_start_ = std::max(nbands_local, nbands);
     this->nbands_complem_ = this->nbands_start_ - nbands_local;
 
     ModuleBase::timer::end("psi_init_nao", "initialize");
@@ -199,10 +212,16 @@ void psi_init_nao<T>::tabulate()
 {
     ModuleBase::timer::start("psi_init_nao", "tabulate");
 
+    if (this->nqx_ <= 0)
+    {
+        ModuleBase::WARNING_QUIT("psi_init_nao<T>::tabulate",
+            "nqx_ must be greater than 0. Did you forget to call prepare_params() with valid nqx?");
+    }
+
     // a uniformed qgrid
-    std::vector<double> qgrid(PARAM.globalv.nqx);
+    std::vector<double> qgrid(this->nqx_);
     std::iota(qgrid.begin(), qgrid.end(), 0);
-    std::for_each(qgrid.begin(), qgrid.end(), [this](double& q) { q = q * PARAM.globalv.dq; });
+    std::for_each(qgrid.begin(), qgrid.end(), [this](double& q) { q = q * this->dq_; });
 
     // only when needed, allocate memory for cubspl_
     if (this->cubspl_.get())
@@ -224,7 +243,7 @@ void psi_init_nao<T>::tabulate()
     ModuleBase::SphericalBesselTransformer sbt_(true); // bool: enable cache
 
     // tabulate the spherical bessel transform of numerical orbital function
-    std::vector<double> Jlfq(PARAM.globalv.nqx, 0.0);
+    std::vector<double> Jlfq(this->nqx_, 0.0);
     int i = 0;
     for (int it = 0; it < this->p_ucell_->ntype; it++)
     {
@@ -237,7 +256,7 @@ void psi_init_nao<T>::tabulate()
                             this->nr_[it][ic],
                             this->rgrid_[it][ic].data(),
                             this->chi_[it][ic].data(),
-                            PARAM.globalv.nqx,
+                            this->nqx_,
                             qgrid.data(),
                             Jlfq.data());
                 this->cubspl_->add(Jlfq.data());
@@ -258,7 +277,7 @@ void psi_init_nao<T>::init_psig(T* psig, const int& ik)
     const int npwk_max = this->pw_wfc_->npwk_max;
     const int total_lm = (this->p_ucell_->lmax + 1) * (this->p_ucell_->lmax + 1);
     ModuleBase::matrix ylm(total_lm, npw);
-    ModuleBase::GlobalFunc::ZEROS(psig, PARAM.globalv.npol * this->nbands_start_ * npwk_max);
+    ModuleBase::GlobalFunc::ZEROS(psig, this->npol_ * this->nbands_start_ * npwk_max);
 
     std::vector<std::complex<double>> aux(npw);
     std::vector<double> qnorm(npw);
@@ -299,7 +318,7 @@ void psi_init_nao<T>::init_psig(T* psig, const int& ik)
                     this->cubspl_->eval(npw, qnorm.data(), Jlfq.data(), nullptr, nullptr, this->projmap_(it, L, N));
 
                     /* FOR EVERY NAO IN EACH ATOM */
-                    if (PARAM.inp.nspin == 4)
+                    if (this->nspin_ == 4)
                     {
                         /* FOR EACH SPIN CHANNEL */
                         for (int is_N = 0; is_N < 2; is_N++) // rotate base
diff --git a/source/source_psi/psi_init_nao.h b/source/source_psi/psi_init_nao.h
index bb05962ec07..5c98b9713c3 100644
--- a/source/source_psi/psi_init_nao.h
+++ b/source/source_psi/psi_init_nao.h
@@ -3,17 +3,29 @@
 #include "source_base/cubic_spline.h"
 #include "source_base/realarray.h"
 #include "source_base/spherical_bessel_transformer.h"
-#include "psi_initializer.h"
+#include "psi_base.h"
+#include "source_basis/module_pw/pw_basis_k.h"
+#include "source_cell/unitcell.h"
+#include "source_pw/module_pwdft/structure_factor.h"
 
 #include <memory>
+#include <string>
+#include <vector>
+#include <complex>
 /*
 Psi (planewave based wavefunction) initializer: numerical atomic orbital method
 */
 template <typename T>
-class psi_init_nao : public psi_initializer<T>
+class psi_init_nao : public psi_base<T>
 {
   private:
     using Real = typename GetTypeReal<T>::type;
+    int nqx_ = 0;
+    double dq_ = 0.0;
+    int nspin_ = 1;
+    std::string orbital_dir_;
+    bool params_prepared_ = false;
+    const Structure_Factor* sf_ = nullptr;
 
   public:
     psi_init_nao()
@@ -22,63 +34,116 @@ class psi_init_nao : public psi_initializer<T>
     };
     ~psi_init_nao(){};
 
-    virtual void init_psig(T* psig, const int& ik) override;
+    /**
+     * @brief Prepare parameters before initialization.
+     * 
+     * This method must be called before initialize(). It sets up the necessary
+     * parameters for the psi initialization process.
+     * 
+     * @param nqx Number of q-points for interpolation
+     * @param dq Spacing between q-points
+     * @param nspin Number of spin components
+     * @param orbital_dir Directory containing orbital files
+     * 
+     * @see initialize()
+     */
+    void prepare_params(const int& nqx,
+                        const double& dq,
+                        const int& nspin,
+                        const std::string& orbital_dir,
+                        const Structure_Factor& sf);
 
-    /// @brief initialize the psi_initializer with external data and methods
-    virtual void initialize(const Structure_Factor*,             //< structure factor
-                            const ModulePW::PW_Basis_K*,         //< planewave basis
-                            const UnitCell*,                     //< unit cell
-                            const K_Vectors*,                    //< kpoints
-                            const int& = 1,                      //< random seed
-                            const pseudopot_cell_vnl* = nullptr, //< nonlocal pseudopotential
-                            const int& = 0) override;            //< MPI rank
+    /**
+     * @brief Initialize the psi_init with external data and methods.
+     * 
+     * This method must be called after prepare_params(). It initializes the
+     * psi initializer with the provided structure factor, planewave basis,
+     * and unit cell information.
+     * 
+     * @param sf Structure factor
+     * @param pw_wfc Planewave basis
+     * @param p_ucell Unit cell
+     * @param ik2iktot Local->global k-point mapping
+     * @param nkstot Total number of k-points
+     * @param random_seed Random seed
+     * @param lmaxkb Max angular momentum for non-local projectors
+     * @param rank MPI rank
+     * @param npol Number of polarization components
+     * @param nbands Number of bands
+     * 
+     * @see prepare_params()
+     */
+    virtual void initialize(const ModulePW::PW_Basis_K* pw_wfc,         //< planewave basis
+                            const UnitCell* p_ucell,                     //< unit cell
+                            const std::vector<int>& ik2iktot,             //< ik2iktot: local->global k-point mapping
+                            const int& nkstot,                          //< nkstot: total number of k-points
+                            const int& random_seed,                      //< random seed
+                            const int& lmaxkb,                          //< lmaxkb: max angular momentum for non-local projectors
+                            const int& rank,                            //< MPI rank
+                            const int& npol,                            //< npol
+                            const int& nbands) override;                //< nbands
 
     void read_external_orbs(const std::string* orbital_files, const int& rank);
+
     virtual void tabulate() override;
+
+    virtual void init_psig(T* psig, const int& ik) override;
+
     std::vector<std::string> external_orbs() const
     {
         return orbital_files_;
     }
+
     std::vector<std::vector<int>> nr() const
     {
         return nr_;
     }
+
     std::vector<int> nr(const int& itype) const
     {
         return nr_[itype];
     }
+
     int nr(const int& itype, const int& ichi) const
     {
         return nr_[itype][ichi];
     }
+
     std::vector<std::vector<std::vector<double>>> chi() const
     {
         return chi_;
     }
+
     std::vector<std::vector<double>> chi(const int& itype) const
     {
         return chi_[itype];
     }
+
     std::vector<double> chi(const int& itype, const int& ichi) const
     {
         return chi_[itype][ichi];
     }
+
     double chi(const int& itype, const int& ichi, const int& ir) const
     {
         return chi_[itype][ichi][ir];
     }
+
     std::vector<std::vector<std::vector<double>>> rgrid() const
     {
         return rgrid_;
     }
+
     std::vector<std::vector<double>> rgrid(const int& itype) const
     {
         return rgrid_[itype];
     }
+
     std::vector<double> rgrid(const int& itype, const int& ichi) const
     {
         return rgrid_[itype][ichi];
     }
+
     double rgrid(const int& itype, const int& ichi, const int& ir) const
     {
         return rgrid_[itype][ichi][ir];
@@ -101,4 +166,4 @@ class psi_init_nao : public psi_initializer<T>
     /// @brief useful for atomic-like methods
     ModuleBase::SphericalBesselTransformer sbt;
 };
-#endif
\ No newline at end of file
+#endif
diff --git a/source/source_psi/psi_init_nao_random.cpp b/source/source_psi/psi_init_nao_random.cpp
index e3e2b8e89cd..aaabbfd8c3a 100644
--- a/source/source_psi/psi_init_nao_random.cpp
+++ b/source/source_psi/psi_init_nao_random.cpp
@@ -3,15 +3,17 @@
 #include "source_io/module_parameter/parameter.h"
 
 template <typename T>
-void psi_init_nao_random<T>::initialize(const Structure_Factor* sf,
-                                               const ModulePW::PW_Basis_K* pw_wfc,
+void psi_init_nao_random<T>::initialize(const ModulePW::PW_Basis_K* pw_wfc,
                                                const UnitCell* p_ucell,
-                                               const K_Vectors* p_kv_in,
+                                               const std::vector<int>& ik2iktot,
+                                               const int& nkstot,
                                                const int& random_seed,
-                                               const pseudopot_cell_vnl* p_pspot_nl,
-                                               const int& rank)
+                                               const int& lmaxkb,
+                                               const int& rank,
+                                               const int& npol,
+                                               const int& nbands)
 {
-    psi_init_nao<T>::initialize(sf, pw_wfc, p_ucell, p_kv_in, random_seed, p_pspot_nl, rank);
+    psi_init_nao<T>::initialize(pw_wfc, p_ucell, ik2iktot, nkstot, random_seed, lmaxkb, rank, npol, nbands);
 }
 
 template <typename T>
@@ -19,7 +21,7 @@ void psi_init_nao_random<T>::init_psig(T* psig, const int& ik)
 {
     double rm = this->mixing_coef_;
     psi_init_nao<T>::init_psig(psig, ik);
-    const int npol = PARAM.globalv.npol;
+    const int npol = this->npol_;
     const int nbasis = this->pw_wfc_->npwk_max * npol;
     psi::Psi<T> psi_random(1, this->nbands_start_, nbasis, nbasis, true);
     psi_random.fix_k(0);
diff --git a/source/source_psi/psi_init_nao_random.h b/source/source_psi/psi_init_nao_random.h
index d6b38bb5568..cbe042e68be 100644
--- a/source/source_psi/psi_init_nao_random.h
+++ b/source/source_psi/psi_init_nao_random.h
@@ -1,6 +1,5 @@
 #ifndef PSI_INIT_NAO_RANDOM_H
 #define PSI_INIT_NAO_RANDOM_H
-#include "source_pw/module_pwdft/vnl_pw.h"
 #include "psi_init_nao.h"
 
 /*
@@ -20,14 +19,33 @@ class psi_init_nao_random : public psi_init_nao<T>
     };
     ~psi_init_nao_random(){};
 
-    /// @brief initialize the psi_init with external data and methods
-    virtual void initialize(const Structure_Factor*,             //< structure factor
-                            const ModulePW::PW_Basis_K*,         //< planewave basis
-                            const UnitCell*,                     //< unit cell
-                            const K_Vectors*,                    //< kpoints
-                            const int& = 1,                      //< random seed
-                            const pseudopot_cell_vnl* = nullptr, //< nonlocal pseudopotential
-                            const int& = 0) override;            //< MPI rank
+    /**
+     * @brief Initialize the psi_init with external data and methods.
+     * 
+     * This method must be called after prepare_params(). It initializes the
+     * psi initializer with the provided planewave basis and unit cell information.
+     * 
+     * @param pw_wfc Planewave basis
+     * @param p_ucell Unit cell
+     * @param ik2iktot Local->global k-point mapping
+     * @param nkstot Total number of k-points
+     * @param random_seed Random seed
+     * @param lmaxkb Max angular momentum for non-local projectors
+     * @param rank MPI rank
+     * @param npol Number of polarization components
+     * @param nbands Number of bands
+     * 
+     * @see prepare_params()
+     */
+    virtual void initialize(const ModulePW::PW_Basis_K* pw_wfc,         //< planewave basis
+                            const UnitCell* p_ucell,                     //< unit cell
+                            const std::vector<int>& ik2iktot,             //< ik2iktot: local->global k-point mapping
+                            const int& nkstot,                          //< nkstot: total number of k-points
+                            const int& random_seed,                      //< random seed
+                            const int& lmaxkb,                          //< lmaxkb: max angular momentum for non-local projectors
+                            const int& rank,                            //< MPI rank
+                            const int& npol,                            //< npol
+                            const int& nbands) override;                //< nbands
 
     virtual void init_psig(T* psig, const int& ik) override;
 };
diff --git a/source/source_psi/psi_init_random.cpp b/source/source_psi/psi_init_random.cpp
index 697a4476adb..b509a57d869 100644
--- a/source/source_psi/psi_init_random.cpp
+++ b/source/source_psi/psi_init_random.cpp
@@ -1,20 +1,22 @@
 #include "psi_init_random.h"
-#include "source_io/module_parameter/parameter.h"
+#include <vector>
 
 template <typename T>
-void psi_init_random<T>::initialize(const Structure_Factor* sf,
-                                               const ModulePW::PW_Basis_K* pw_wfc,
+void psi_init_random<T>::initialize(const ModulePW::PW_Basis_K* pw_wfc,
                                                const UnitCell* p_ucell,
-                                               const K_Vectors* p_kv_in,
+                                               const std::vector<int>& ik2iktot,
+                                               const int& nkstot,
                                                const int& random_seed,
-                                               const pseudopot_cell_vnl* p_pspot_nl,
-                                               const int& rank)
+                                               const int& lmaxkb,
+                                               const int& rank,
+                                               const int& npol,
+                                               const int& nbands)
 {
-    psi_initializer<T>::initialize(sf, pw_wfc, p_ucell, p_kv_in, random_seed, p_pspot_nl, rank);
+    psi_base<T>::initialize(pw_wfc, p_ucell, ik2iktot, nkstot, random_seed, lmaxkb, rank, npol, nbands);
     this->ixy2is_.clear();
     this->ixy2is_.resize(this->pw_wfc_->fftnxy);
     this->pw_wfc_->getfftixy2is(this->ixy2is_.data());
-    this->nbands_start_ = PARAM.inp.nbands;
+    this->nbands_start_ = nbands;
     this->nbands_complem_ = 0;
 }
 
diff --git a/source/source_psi/psi_init_random.h b/source/source_psi/psi_init_random.h
index 8d66035ec60..c4c9e979ac0 100644
--- a/source/source_psi/psi_init_random.h
+++ b/source/source_psi/psi_init_random.h
@@ -1,14 +1,15 @@
 #ifndef PSI_INIT_RANDOM_H
 #define PSI_INIT_RANDOM_H
 
+#include <vector>
 #include "source_pw/module_pwdft/vnl_pw.h"
-#include "psi_initializer.h"
+#include "psi_base.h"
 
 /*
 Psi (planewave based wavefunction) initializer: random method
 */
 template <typename T>
-class psi_init_random : public psi_initializer<T>
+class psi_init_random : public psi_base<T>
 {
   private:
     using Real = typename GetTypeReal<T>::type;
@@ -24,12 +25,14 @@ class psi_init_random : public psi_initializer<T>
     /// @return initialized planewave wavefunction (psi::Psi<std::complex<double>>*)
     virtual void init_psig(T* psig, const int& ik) override;
     /// @brief initialize the psi_init with external data and methods
-    virtual void initialize(const Structure_Factor*,             //< structure factor
-                            const ModulePW::PW_Basis_K*,         //< planewave basis
-                            const UnitCell*,                     //< unit cell
-                            const K_Vectors*,                    //< kpoints
-                            const int& = 1,                      //< random seed
-                            const pseudopot_cell_vnl* = nullptr, //< nonlocal pseudopotential
-                            const int& = 0) override;            //< MPI rank
+    virtual void initialize(const ModulePW::PW_Basis_K* pw_wfc,         //< planewave basis
+                            const UnitCell* p_ucell,                     //< unit cell
+                            const std::vector<int>& ik2iktot,             //< ik2iktot: local->global k-point mapping
+                            const int& nkstot,                          //< nkstot: total number of k-points
+                            const int& random_seed,                      //< random seed
+                            const int& lmaxkb,                          //< lmaxkb: max angular momentum for non-local projectors
+                            const int& rank,                            //< MPI rank
+                            const int& npol,                            //< npol
+                            const int& nbands) override;                //< nbands
 };
 #endif
\ No newline at end of file
diff --git a/source/source_psi/psi_prepare.cpp b/source/source_psi/psi_prepare.cpp
index adc10eeb3fa..3ab5c4d5756 100644
--- a/source/source_psi/psi_prepare.cpp
+++ b/source/source_psi/psi_prepare.cpp
@@ -7,13 +7,14 @@
 #include "source_base/timer.h"
 #include "source_base/tool_quit.h"
 #include "source_hsolver/diago_iter_assist.h"
-#include "source_io/module_parameter/parameter.h"
+
 #include "source_psi/psi_init_atomic.h"
 #include "source_psi/psi_init_atomic_random.h"
 #include "source_psi/psi_init_file.h"
 #include "source_psi/psi_init_nao.h"
 #include "source_psi/psi_init_nao_random.h"
 #include "source_psi/psi_init_random.h"
+
 namespace psi
 {
 
@@ -24,10 +25,11 @@ PSIPrepare<T, Device>::PSIPrepare(const std::string& init_wfc_in,
                             const int& rank_in,
                             const UnitCell& ucell_in,
                             const Structure_Factor& sf_in,
-                            const K_Vectors& kv_in,
-                            const pseudopot_cell_vnl& nlpp_in,
+                            const std::vector<int>& ik2iktot_in,
+                            const int& nkstot_in,
+                            const int& lmaxkb_in,
                             const ModulePW::PW_Basis_K& pw_wfc_in)
-    : ucell(ucell_in), sf(sf_in), nlpp(nlpp_in), kv(kv_in), pw_wfc(pw_wfc_in), rank(rank_in)
+    : ucell(ucell_in), sf(sf_in), lmaxkb(lmaxkb_in), pw_wfc(pw_wfc_in), rank(rank_in), ik2iktot_(ik2iktot_in), nkstot_(nkstot_in)
 {
     this->init_wfc = init_wfc_in;
     this->ks_solver = ks_solver_in;
@@ -35,7 +37,18 @@ PSIPrepare<T, Device>::PSIPrepare(const std::string& init_wfc_in,
 }
 
 template <typename T, typename Device>
-void PSIPrepare<T, Device>::prepare_init(const int& random_seed)
+void PSIPrepare<T, Device>::prepare_init(const int& random_seed,
+                                          const int& nbands,
+                                          const int& nspin,
+                                          const int& nqx,
+                                          const double& dq,
+                                          const bool& domag,
+                                          const bool& domag_z,
+                                          const bool& pseudo_mesh,
+                                          const std::string& orbital_dir,
+                                          const std::string& global_readin_dir,
+                                          const bool& ks_run,
+                                          const int& npol_in)
 {
 
     // under restriction of C++11, std::unique_ptr can not be allocate via std::make_unique
@@ -44,12 +57,19 @@ void PSIPrepare<T, Device>::prepare_init(const int& random_seed)
     this->psi_initer.reset();
     if (this->init_wfc == "random")
     {
-        this->psi_initer = std::unique_ptr<psi_initializer<T>>(new psi_init_random<T>());
-        GlobalV::ofs_running << "\n Using RANDOM starting wave functions for all " << PARAM.inp.nbands << " bands\n";
+        this->psi_initer = std::unique_ptr<psi_base<T>>(new psi_init_random<T>());
+        GlobalV::ofs_running << "\n Using RANDOM starting wave functions for all " << nbands << " bands\n";
     }
     else if (this->init_wfc == "file")
     {
-        this->psi_initer = std::unique_ptr<psi_initializer<T>>(new psi_init_file<T>());
+        psi_init_file<T>* file_initer = new psi_init_file<T>();
+        file_initer->prepare_params(
+            nspin,
+            global_readin_dir,
+            GlobalV::RANK_IN_POOL,
+            GlobalV::NPROC_IN_POOL
+        );
+        this->psi_initer = std::unique_ptr<psi_base<T>>(file_initer);
         GlobalV::ofs_running << "\n Using FILE starting wave functions\n";
     }
     else if ((this->init_wfc.substr(0, 6) == "atomic") && (this->ucell.natomwfc == 0))
@@ -57,7 +77,7 @@ void PSIPrepare<T, Device>::prepare_init(const int& random_seed)
         std::cout << " WARNING: init_wfc = " + this->init_wfc +
             " requires atomic pseudo wavefunctions(PP_PSWFC),\n but none available."
             " Automatically switch to random initialization." << std::endl;
-        GlobalV::ofs_running << "\n Using RANDOM starting wave functions for all " << PARAM.inp.nbands << " bands\n";
+        GlobalV::ofs_running << "\n Using RANDOM starting wave functions for all " << nbands << " bands\n";
         GlobalV::ofs_running << "\n WARNING:\n init_wfc = " + this->init_wfc + " requires atomic pseudo wavefunctions(PP_PSWFC), but none available. \n"
             " Automatically switch to random initialization.\n"
             " Note: Random starting wavefunctions may slow down convergence.\n"
@@ -65,39 +85,75 @@ void PSIPrepare<T, Device>::prepare_init(const int& random_seed)
             "      1) A pseudopotential file that includes atomic wavefunctions (with PP_PSWFC), or\n"
             "      2) Numerical atomic orbitals with 'init_wfc = nao' or 'nao+random' if available.\n"
             << std::endl;
-        this->psi_initer = std::unique_ptr<psi_initializer<T>>(new psi_init_random<T>());
+        this->psi_initer = std::unique_ptr<psi_base<T>>(new psi_init_random<T>());
     }
     else if (this->init_wfc == "atomic"
-             || (this->init_wfc == "atomic+random" && this->ucell.natomwfc < PARAM.inp.nbands))
+             || (this->init_wfc == "atomic+random" && this->ucell.natomwfc < nbands))
     {
-        if (this->ucell.natomwfc < PARAM.inp.nbands)
+        if (this->ucell.natomwfc < nbands)
         {
-            int nrandom = PARAM.inp.nbands - this->ucell.natomwfc;
+            int nrandom = nbands - this->ucell.natomwfc;
             GlobalV::ofs_running << "\n Using ATOMIC starting wave functions with " << this->ucell.natomwfc << " atomic orbitals"
             << " + " << nrandom << " random orbitals"
-            << " (total " << PARAM.inp.nbands << " bands)\n";
+            << " (total " << nbands << " bands)\n";
         }
         else
         {
             GlobalV::ofs_running << "\n Using ATOMIC starting wave functions for all " << this->ucell.natomwfc << " atomic orbitals"
-                << " (covers " << PARAM.inp.nbands << " bands)\n";
+                << " (covers " << nbands << " bands)\n";
         }
-        this->psi_initer = std::unique_ptr<psi_initializer<T>>(new psi_init_atomic<T>());
+        psi_init_atomic<T>* atomic_initer = new psi_init_atomic<T>();
+        atomic_initer->prepare_params(
+            nqx,
+            dq,
+            nspin,
+            domag,
+            domag_z,
+            pseudo_mesh,
+            this->sf
+        );
+        this->psi_initer = std::unique_ptr<psi_base<T>>(atomic_initer);
     }
     else if (this->init_wfc == "atomic+random")
     {
-        this->psi_initer = std::unique_ptr<psi_initializer<T>>(new psi_init_atomic_random<T>());
+        psi_init_atomic_random<T>* atomic_rand_initer = new psi_init_atomic_random<T>();
+        atomic_rand_initer->prepare_params(
+            nqx,
+            dq,
+            nspin,
+            domag,
+            domag_z,
+            pseudo_mesh,
+            this->sf
+        );
+        this->psi_initer = std::unique_ptr<psi_base<T>>(atomic_rand_initer);
         GlobalV::ofs_running << "\n Using ATOMIC+RANDOM starting wave functions with "
                              << this->ucell.natomwfc << " atomic orbitals\n";
     }
     else if (this->init_wfc == "nao")
     {
-        this->psi_initer = std::unique_ptr<psi_initializer<T>>(new psi_init_nao<T>());
+        psi_init_nao<T>* nao_initer = new psi_init_nao<T>();
+        nao_initer->prepare_params(
+            nqx,
+            dq,
+            nspin,
+            orbital_dir,
+            this->sf
+        );
+        this->psi_initer = std::unique_ptr<psi_base<T>>(nao_initer);
         GlobalV::ofs_running << "\n Using NAO starting wave functions\n";
     }
     else if (this->init_wfc == "nao+random")
     {
-        this->psi_initer = std::unique_ptr<psi_initializer<T>>(new psi_init_nao_random<T>());
+        psi_init_nao_random<T>* nao_rand_initer = new psi_init_nao_random<T>();
+        nao_rand_initer->prepare_params(
+            nqx,
+            dq,
+            nspin,
+            orbital_dir,
+            this->sf
+        );
+        this->psi_initer = std::unique_ptr<psi_base<T>>(nao_rand_initer);
         GlobalV::ofs_running << "\n Using NAO+RANDOM starting wave functions\n";
     }
     else
@@ -105,7 +161,8 @@ void PSIPrepare<T, Device>::prepare_init(const int& random_seed)
         ModuleBase::WARNING_QUIT("PSIInit::prepare_init", "for new psi initializer, init_wfc type not supported");
     }
 
-    this->psi_initer->initialize(&sf, &pw_wfc, &ucell, &kv, random_seed, &nlpp, rank);
+    this->psi_initer->initialize(&pw_wfc, &ucell, ik2iktot_, nkstot_, random_seed, lmaxkb, rank,
+                                     npol_in, nbands);
     this->psi_initer->tabulate();
 
     ModuleBase::timer::end("PSIPrepare", "prepare_init");
@@ -115,9 +172,18 @@ template <typename T, typename Device>
 void PSIPrepare<T, Device>::initialize_psi(Psi<std::complex<double>>* psi,
                                         psi::Psi<T, Device>* kspw_psi,
                                         hamilt::Hamilt<T, Device>* p_hamilt,
-                                        std::ofstream& ofs_running)
+                                        std::ofstream& ofs_running,
+                                        const std::string& device,
+                                        const std::string& precision,
+                                        const std::string& ks_solver_in,
+                                        const int& bndpar,
+                                        const bool& use_k_continuity,
+                                        const std::string& calculation,
+                                        const int& mem_saver,
+                                        const int& npol_in,
+                                        const bool& ks_run_in)
 {
-    if (kspw_psi->get_nbands() == 0 || (!PARAM.globalv.ks_run))
+    if (kspw_psi->get_nbands() == 0 || (!ks_run_in))
     {
         return;
     }
@@ -135,18 +201,18 @@ void PSIPrepare<T, Device>::initialize_psi(Psi<std::complex<double>>* psi,
     Psi<T>* psi_cpu = reinterpret_cast<psi::Psi<T>*>(psi);
     Psi<T, Device>* psi_device = kspw_psi;
 
-    bool fill = PARAM.inp.ks_solver != "bpcg" || GlobalV::MY_BNDGROUP == 0;
+    bool fill = ks_solver_in != "bpcg" || GlobalV::MY_BNDGROUP == 0;
     if (fill)
     {
         if (not_equal)
         {
             psi_cpu = new Psi<T>(1, nbands_start, nbasis, nbasis, true);
-            psi_device = PARAM.inp.device == "gpu" ? new psi::Psi<T, Device>(psi_cpu[0])
+            psi_device = device == "gpu" ? new psi::Psi<T, Device>(psi_cpu[0])
                                                    : reinterpret_cast<psi::Psi<T, Device>*>(psi_cpu);
         }
-        else if (PARAM.inp.precision == "single")
+        else if (precision == "single")
         {
-            if (PARAM.inp.device == "cpu")
+            if (device == "cpu")
             {
                 psi_cpu = reinterpret_cast<psi::Psi<T>*>(kspw_psi);
                 psi_device = kspw_psi;
@@ -163,7 +229,7 @@ void PSIPrepare<T, Device>::initialize_psi(Psi<std::complex<double>>* psi,
     // like (1, nbands, npwx), in which npwx is the maximal npw of all kpoints
     for (int ik = 0; ik < this->pw_wfc.nks; ik++)
     {
-        if(PARAM.inp.use_k_continuity && ik > 0) continue;
+        if(use_k_continuity && ik > 0) continue;
         //! Fix the wavefunction to initialize at given kpoint
         psi->fix_k(ik);
         kspw_psi->fix_k(ik);
@@ -213,21 +279,21 @@ void PSIPrepare<T, Device>::initialize_psi(Psi<std::complex<double>>* psi,
             }
         }
 #ifdef __MPI
-        if (PARAM.inp.ks_solver == "bpcg" && PARAM.inp.bndpar > 1)
+        if (ks_solver_in == "bpcg" && bndpar > 1)
         {
-            std::vector<int> sendcounts(PARAM.inp.bndpar);
-            std::vector<int> displs(PARAM.inp.bndpar);
+            std::vector<int> sendcounts(bndpar);
+            std::vector<int> displs(bndpar);
             MPI_Allgather(&nbands_l, 1, MPI_INT, sendcounts.data(), 1, MPI_INT, BP_WORLD);
             displs[0] = 0;
             sendcounts[0] *= nbasis;
-            for (int i = 1; i < PARAM.inp.bndpar; i++)
+            for (int i = 1; i < bndpar; i++)
             {
                 sendcounts[i] *= nbasis;
                 displs[i] = displs[i - 1] + sendcounts[i - 1];
             }
             if (GlobalV::MY_BNDGROUP == 0)
             {
-                for (int ip = 1; ip < PARAM.inp.bndpar; ++ip)
+                for (int ip = 1; ip < bndpar; ++ip)
                 {
                     Parallel_Common::send_data(psi_cpu->get_pointer() + displs[ip], sendcounts[ip], ip, 0, BP_WORLD);
                 }
@@ -246,12 +312,12 @@ void PSIPrepare<T, Device>::initialize_psi(Psi<std::complex<double>>* psi,
         if (not_equal)
         {
             delete psi_cpu;
-            if (PARAM.inp.device == "gpu")
+            if (device == "gpu")
             {
                 delete psi_device;
             }
         }
-        else if (PARAM.inp.precision == "single" && PARAM.inp.device == "gpu")
+        else if (precision == "single" && device == "gpu")
         {
             delete psi_cpu;
         }
@@ -276,7 +342,10 @@ void allocate_psi(Psi<std::complex<double>>*& psi,
                   const int& nks,
                   const std::vector<int>& ngk,
                   const int& nbands,
-                  const int& npwx)
+                  const int& npwx,
+                  const std::string& calculation,
+                  const int& mem_saver,
+                  const int& npol)
 {
     assert(npwx > 0);
     assert(nks > 0);
@@ -284,12 +353,12 @@ void allocate_psi(Psi<std::complex<double>>*& psi,
 
     delete psi;
     int nks2 = nks;
-    if (PARAM.inp.calculation == "nscf" && PARAM.inp.mem_saver == 1)
+    if (calculation == "nscf" && mem_saver == 1)
     {
         nks2 = 1;
     }
-    psi = new psi::Psi<std::complex<double>>(nks2, nbands, npwx * PARAM.globalv.npol, ngk, true);
-    const size_t memory_cost = sizeof(std::complex<double>) * nks2 * nbands * (PARAM.globalv.npol * npwx);
+    psi = new psi::Psi<std::complex<double>>(nks2, nbands, npwx * npol, ngk, true);
+    const size_t memory_cost = sizeof(std::complex<double>) * nks2 * nbands * (npol * npwx);
     std::cout << " MEMORY FOR PSI (MB)  : " << static_cast<double>(memory_cost) / 1024.0 / 1024.0 << std::endl;
     ModuleBase::Memory::record("Psi_PW", memory_cost);
 }
diff --git a/source/source_psi/psi_prepare.h b/source/source_psi/psi_prepare.h
index 4b35f545214..7136e15171e 100644
--- a/source/source_psi/psi_prepare.h
+++ b/source/source_psi/psi_prepare.h
@@ -1,13 +1,17 @@
 #ifndef PSI_PREPARE_H
 #define PSI_PREPARE_H
+#include <string>
 #include "source_hamilt/hamilt.h"
-#include "source_psi/psi_initializer.h"
+#include "source_psi/psi_base.h"
 #include "source_psi/psi_prepare_base.h"
 
+class UnitCell;
+class Structure_Factor;
+namespace ModulePW { class PW_Basis_K; }
+
 namespace psi
 {
 
-// This class is used to prepare the wavefunction
 template <typename T, typename Device = base_device::DEVICE_CPU>
 class PSIPrepare : public PSIPrepareBase
 {
@@ -18,76 +22,72 @@ class PSIPrepare : public PSIPrepareBase
             const int& rank,
             const UnitCell& ucell,
             const Structure_Factor& sf,
-            const K_Vectors& kv_in,
-            const pseudopot_cell_vnl& nlpp,
+            const std::vector<int>& ik2iktot,
+            const int& nkstot,
+            const int& lmaxkb,
             const ModulePW::PW_Basis_K& pw_wfc);
+
     ~PSIPrepare(){};
 
-    ///@brief prepare the wavefunction initialization
-    void prepare_init(const int& random_seed);
-
-    //------------------------ only for psi_initializer --------------------
-    /**
-     * @brief initialize the wavefunction
-     *
-     * @param psi store the wavefunction
-     * @param p_hamilt Hamiltonian operator
-     * @param ofs_running output stream for running information
-     * @param is_already_initpsi whether psi has been initialized
-     */
+    void prepare_init(const int& random_seed,
+                       const int& nbands,
+                       const int& nspin,
+                       const int& nqx,
+                       const double& dq,
+                       const bool& domag,
+                       const bool& domag_z,
+                       const bool& pseudo_mesh,
+                       const std::string& orbital_dir,
+                       const std::string& global_readin_dir,
+                       const bool& ks_run,
+                       const int& npol);
+
     void initialize_psi(Psi<std::complex<double>>* psi,
                         psi::Psi<T, Device>* kspw_psi,
                         hamilt::Hamilt<T, Device>* p_hamilt,
-                        std::ofstream& ofs_running);
+                        std::ofstream& ofs_running,
+                        const std::string& device,
+                        const std::string& precision,
+                        const std::string& ks_solver_in,
+                        const int& bndpar,
+                        const bool& use_k_continuity,
+                        const std::string& calculation,
+                        const int& mem_saver,
+                        const int& npol_in,
+                        const bool& ks_run_in);
 
-    /**
-     * @brief initialize NAOs in plane wave basis, only for LCAO_IN_PW
-     *
-     */
     void initialize_lcao_in_pw(Psi<T>* psi_local, std::ofstream& ofs_running);
 
-    // psi_initializer<T, Device>* psi_initer = nullptr;
-    // change to use smart pointer to manage the memory, and avoid memory leak
-    // while the std::make_unique() is not supported till C++14,
-    // so use the new and std::unique_ptr to manage the memory, but this makes new-delete not symmetric
-    std::unique_ptr<psi_initializer<T>> psi_initer;
+    std::unique_ptr<psi_base<T>> psi_initer;
 
   private:
-    // wavefunction initialization type
+
     std::string init_wfc = "none";
 
-    // Kohn-Sham solver type
     std::string ks_solver = "none";
 
-    // basis type
     std::string basis_type = "none";
 
-    // pw basis
     const ModulePW::PW_Basis_K& pw_wfc;
 
-    // parallel kpoints
-    const K_Vectors& kv;
+    const std::vector<int>& ik2iktot_;
+    const int nkstot_;
 
-    // unit cell
     const UnitCell& ucell;
 
-    // structure factor
     const Structure_Factor& sf;
 
-    // nonlocal pseudopotential
-    const pseudopot_cell_vnl& nlpp;
+    const int lmaxkb;
 
-    Device* ctx = {};                      ///< device
-    base_device::DEVICE_CPU* cpu_ctx = {}; ///< CPU device
-    const int rank;                        ///< MPI rank
+    Device* ctx = {};
+    base_device::DEVICE_CPU* cpu_ctx = {};
+    const int rank;
 
-    //-------------------------OP--------------------------------------------
     using syncmem_complex_op = base_device::memory::synchronize_memory_op<T, Device, Device>;
     using syncmem_h2d_op = base_device::memory::synchronize_memory_op<T, Device, base_device::DEVICE_CPU>;
 };
 
-///@brief allocate the wavefunction
-void allocate_psi(Psi<std::complex<double>>*& psi, const int& nks, const std::vector<int>& ngk, const int& nbands, const int& npwx);
+void allocate_psi(Psi<std::complex<double>>*& psi, const int& nks, const std::vector<int>& ngk, const int& nbands, const int& npwx, const std::string& calculation, const int& mem_saver, const int& npol);
 
 } // namespace psi
-#endif
\ No newline at end of file
+#endif
diff --git a/source/source_psi/psi_prepare_base.h b/source/source_psi/psi_prepare_base.h
index c7a10718bf9..d7f5f953fbb 100644
--- a/source/source_psi/psi_prepare_base.h
+++ b/source/source_psi/psi_prepare_base.h
@@ -1,22 +1,28 @@
 #ifndef PSI_PREPARE_BASE_H
 #define PSI_PREPARE_BASE_H
 
+#include <string>
+
 namespace psi
 {
 
-/**
- * @brief Base class for PSIPrepare without template parameters.
- *
- * This class provides a non-template base class for PSIPrepare<T, Device>,
- * allowing Setup_Psi_pw to store a base class pointer instead of a template pointer.
- * This is part of the gradual refactoring to remove template parameters from Setup_Psi_pw.
- */
 class PSIPrepareBase
 {
   public:
     PSIPrepareBase() = default;
     virtual ~PSIPrepareBase() = default;
-    virtual void prepare_init(const int& random_seed) = 0;
+    virtual void prepare_init(const int& random_seed,
+                              const int& nbands,
+                              const int& nspin,
+                              const int& nqx,
+                              const double& dq,
+                              const bool& domag,
+                              const bool& domag_z,
+                              const bool& pseudo_mesh,
+                              const std::string& orbital_dir,
+                              const std::string& global_readin_dir,
+                              const bool& ks_run,
+                              const int& npol) = 0;
 };
 
 } // namespace psi
diff --git a/source/source_psi/setup_psi.cpp b/source/source_psi/setup_psi.cpp
index ba658a02de4..cdb7a0c247b 100644
--- a/source/source_psi/setup_psi.cpp
+++ b/source/source_psi/setup_psi.cpp
@@ -1,4 +1,5 @@
 #include "source_psi/setup_psi.h"
+#include "source_cell/klist.h"
 #include "source_io/module_parameter/parameter.h" // use parameter
 
 template <typename T>
@@ -12,10 +13,10 @@ Setup_Psi<T>::~Setup_Psi(){}
 // In that case, psi may change its size multiple times during SCF
 template <typename T>
 void Setup_Psi<T>::allocate_psi(
-		psi::Psi<T>* &psi,
-		const K_Vectors &kv,
-        const Parallel_Orbitals &para_orb,
-		const Input_para &inp)
+    psi::Psi<T>* &psi,
+    const K_Vectors &kv,
+    const Parallel_Orbitals &para_orb,
+    const Input_para &inp)
 {
     // init electronic wave function psi
     if (psi == nullptr)
@@ -50,10 +51,10 @@ void Setup_Psi<T>::allocate_psi(
 template <typename T>
 void Setup_Psi<T>::deallocate_psi(psi::Psi<T>* &psi)
 {
-	if(psi!=nullptr)
-	{
-		delete psi;
-	}
+    if(psi!=nullptr)
+    {
+        delete psi;
+    }
 }
 
 template class Setup_Psi<double>;
diff --git a/source/source_psi/setup_psi.h b/source/source_psi/setup_psi.h
index a4c0f11d3e3..5ce3f1f2fe8 100644
--- a/source/source_psi/setup_psi.h
+++ b/source/source_psi/setup_psi.h
@@ -14,11 +14,11 @@ class Setup_Psi
     Setup_Psi();
     ~Setup_Psi();
 
-	static void allocate_psi(
-		psi::Psi<T>* &psi,
-		const K_Vectors &kv,
+    static void allocate_psi(
+        psi::Psi<T>* &psi,
+        const K_Vectors &kv,
         const Parallel_Orbitals &para_orb,
-		const Input_para &inp);
+        const Input_para &inp);
 
     static void deallocate_psi(psi::Psi<T>* &psi);
 
diff --git a/source/source_psi/setup_psi_pw.cpp b/source/source_psi/setup_psi_pw.cpp
index f5bc240292d..eab1adb361a 100644
--- a/source/source_psi/setup_psi_pw.cpp
+++ b/source/source_psi/setup_psi_pw.cpp
@@ -1,5 +1,9 @@
 #include "source_psi/setup_psi_pw.h"
-#include "source_io/module_parameter/parameter.h" // use parameter
+#include "source_cell/klist.h"
+#include "source_cell/unitcell.h"
+#include "source_pw/module_pwdft/structure_factor.h"
+#include "source_basis/module_pw/pw_basis_k.h"
+#include "source_io/module_parameter/parameter.h"
 
 Setup_Psi_pw::Setup_Psi_pw(){}
 
@@ -11,37 +15,62 @@ void Setup_Psi_pw::before_runner_impl(
         const K_Vectors &kv,
         const Structure_Factor &sf,
         const ModulePW::PW_Basis_K &pw_wfc, 
-        const pseudopot_cell_vnl &ppcell,
+        const int &lmaxkb,
         const Input_para &inp)
 {
     this->p_psi_init = new psi::PSIPrepare<T, Device>(inp.init_wfc,
       inp.ks_solver, inp.basis_type, GlobalV::MY_RANK, ucell,
-      sf, kv, ppcell, pw_wfc);
+      sf, kv.ik2iktot, kv.get_nkstot(), lmaxkb, pw_wfc);
 
-    allocate_psi(this->psi_cpu, kv.get_nks(), kv.ngk, PARAM.globalv.nbands_l, pw_wfc.npwk_max);
+    allocate_psi(this->psi_cpu, kv.get_nks(), kv.ngk, PARAM.globalv.nbands_l, pw_wfc.npwk_max,
+                  inp.calculation, inp.mem_saver, PARAM.globalv.npol);
 
     auto* p_psi_init = static_cast<psi::PSIPrepare<T, Device>*>(this->p_psi_init);
-    p_psi_init->prepare_init(inp.pw_seed);
-
-    if (std::is_same<T, float>::value) {
+    p_psi_init->prepare_init(inp.pw_seed,
+                             PARAM.globalv.nbands_l,
+                             inp.nspin,
+                             PARAM.globalv.nqx,
+                             PARAM.globalv.dq,
+                             PARAM.globalv.domag,
+                             PARAM.globalv.domag_z,
+                             inp.pseudo_mesh,
+                             inp.orbital_dir,
+                             PARAM.globalv.global_readin_dir,
+                             PARAM.globalv.ks_run,
+                             PARAM.globalv.npol);
+
+    if (std::is_same<T, float>::value) 
+    {
         precision_type_ = PrecisionType::Float;
-    } else if (std::is_same<T, double>::value) {
+    } 
+    else if (std::is_same<T, double>::value) 
+    {
         precision_type_ = PrecisionType::Double;
-    } else if (std::is_same<T, std::complex<float>>::value) {
+    } 
+    else if (std::is_same<T, std::complex<float>>::value) 
+    {
         precision_type_ = PrecisionType::ComplexFloat;
-    } else {
+    } 
+    else 
+    {
         precision_type_ = PrecisionType::ComplexDouble;
     }
     
-    if (std::is_same<Device, base_device::DEVICE_GPU>::value) {
+    if (std::is_same<Device, base_device::DEVICE_GPU>::value) 
+    {
         device_type_ = base_device::GpuDevice;
-    } else {
+    } 
+    else 
+    {
         device_type_ = base_device::CpuDevice;
     }
 
-    if (inp.device == "gpu" || inp.precision == "single") {
+    if (inp.device == "gpu" || inp.precision == "single") 
+    {
         this->psi_t = static_cast<void*>(new psi::Psi<T, Device>(this->psi_cpu[0]));
-    } else {
+    } 
+    else 
+    {
         this->psi_t = static_cast<void*>(reinterpret_cast<psi::Psi<T, Device>*>(this->psi_cpu));
     }
 }
@@ -51,30 +80,37 @@ void Setup_Psi_pw::before_runner(
         const K_Vectors &kv,
         const Structure_Factor &sf,
         const ModulePW::PW_Basis_K &pw_wfc, 
-        const pseudopot_cell_vnl &ppcell,
+        const int &lmaxkb,
         const Input_para &inp)
 {
     const bool is_gpu = (inp.device == "gpu");
     const bool is_single = (inp.precision == "single");
 
 #if ((defined __CUDA) || (defined __ROCM))
-    if (is_gpu) {
-        if (is_single) {
+    if (is_gpu) 
+    {
+        if (is_single) 
+	{
             before_runner_impl<std::complex<float>, base_device::DEVICE_GPU>(
-                ucell, kv, sf, pw_wfc, ppcell, inp);
-        } else {
+                ucell, kv, sf, pw_wfc, lmaxkb, inp);
+        } 
+	else 
+	{
             before_runner_impl<std::complex<double>, base_device::DEVICE_GPU>(
-                ucell, kv, sf, pw_wfc, ppcell, inp);
+                ucell, kv, sf, pw_wfc, lmaxkb, inp);
         }
     } else
 #endif
     {
-        if (is_single) {
+        if (is_single) 
+	{
             before_runner_impl<std::complex<float>, base_device::DEVICE_CPU>(
-                ucell, kv, sf, pw_wfc, ppcell, inp);
-        } else {
+                ucell, kv, sf, pw_wfc, lmaxkb, inp);
+        } 
+	else 
+	{
             before_runner_impl<std::complex<double>, base_device::DEVICE_CPU>(
-                ucell, kv, sf, pw_wfc, ppcell, inp);
+                ucell, kv, sf, pw_wfc, lmaxkb, inp);
         }
     }
 }
@@ -88,10 +124,12 @@ void Setup_Psi_pw::update_psi_d_impl()
         delete this->get_psi_d<T, Device>();
     }
 
-    // Refresh this->psi_d
-    if (this->precision_type_ == PrecisionType::ComplexFloat) {
+    if (this->precision_type_ == PrecisionType::ComplexFloat) 
+    {
         this->psi_d = static_cast<void*>(new psi::Psi<std::complex<double>, Device>(*this->get_psi_t<T, Device>()));
-    } else {
+    } 
+    else 
+    {
         this->psi_d = static_cast<void*>(reinterpret_cast<psi::Psi<std::complex<double>, Device>*>(this->psi_t));
     }
 }
@@ -125,17 +163,37 @@ void Setup_Psi_pw::update_psi_d()
 }
 
 template <typename T, typename Device>
-void Setup_Psi_pw::init_impl(hamilt::Hamilt<T, Device>* p_hamilt)
+void Setup_Psi_pw::init_impl(hamilt::Hamilt<T, Device>* p_hamilt,
+                             const std::string& device,
+                             const std::string& precision,
+                             const std::string& ks_solver,
+                             const int& bndpar,
+                             const bool& use_k_continuity,
+                             const std::string& calculation,
+                             const int& mem_saver,
+                             const int& npol,
+                             const bool& ks_run)
 {
     if (!this->already_initpsi)
     {
         auto* p_psi_init = static_cast<psi::PSIPrepare<T, Device>*>(this->p_psi_init);
-        p_psi_init->initialize_psi(this->psi_cpu, this->get_psi_t<T, Device>(), p_hamilt, GlobalV::ofs_running);
+        p_psi_init->initialize_psi(this->psi_cpu, this->get_psi_t<T, Device>(), p_hamilt, GlobalV::ofs_running,
+                                    device, precision, ks_solver, bndpar, use_k_continuity,
+                                    calculation, mem_saver, npol, ks_run);
         this->already_initpsi = true;
     }
 }
 
-void Setup_Psi_pw::init(hamilt::HamiltBase* p_hamilt)
+void Setup_Psi_pw::init(hamilt::HamiltBase* p_hamilt,
+                         const std::string& device,
+                         const std::string& precision,
+                         const std::string& ks_solver,
+                         const int& bndpar,
+                         const bool& use_k_continuity,
+                         const std::string& calculation,
+                         const int& mem_saver,
+                         const int& npol,
+                         const bool& ks_run)
 {
     if (this->already_initpsi)
     {
@@ -148,12 +206,16 @@ void Setup_Psi_pw::init(hamilt::HamiltBase* p_hamilt)
         if (this->precision_type_ == PrecisionType::ComplexFloat)
         {
             init_impl<std::complex<float>, base_device::DEVICE_GPU>(
-                static_cast<hamilt::Hamilt<std::complex<float>, base_device::DEVICE_GPU>*>(p_hamilt));
+                static_cast<hamilt::Hamilt<std::complex<float>, base_device::DEVICE_GPU>*>(p_hamilt),
+                device, precision, ks_solver, bndpar, use_k_continuity,
+                calculation, mem_saver, npol, ks_run);
         }
         else
         {
             init_impl<std::complex<double>, base_device::DEVICE_GPU>(
-                static_cast<hamilt::Hamilt<std::complex<double>, base_device::DEVICE_GPU>*>(p_hamilt));
+                static_cast<hamilt::Hamilt<std::complex<double>, base_device::DEVICE_GPU>*>(p_hamilt),
+                device, precision, ks_solver, bndpar, use_k_continuity,
+                calculation, mem_saver, npol, ks_run);
         }
     }
     else
@@ -162,12 +224,16 @@ void Setup_Psi_pw::init(hamilt::HamiltBase* p_hamilt)
         if (this->precision_type_ == PrecisionType::ComplexFloat)
         {
             init_impl<std::complex<float>, base_device::DEVICE_CPU>(
-                static_cast<hamilt::Hamilt<std::complex<float>, base_device::DEVICE_CPU>*>(p_hamilt));
+                static_cast<hamilt::Hamilt<std::complex<float>, base_device::DEVICE_CPU>*>(p_hamilt),
+                device, precision, ks_solver, bndpar, use_k_continuity,
+                calculation, mem_saver, npol, ks_run);
         }
         else
         {
             init_impl<std::complex<double>, base_device::DEVICE_CPU>(
-                static_cast<hamilt::Hamilt<std::complex<double>, base_device::DEVICE_CPU>*>(p_hamilt));
+                static_cast<hamilt::Hamilt<std::complex<double>, base_device::DEVICE_CPU>*>(p_hamilt),
+                device, precision, ks_solver, bndpar, use_k_continuity,
+                calculation, mem_saver, npol, ks_run);
         }
     }
 }
@@ -203,13 +269,13 @@ void Setup_Psi_pw::copy_d2h()
 }
 
 template <typename T, typename Device>
-void Setup_Psi_pw::castmem_d2h_impl(std::complex<double>* dst, const std::complex<double>* src, const size_t size)
+void Setup_Psi_pw::castmem_d2h_impl(std::complex<double>* dst, const std::complex<double>* src, const std::size_t size)
 {
     base_device::memory::cast_memory_op<std::complex<double>, std::complex<double>, base_device::DEVICE_CPU, Device>()(dst, src, size);
 }
 
 template <typename T, typename Device>
-void Setup_Psi_pw::castmem_d2h_impl(std::complex<double>* dst, const std::complex<float>* src, const size_t size)
+void Setup_Psi_pw::castmem_d2h_impl(std::complex<double>* dst, const std::complex<float>* src, const std::size_t size)
 {
     base_device::memory::cast_memory_op<std::complex<double>, std::complex<float>, base_device::DEVICE_CPU, Device>()(dst, src, size);
 }
@@ -263,17 +329,21 @@ template class psi::PSIPrepare<std::complex<double>, base_device::DEVICE_CPU>;
 
 template void Setup_Psi_pw::before_runner_impl<std::complex<float>, base_device::DEVICE_CPU>(
     const UnitCell&, const K_Vectors&, const Structure_Factor&,
-    const ModulePW::PW_Basis_K&, const pseudopot_cell_vnl&, const Input_para&);
+    const ModulePW::PW_Basis_K&, const int&, const Input_para&);
 
 template void Setup_Psi_pw::before_runner_impl<std::complex<double>, base_device::DEVICE_CPU>(
     const UnitCell&, const K_Vectors&, const Structure_Factor&,
-    const ModulePW::PW_Basis_K&, const pseudopot_cell_vnl&, const Input_para&);
+    const ModulePW::PW_Basis_K&, const int&, const Input_para&);
 
 template void Setup_Psi_pw::init_impl<std::complex<float>, base_device::DEVICE_CPU>(
-    hamilt::Hamilt<std::complex<float>, base_device::DEVICE_CPU>*);
+    hamilt::Hamilt<std::complex<float>, base_device::DEVICE_CPU>*,
+    const std::string&, const std::string&, const std::string&,
+    const int&, const bool&, const std::string&, const int&, const int&, const bool&);
 
 template void Setup_Psi_pw::init_impl<std::complex<double>, base_device::DEVICE_CPU>(
-    hamilt::Hamilt<std::complex<double>, base_device::DEVICE_CPU>*);
+    hamilt::Hamilt<std::complex<double>, base_device::DEVICE_CPU>*,
+    const std::string&, const std::string&, const std::string&,
+    const int&, const bool&, const std::string&, const int&, const int&, const bool&);
 
 template void Setup_Psi_pw::update_psi_d_impl<std::complex<float>, base_device::DEVICE_CPU>();
 
@@ -284,16 +354,16 @@ template void Setup_Psi_pw::clean_impl<std::complex<float>, base_device::DEVICE_
 template void Setup_Psi_pw::clean_impl<std::complex<double>, base_device::DEVICE_CPU>();
 
 template void Setup_Psi_pw::castmem_d2h_impl<std::complex<float>, base_device::DEVICE_CPU>(
-    std::complex<double>*, const std::complex<float>*, const size_t);
+    std::complex<double>*, const std::complex<float>*, const std::size_t);
 
 template void Setup_Psi_pw::castmem_d2h_impl<std::complex<float>, base_device::DEVICE_CPU>(
-    std::complex<double>*, const std::complex<double>*, const size_t);
+    std::complex<double>*, const std::complex<double>*, const std::size_t);
 
 template void Setup_Psi_pw::castmem_d2h_impl<std::complex<double>, base_device::DEVICE_CPU>(
-    std::complex<double>*, const std::complex<float>*, const size_t);
+    std::complex<double>*, const std::complex<float>*, const std::size_t);
 
 template void Setup_Psi_pw::castmem_d2h_impl<std::complex<double>, base_device::DEVICE_CPU>(
-    std::complex<double>*, const std::complex<double>*, const size_t);
+    std::complex<double>*, const std::complex<double>*, const std::size_t);
 
 #if ((defined __CUDA) || (defined __ROCM))
 template class psi::PSIPrepare<std::complex<float>, base_device::DEVICE_GPU>;
@@ -301,17 +371,21 @@ template class psi::PSIPrepare<std::complex<double>, base_device::DEVICE_GPU>;
 
 template void Setup_Psi_pw::before_runner_impl<std::complex<float>, base_device::DEVICE_GPU>(
     const UnitCell&, const K_Vectors&, const Structure_Factor&,
-    const ModulePW::PW_Basis_K&, const pseudopot_cell_vnl&, const Input_para&);
+    const ModulePW::PW_Basis_K&, const int&, const Input_para&);
 
 template void Setup_Psi_pw::before_runner_impl<std::complex<double>, base_device::DEVICE_GPU>(
     const UnitCell&, const K_Vectors&, const Structure_Factor&,
-    const ModulePW::PW_Basis_K&, const pseudopot_cell_vnl&, const Input_para&);
+    const ModulePW::PW_Basis_K&, const int&, const Input_para&);
 
 template void Setup_Psi_pw::init_impl<std::complex<float>, base_device::DEVICE_GPU>(
-    hamilt::Hamilt<std::complex<float>, base_device::DEVICE_GPU>*);
+    hamilt::Hamilt<std::complex<float>, base_device::DEVICE_GPU>*,
+    const std::string&, const std::string&, const std::string&,
+    const int&, const bool&, const std::string&, const int&, const int&, const bool&);
 
 template void Setup_Psi_pw::init_impl<std::complex<double>, base_device::DEVICE_GPU>(
-    hamilt::Hamilt<std::complex<double>, base_device::DEVICE_GPU>*);
+    hamilt::Hamilt<std::complex<double>, base_device::DEVICE_GPU>*,
+    const std::string&, const std::string&, const std::string&,
+    const int&, const bool&, const std::string&, const int&, const int&, const bool&);
 
 template void Setup_Psi_pw::update_psi_d_impl<std::complex<float>, base_device::DEVICE_GPU>();
 
@@ -326,14 +400,14 @@ template void Setup_Psi_pw::clean_impl<std::complex<float>, base_device::DEVICE_
 template void Setup_Psi_pw::clean_impl<std::complex<double>, base_device::DEVICE_GPU>();
 
 template void Setup_Psi_pw::castmem_d2h_impl<std::complex<float>, base_device::DEVICE_GPU>(
-    std::complex<double>*, const std::complex<float>*, const size_t);
+    std::complex<double>*, const std::complex<float>*, const std::size_t);
 
 template void Setup_Psi_pw::castmem_d2h_impl<std::complex<float>, base_device::DEVICE_GPU>(
-    std::complex<double>*, const std::complex<double>*, const size_t);
+    std::complex<double>*, const std::complex<double>*, const std::size_t);
 
 template void Setup_Psi_pw::castmem_d2h_impl<std::complex<double>, base_device::DEVICE_GPU>(
-    std::complex<double>*, const std::complex<float>*, const size_t);
+    std::complex<double>*, const std::complex<float>*, const std::size_t);
 
 template void Setup_Psi_pw::castmem_d2h_impl<std::complex<double>, base_device::DEVICE_GPU>(
-    std::complex<double>*, const std::complex<double>*, const size_t);
+    std::complex<double>*, const std::complex<double>*, const std::size_t);
 #endif
diff --git a/source/source_psi/setup_psi_pw.h b/source/source_psi/setup_psi_pw.h
index 88e9d42bf1c..f1c338a4338 100644
--- a/source/source_psi/setup_psi_pw.h
+++ b/source/source_psi/setup_psi_pw.h
@@ -6,10 +6,10 @@
 #include "source_cell/klist.h"
 #include "source_pw/module_pwdft/structure_factor.h"
 #include "source_basis/module_pw/pw_basis_k.h"
-#include "source_pw/module_pwdft/vnl_pw.h"
 #include "source_io/module_parameter/input_parameter.h"
 #include "source_base/module_device/device.h"
 #include "source_hamilt/hamilt.h"
+#include "source_psi/psi_base.h"
 
 class Setup_Psi_pw
 {
@@ -39,7 +39,7 @@ class Setup_Psi_pw
     // for PW, we have psi_cpu
     psi::Psi<std::complex<double>, base_device::DEVICE_CPU>* psi_cpu = nullptr;
 
-    // psi_initializer controller
+    // psi_base controller
     psi::PSIPrepareBase* p_psi_init = nullptr;
 
     //------------
@@ -51,10 +51,19 @@ class Setup_Psi_pw
         const K_Vectors &kv,
         const Structure_Factor &sf,
         const ModulePW::PW_Basis_K &pw_wfc, 
-        const pseudopot_cell_vnl &ppcell,
+        const int &lmaxkb,
         const Input_para &inp);
 
-    void init(hamilt::HamiltBase* p_hamilt);
+    void init(hamilt::HamiltBase* p_hamilt,
+              const std::string& device,
+              const std::string& precision,
+              const std::string& ks_solver,
+              const int& bndpar,
+              const bool& use_k_continuity,
+              const std::string& calculation,
+              const int& mem_saver,
+              const int& npol,
+              const bool& ks_run);
 
     void update_psi_d();
 
@@ -71,7 +80,7 @@ class Setup_Psi_pw
     int get_nbands() const { return this->psi_cpu->get_nbands(); }
     int get_nk() const { return this->psi_cpu->get_nk(); }
     int get_nbasis() const { return this->psi_cpu->get_nbasis(); }
-    size_t size() const { return this->psi_cpu->size(); }
+    std::size_t size() const { return this->psi_cpu->size(); }
     
     // Get runtime type information
     base_device::AbacusDevice_t get_device_type() const { return device_type_; }
@@ -126,11 +135,20 @@ class Setup_Psi_pw
         const K_Vectors &kv,
         const Structure_Factor &sf,
         const ModulePW::PW_Basis_K &pw_wfc, 
-        const pseudopot_cell_vnl &ppcell,
+        const int &lmaxkb,
         const Input_para &inp);
 
     template <typename T, typename Device>
-    void init_impl(hamilt::Hamilt<T, Device>* p_hamilt);
+    void init_impl(hamilt::Hamilt<T, Device>* p_hamilt,
+                   const std::string& device,
+                   const std::string& precision,
+                   const std::string& ks_solver,
+                   const int& bndpar,
+                   const bool& use_k_continuity,
+                   const std::string& calculation,
+                   const int& mem_saver,
+                   const int& npol,
+                   const bool& ks_run);
 
     template <typename T, typename Device>
     void update_psi_d_impl();
@@ -142,10 +160,10 @@ class Setup_Psi_pw
     void copy_d2h_impl();
 
     template <typename T, typename Device>
-    void castmem_d2h_impl(std::complex<double>* dst, const std::complex<double>* src, const size_t size);
+    void castmem_d2h_impl(std::complex<double>* dst, const std::complex<double>* src, const std::size_t size);
     
     template <typename T, typename Device>
-    void castmem_d2h_impl(std::complex<double>* dst, const std::complex<float>* src, const size_t size);
+    void castmem_d2h_impl(std::complex<double>* dst, const std::complex<float>* src, const std::size_t size);
 
 };
 
diff --git a/source/source_pw/module_pwdft/soc.cpp b/source/source_psi/soc_base.cpp
similarity index 93%
rename from source/source_pw/module_pwdft/soc.cpp
rename to source/source_psi/soc_base.cpp
index a96c7ad2351..92a1b15e066 100644
--- a/source/source_pw/module_pwdft/soc.cpp
+++ b/source/source_psi/soc_base.cpp
@@ -1,4 +1,4 @@
-#include "soc.h"
+#include "soc_base.h"
 
 Fcoef::~Fcoef()
 {
@@ -15,15 +15,21 @@ void Fcoef::create(const int i1, const int i2, const int i3)
         ind1 = i1;
         ind4 = i2;
         ind5 = i3;
+
         if (this->p != nullptr)
         {
             delete[] this->p;
             this->p = nullptr;
         }
+
         int tot = ind1 * ind2 * ind3 * ind4 * ind5;
+
         this->p = new std::complex<double>[tot];
+
         for (int i = 0; i < tot; i++)
+	{
             this->p[i] = std::complex<double>(0.0, 0.0);
+	}
     }
     else
     {
@@ -46,7 +52,7 @@ void Fcoef::create(const int i1, const int i2, const int i3)
     return;
 }
 
-Soc::~Soc()
+SocBase::~SocBase()
 {
     if (this->p_rot != nullptr)
     {
@@ -54,12 +60,16 @@ Soc::~Soc()
     }
 }
 
-double Soc::spinor(const int l, const double j, const int m, const int spin) const
+double SocBase::spinor(const int l, const double j, const int m, const int spin) const
 {
     if (spin != 0 && spin != 1)
+    {
         ModuleBase::WARNING_QUIT("spinor", "spin direction unknown");
+    }
     if (m < -l - 1 || m > l)
+    {
         ModuleBase::WARNING_QUIT("spinor", "m not allowed");
+    }
 
     double den = 1.0 / (2.0 * l + 1.0); // denominator
 
@@ -86,9 +96,13 @@ double Soc::spinor(const int l, const double j, const int m, const int spin) con
         else
         {
             if (spin == 0)
+	    {
                 spinor0 = sqrt((l - m + 1.0) * den);
+	    }
             if (spin == 1)
+	    {
                 spinor0 = -sqrt((l + m) * den);
+	    }
         }
     }
     else
@@ -99,7 +113,7 @@ double Soc::spinor(const int l, const double j, const int m, const int spin) con
     return spinor0;
 }
 
-void Soc::rot_ylm(const int lmax)
+void SocBase::rot_ylm(const int lmax)
 {
     // initialize the l_max_ and l2plus1_
     this->l_max_ = lmax;
@@ -129,7 +143,7 @@ void Soc::rot_ylm(const int lmax)
     return;
 }
 
-int Soc::sph_ind(const int l, const double j, const int m, const int spin) const
+int SocBase::sph_ind(const int l, const double j, const int m, const int spin) const
 {
     // This function calculates the m index of the spherical harmonic
     // in a spinor with orbital angular momentum l, total angular
@@ -147,9 +161,13 @@ int Soc::sph_ind(const int l, const double j, const int m, const int spin) const
     if (fabs(j - l - 0.5) < 1e-8)
     {
         if (spin == 0)
+	{
             sph_ind0 = m;
+	}
         if (spin == 1)
+	{
             sph_ind0 = m + 1;
+	}
     }
     else if (fabs(j - l + 0.5) < 1e-8)
     {
@@ -181,7 +199,7 @@ int Soc::sph_ind(const int l, const double j, const int m, const int spin) const
     return sph_ind0;
 }
 
-void Soc::set_fcoef(const int &l1,
+void SocBase::set_fcoef(const int &l1,
                     const int &l2,
                     const int &is1,
                     const int &is2,
@@ -201,4 +219,4 @@ void Soc::set_fcoef(const int &l1,
         coeff += rotylm(m1, mi) * spinor(l1, j1, m, is1) * conj(rotylm(m2, mj)) * spinor(l2, j2, m, is2);
     }
     this->fcoef(it, is1, is2, ip1, ip2) = coeff;
-}
\ No newline at end of file
+}
diff --git a/source/source_pw/module_pwdft/soc.h b/source/source_psi/soc_base.h
similarity index 89%
rename from source/source_pw/module_pwdft/soc.h
rename to source/source_psi/soc_base.h
index 40a3b2e41ac..a29930198dd 100644
--- a/source/source_pw/module_pwdft/soc.h
+++ b/source/source_psi/soc_base.h
@@ -1,7 +1,7 @@
-#ifndef SOC_H
-#define SOC_H
+#ifndef SOC_BASE_H
+#define SOC_BASE_H
 
-#include "source_base/global_function.h"
+#include "../source_base/global_function.h"
 
 #include <complex>
 #include <iostream>
@@ -42,22 +42,18 @@ class Fcoef
     int ind3 = 2;
 };
 
-//-----------------------
-// spin-orbital coupling
-//-----------------------
-class Soc
+class SocBase
 {
 
   public:
-    Soc(){};
-    ~Soc();
+    SocBase(){};
+    ~SocBase();
 
     double spinor(const int l, const double j, const int m, const int spin) const;
 
     int sph_ind(const int l, const double j, const int m, const int spin) const;
 
     void rot_ylm(const int lmax);
-    // std::complex<double> **rotylm;
 
     const std::complex<double> &rotylm(const int &i1, const int &i2) const
     {
@@ -65,6 +61,7 @@ class Soc
     }
 
     Fcoef fcoef;
+
     void set_fcoef(const int &l1,
                    const int &l2,
                    const int &is1,
@@ -77,10 +74,12 @@ class Soc
                    const int &ip1,
                    const int &ip2);
 
-    // int npol;
   private:
+
     int l_max_ = -1; // maximum orbital index, must initialized before set_fcoef
+
     int l2plus1_ = -1;
+
     std::complex<double> *p_rot = nullptr;
 };
 #endif
diff --git a/source/source_psi/test/CMakeLists.txt b/source/source_psi/test/CMakeLists.txt
index 63af5799e1a..87c89330382 100644
--- a/source/source_psi/test/CMakeLists.txt
+++ b/source/source_psi/test/CMakeLists.txt
@@ -8,13 +8,12 @@ AddTest(
 
 if(ENABLE_LCAO)
 AddTest(
-    TARGET MODULE_PSI_initializer_unit_test
-    LIBS parameter base device psi psi_initializer planewave
+    TARGET MODULE_PSI_init_test
+    LIBS parameter base device psi psi_init planewave
     SOURCES
-        psi_initializer_unit_test.cpp
-        ../../source_pw/module_pwdft/soc.cpp
+        psi_init_test.cpp
+        soc_base_test.cpp
         ../../source_cell/atom_spec.cpp
-        ../../source_cell/parallel_kpoints.cpp
         ../../source_cell/test/support/mock_unitcell.cpp
         ../../source_io/module_output/orb_io.cpp
         ../../source_io/module_output/write_pao.cpp
diff --git a/source/source_psi/test/psi_initializer_unit_test.cpp b/source/source_psi/test/psi_init_test.cpp
similarity index 52%
rename from source/source_psi/test/psi_initializer_unit_test.cpp
rename to source/source_psi/test/psi_init_test.cpp
index 342e04a5a60..b7c653b457e 100644
--- a/source/source_psi/test/psi_initializer_unit_test.cpp
+++ b/source/source_psi/test/psi_init_test.cpp
@@ -1,55 +1,61 @@
 #include <gtest/gtest.h>
-#define private public
-#include "source_io/module_parameter/parameter.h"
-#undef private
-#include "../psi_initializer.h"
+#include <complex>
+#include <functional>
+
+#include "source_pw/module_pwdft/vl_pw.h"
+#include "source_pw/module_pwdft/structure_factor.h"
+#include "source_pw/module_pwdft/parallel_grid.h"
+#include "source_basis/module_pw/pw_basis.h"
+#include "source_cell/pseudo.h"
+#include "source_cell/atom_pseudo.h"
+#include "source_cell/magnetism.h"
+#include "source_cell/unitcell.h"
+#include "../psi_base.h"
 #include "../psi_init_atomic.h"
 #include "../psi_init_atomic_random.h"
 #include "../psi_init_nao.h"
 #include "../psi_init_nao_random.h"
 #include "../psi_init_random.h"
-#include "source_pw/module_pwdft/vl_pw.h"
-#include "source_cell/klist.h"
 #include "source_base/output.h"
 
 /*
 =========================
-psi initializer unit test
+psi base unit test
 =========================
 - Tested functions:
-    - psi_initializer_random::psi_initializer_random
-      - constructor of psi_initializer_random
-    - psi_initializer_atomic::psi_initializer_atomic
-      - constructor of psi_initializer_atomic
-    - psi_initializer_atomic_random::psi_initializer_atomic_random
-      - constructor of psi_initializer_atomic_random
-    - psi_initializer_nao::psi_initializer_nao
-      - constructor of psi_initializer_nao
-    - psi_initializer_nao_random::psi_initializer_nao_random
-      - constructor of psi_initializer_nao_random
-    - psi_initializer::cast_to_T (psi_initializer specialized as random)
+    - psi_init_random::psi_init_random
+      - constructor of psi_init_random
+    - psi_init_atomic::psi_init_atomic
+      - constructor of psi_init_atomic
+    - psi_init_atomic_random::psi_init_atomic_random
+      - constructor of psi_init_atomic_random
+    - psi_init_nao::psi_init_nao
+      - constructor of psi_init_nao
+    - psi_init_nao_random::psi_init_nao_random
+      - constructor of psi_init_nao_random
+    - psi_base::cast_to_T (psi_base specialized as random)
       - function cast std::complex<double> to float, double, std::complex<float>, std::complex<double>
-    - psi_initializer_random::allocate
+    - psi_init_random::allocate
       - allocate wavefunctions with random-specific method
-    - psi_initializer_atomic::allocate
+    - psi_init_atomic::allocate
       - allocate wavefunctions with atomic-specific method
-    - psi_initializer_atomic_random::allocate
+    - psi_init_atomic_random::allocate
       - allocate wavefunctions with atomic-specific method
-    - psi_initializer_nao::allocate
+    - psi_init_nao::allocate
       - allocate wavefunctions with nao-specific method
-    - psi_initializer_nao_random::allocate
+    - psi_init_nao_random::allocate
       - allocate wavefunctions with nao-specific method
-    - psi_initializer_random::proj_ao_onkG
+    - psi_init_random::proj_ao_onkG
       - calculate wavefunction initial guess (before diagonalization) by randomly generating numbers
-    - psi_initializer_atomic::proj_ao_onkG
+    - psi_init_atomic::proj_ao_onkG
       - calculate wavefunction initial guess (before diagonalization) with atomic pseudo wavefunctions
       - nspin = 4 case
       - nspin = 4 with has_so case
-    - psi_initializer_atomic_random::proj_ao_onkG
+    - psi_init_atomic_random::proj_ao_onkG
       - calculate wavefunction initial guess (before diagonalization) with atomic pseudo wavefunctions and random numbers
-    - psi_initializer_nao::proj_ao_onkG
+    - psi_init_nao::proj_ao_onkG
       - calculate wavefunction initial guess (before diagonalization) with numerical atomic orbital wavefunctions
-    - psi_initializer_nao_random::proj_ao_onkG
+    - psi_init_nao_random::proj_ao_onkG
       - calculate wavefunction initial guess (before diagonalization) with numerical atomic orbital wavefunctions and random numbers
 */
 
@@ -64,10 +70,6 @@ void Atom_pseudo::bcast_atom_pseudo() {}
 pseudo::pseudo() {}
 pseudo::~pseudo() {}
 
-pseudopot_cell_vnl::pseudopot_cell_vnl() {}
-pseudopot_cell_vnl::~pseudopot_cell_vnl()
-{
-}
 pseudopot_cell_vl::pseudopot_cell_vl() {}
 pseudopot_cell_vl::~pseudopot_cell_vl() {}
 Magnetism::Magnetism() {}
@@ -86,8 +88,10 @@ std::complex<double>* Structure_Factor::get_sk(int ik, int it, int ia, ModulePW:
 {
     int npw = wfc_basis->npwk[ik];
     std::complex<double> *sk = new std::complex<double>[npw];
-    for(int ipw = 0; ipw < npw; ++ipw) { sk[ipw] = std::complex<double>(0.0, 0.0);
-}
+    for(int ipw = 0; ipw < npw; ++ipw)
+    {
+        sk[ipw] = std::complex<double>(0.0, 0.0);
+    }
     return sk;
 }
 
@@ -96,13 +100,23 @@ class PsiIntializerUnitTest : public ::testing::Test {
         Structure_Factor* p_sf = nullptr;
         ModulePW::PW_Basis_K* p_pw_wfc = nullptr;
         UnitCell* p_ucell = nullptr;
-        pseudopot_cell_vnl* p_pspot_vnl = nullptr;
-        K_Vectors* p_kv = nullptr;
+        int lmaxkb = 0;
+        std::vector<int> ik2iktot_;
+        int nkstot_ = 0;
         int random_seed = 1;
 
-        psi_initializer<std::complex<double>>* psi_init;
+        psi_base<std::complex<double>>* psi_init;
+
+        int nbands_ = 1;
+        int nspin_ = 1;
+        int npol_ = 1;
+        bool domag_ = false;
+        bool domag_z_ = false;
+        std::string orbital_dir_ = "./support/";
+        int nqx_ = 100;
+        double dq_ = 0.01;
+        bool pseudo_mesh_ = false;
 
-      private:
       protected:
         void SetUp() override
         {
@@ -110,19 +124,6 @@ class PsiIntializerUnitTest : public ::testing::Test {
             this->p_sf = new Structure_Factor();
             this->p_pw_wfc = new ModulePW::PW_Basis_K();
             this->p_ucell = new UnitCell();
-            this->p_pspot_vnl = new pseudopot_cell_vnl();
-            this->p_kv = new K_Vectors();
-            // mock
-            PARAM.input.nbands = 1;
-            PARAM.input.nspin = 1;
-            PARAM.input.orbital_dir = "./support/";
-            PARAM.input.pseudo_dir = "./support/";
-            PARAM.sys.npol = 1;
-            PARAM.input.calculation = "scf";
-            PARAM.input.init_wfc = "random";
-            PARAM.input.ks_solver = "cg";
-            PARAM.sys.domag = false;
-            PARAM.sys.domag_z = false;
             // lattice
             this->p_ucell->a1 = {10.0, 0.0, 0.0};
             this->p_ucell->a2 = {0.0, 10.0, 0.0};
@@ -165,19 +166,28 @@ class PsiIntializerUnitTest : public ::testing::Test {
             this->p_ucell->atoms[0].ncpp.mesh = 11;
             this->p_ucell->atoms[0].ncpp.msh = 11;
             this->p_ucell->atoms[0].ncpp.lmax = 2;
-            //if(this->p_ucell->atoms[0].ncpp.rab != nullptr) delete[] this->p_ucell->atoms[0].ncpp.rab;
+
             this->p_ucell->atoms[0].ncpp.rab = std::vector<double>(11, 0.0);
-            for(int i = 0; i < 11; ++i) { this->p_ucell->atoms[0].ncpp.rab[i] = 0.01;
-}
-            //if(this->p_ucell->atoms[0].ncpp.r != nullptr) delete[] this->p_ucell->atoms[0].ncpp.r;
+            for(int i = 0; i < 11; ++i)
+            {
+                this->p_ucell->atoms[0].ncpp.rab[i] = 0.01;
+            }
+
             this->p_ucell->atoms[0].ncpp.r = std::vector<double>(11, 0.0);
-            for(int i = 0; i < 11; ++i) { this->p_ucell->atoms[0].ncpp.r[i] = 0.01*i;
-}
+            for(int i = 0; i < 11; ++i)
+            {
+                this->p_ucell->atoms[0].ncpp.r[i] = 0.01*i;
+            }
+
             this->p_ucell->atoms[0].ncpp.chi.create(2, 11);
-            for(int i = 0; i < 2; ++i) { for(int j = 0; j < 11; ++j) { this->p_ucell->atoms[0].ncpp.chi(i, j) = 0.01;
-}
-}
-            //if(this->p_ucell->atoms[0].ncpp.lchi != nullptr) delete[] this->p_ucell->atoms[0].ncpp.lchi;
+            for(int i = 0; i < 2; ++i)
+            {
+                for(int j = 0; j < 11; ++j)
+                {
+                    this->p_ucell->atoms[0].ncpp.chi(i, j) = 0.01;
+                }
+            }
+
             this->p_ucell->atoms[0].ncpp.lchi = std::vector<int>(2, 0);
             this->p_ucell->atoms[0].ncpp.lchi[0] = 0;
             this->p_ucell->atoms[0].ncpp.lchi[1] = 1;
@@ -190,6 +200,7 @@ class PsiIntializerUnitTest : public ::testing::Test {
             this->p_ucell->atoms[0].ncpp.jchi = std::vector<double>(2, 0.0);
             this->p_ucell->atoms[0].ncpp.jchi[0] = 0.5;
             this->p_ucell->atoms[0].ncpp.jchi[1] = 1.5;
+
             // atom numerical orbital
             this->p_ucell->lmax = 2;
             p_ucell->orbital_fn.shrink_to_fit();
@@ -201,64 +212,93 @@ class PsiIntializerUnitTest : public ::testing::Test {
             this->p_ucell->atoms[0].l_nchi[1] = 2;
             this->p_ucell->atoms[0].l_nchi[2] = 1;
 
-            
+
             // can support function PW_Basis::getfftixy2is
             this->p_pw_wfc->nks = 1;
             this->p_pw_wfc->npwk_max = 1;
-            if(this->p_pw_wfc->npwk != nullptr) { delete[] this->p_pw_wfc->npwk;
-}
+            if(this->p_pw_wfc->npwk != nullptr)
+            {
+                delete[] this->p_pw_wfc->npwk;
+            }
+
             this->p_pw_wfc->npwk = new int[1];
             this->p_pw_wfc->npwk[0] = 1;
             this->p_pw_wfc->fftnxy = 1;
             this->p_pw_wfc->fftnz = 1;
             this->p_pw_wfc->nst = 1;
             this->p_pw_wfc->nz = 1;
-            if(this->p_pw_wfc->is2fftixy != nullptr) { delete[] this->p_pw_wfc->is2fftixy;
-}
+            if(this->p_pw_wfc->is2fftixy != nullptr)
+            {
+                delete[] this->p_pw_wfc->is2fftixy;
+            }
+
             this->p_pw_wfc->is2fftixy = new int[1];
             this->p_pw_wfc->is2fftixy[0] = 0;
-            if(this->p_pw_wfc->fftixy2ip != nullptr) { delete[] this->p_pw_wfc->fftixy2ip;
-}
+
+            if(this->p_pw_wfc->fftixy2ip != nullptr)
+            {
+                delete[] this->p_pw_wfc->fftixy2ip;
+            }
+
             this->p_pw_wfc->fftixy2ip = new int[1];
             this->p_pw_wfc->fftixy2ip[0] = 0;
-            if(this->p_pw_wfc->igl2isz_k != nullptr) { delete[] this->p_pw_wfc->igl2isz_k;
-}
+            if(this->p_pw_wfc->igl2isz_k != nullptr)
+            {
+                delete[] this->p_pw_wfc->igl2isz_k;
+            }
             this->p_pw_wfc->igl2isz_k = new int[1];
             this->p_pw_wfc->igl2isz_k[0] = 0;
-            if(this->p_pw_wfc->gcar != nullptr) { delete[] this->p_pw_wfc->gcar;
-}
+            if(this->p_pw_wfc->igl2ig_k != nullptr)
+            {
+                delete[] this->p_pw_wfc->igl2ig_k;
+            }
+            this->p_pw_wfc->igl2ig_k = new int[1];
+            this->p_pw_wfc->igl2ig_k[0] = 0;
+            if(this->p_pw_wfc->gcar != nullptr)
+            {
+                delete[] this->p_pw_wfc->gcar;
+            }
             this->p_pw_wfc->gcar = new ModuleBase::Vector3<double>[1];
             this->p_pw_wfc->gcar[0] = {0.0, 0.0, 0.0};
-            if(this->p_pw_wfc->igl2isz_k != nullptr) { delete[] this->p_pw_wfc->igl2isz_k;
-}
-            this->p_pw_wfc->igl2isz_k = new int[1];
-            this->p_pw_wfc->igl2isz_k[0] = 0;
-            if(this->p_pw_wfc->gk2 != nullptr) { delete[] this->p_pw_wfc->gk2;
-}
+            if(this->p_pw_wfc->gk2 != nullptr)
+            {
+                delete[] this->p_pw_wfc->gk2;
+            }
             this->p_pw_wfc->gk2 = new double[1];
             this->p_pw_wfc->gk2[0] = 0.0;
-            this->p_pw_wfc->latvec.e11 = this->p_ucell->latvec.e11; this->p_pw_wfc->latvec.e12 = this->p_ucell->latvec.e12; this->p_pw_wfc->latvec.e13 = this->p_ucell->latvec.e13;
-            this->p_pw_wfc->latvec.e21 = this->p_ucell->latvec.e21; this->p_pw_wfc->latvec.e22 = this->p_ucell->latvec.e22; this->p_pw_wfc->latvec.e23 = this->p_ucell->latvec.e23;
-            this->p_pw_wfc->latvec.e31 = this->p_ucell->latvec.e31; this->p_pw_wfc->latvec.e32 = this->p_ucell->latvec.e32; this->p_pw_wfc->latvec.e33 = this->p_ucell->latvec.e33;
+            this->p_pw_wfc->latvec.e11 = this->p_ucell->latvec.e11; 
+	    this->p_pw_wfc->latvec.e12 = this->p_ucell->latvec.e12; 
+	    this->p_pw_wfc->latvec.e13 = this->p_ucell->latvec.e13;
+            this->p_pw_wfc->latvec.e21 = this->p_ucell->latvec.e21; 
+	    this->p_pw_wfc->latvec.e22 = this->p_ucell->latvec.e22; 
+	    this->p_pw_wfc->latvec.e23 = this->p_ucell->latvec.e23;
+            this->p_pw_wfc->latvec.e31 = this->p_ucell->latvec.e31; 
+	    this->p_pw_wfc->latvec.e32 = this->p_ucell->latvec.e32; 
+	    this->p_pw_wfc->latvec.e33 = this->p_ucell->latvec.e33;
             this->p_pw_wfc->G = this->p_ucell->G;
             this->p_pw_wfc->GT = this->p_ucell->GT;
             this->p_pw_wfc->GGT = this->p_ucell->GGT;
             this->p_pw_wfc->lat0 = this->p_ucell->lat0;
             this->p_pw_wfc->tpiba = 2.0 * M_PI / this->p_ucell->lat0;
             this->p_pw_wfc->tpiba2 = this->p_pw_wfc->tpiba * this->p_pw_wfc->tpiba;
-            if(this->p_pw_wfc->kvec_c != nullptr) { delete[] this->p_pw_wfc->kvec_c;
-}
+            if(this->p_pw_wfc->kvec_c != nullptr)
+            {
+                delete[] this->p_pw_wfc->kvec_c;
+            }
             this->p_pw_wfc->kvec_c = new ModuleBase::Vector3<double>[1];
             this->p_pw_wfc->kvec_c[0] = {0.0, 0.0, 0.0};
-            if(this->p_pw_wfc->kvec_d != nullptr) { delete[] this->p_pw_wfc->kvec_d;
-}
+            if(this->p_pw_wfc->kvec_d != nullptr)
+            {
+                delete[] this->p_pw_wfc->kvec_d;
+            }
             this->p_pw_wfc->kvec_d = new ModuleBase::Vector3<double>[1];
             this->p_pw_wfc->kvec_d[0] = {0.0, 0.0, 0.0};
 
-            this->p_pspot_vnl->lmaxkb = 1;
+            this->lmaxkb = 1;
 
-            this->p_kv->ik2iktot.resize(1);
-            this->p_kv->ik2iktot[0] = 0;
+            this->ik2iktot_.resize(1);
+            this->ik2iktot_[0] = 0;
+            this->nkstot_ = 1;
 
         }
         void TearDown() override
@@ -267,37 +307,41 @@ class PsiIntializerUnitTest : public ::testing::Test {
             delete this->p_sf;
             delete this->p_pw_wfc;
             delete this->p_ucell;
-            delete this->p_pspot_vnl;
-            delete this->p_kv;
          }
 };
 
-TEST_F(PsiIntializerUnitTest, ConstructorRandom) {
+TEST_F(PsiIntializerUnitTest, ConstructorRandom)
+{
     this->psi_init = new psi_init_random<std::complex<double>>();
     EXPECT_EQ("random", this->psi_init->method());
 }
 
-TEST_F(PsiIntializerUnitTest, ConstructorAtomic) {
+TEST_F(PsiIntializerUnitTest, ConstructorAtomic)
+{
     this->psi_init = new psi_init_atomic<std::complex<double>>();
     EXPECT_EQ("atomic", this->psi_init->method());
 }
 
-TEST_F(PsiIntializerUnitTest, ConstructorAtomicRandom) {
+TEST_F(PsiIntializerUnitTest, ConstructorAtomicRandom)
+{
     this->psi_init = new psi_init_atomic_random<std::complex<double>>();
     EXPECT_EQ("atomic+random", this->psi_init->method());
 }
 
-TEST_F(PsiIntializerUnitTest, ConstructorNao) {
+TEST_F(PsiIntializerUnitTest, ConstructorNao)
+{
     this->psi_init = new psi_init_nao<std::complex<double>>();
     EXPECT_EQ("nao", this->psi_init->method());
 }
 
-TEST_F(PsiIntializerUnitTest, ConstructorNaoRandom) {
+TEST_F(PsiIntializerUnitTest, ConstructorNaoRandom)
+{
     this->psi_init = new psi_init_nao_random<std::complex<double>>();
     EXPECT_EQ("nao+random", this->psi_init->method());
 }
 
-TEST_F(PsiIntializerUnitTest, CastToT) {
+TEST_F(PsiIntializerUnitTest, CastToT)
+{
     this->psi_init = new psi_init_random<std::complex<double>>();
     std::complex<double> cd = {1.0, 2.0};
     std::complex<float> cf = {1.0, 2.0};
@@ -309,224 +353,277 @@ TEST_F(PsiIntializerUnitTest, CastToT) {
     EXPECT_EQ(this->psi_init->template cast_to_T<float>(cd), f);
 }
 
-TEST_F(PsiIntializerUnitTest, CalPsigRandom) {
-    PARAM.input.init_wfc = "random";
+TEST_F(PsiIntializerUnitTest, CalPsigRandom)
+{
     this->psi_init = new psi_init_random<std::complex<double>>();
-    this->psi_init->initialize(this->p_sf, 
-                               this->p_pw_wfc, 
-                               this->p_ucell, 
-                               this->p_kv, 
+    this->psi_init->initialize(this->p_pw_wfc,
+                               this->p_ucell,
+                               this->ik2iktot_,
+                               this->nkstot_,
                                this->random_seed,
-                               this->p_pspot_vnl,
-                               GlobalV::MY_RANK);
-    this->psi_init->tabulate(); // always: new, initialize, tabulate, allocate, proj_ao_onkG
+                               this->lmaxkb,
+                               GlobalV::MY_RANK,
+                               this->npol_,
+                               this->nbands_);
+    this->psi_init->tabulate();
     const int nbands_start = this->psi_init->nbands_start();
-    const int nbasis = this->p_pw_wfc->npwk_max * PARAM.globalv.npol;
+    const int nbasis = this->p_pw_wfc->npwk_max * this->npol_;
     psi::Psi<std::complex<double>>* psi = new psi::Psi<std::complex<double>>(1, nbands_start, nbasis, nbasis, true);
     this->psi_init->init_psig(psi->get_pointer(), 0);
     EXPECT_NEAR(-0.66187696761064307, psi->operator()(0,0,0).real(), 1e-4);
     delete psi;
 }
 
-TEST_F(PsiIntializerUnitTest, CalPsigAtomic) {
-    PARAM.input.init_wfc = "atomic";
-    this->psi_init = new psi_init_atomic<std::complex<double>>();
-    this->psi_init->initialize(this->p_sf, 
-                               this->p_pw_wfc, 
-                               this->p_ucell, 
-                               this->p_kv, 
+TEST_F(PsiIntializerUnitTest, CalPsigAtomic)
+{
+    psi_init_atomic<std::complex<double>>* atomic_initer = new psi_init_atomic<std::complex<double>>();
+    atomic_initer->prepare_params(this->nqx_, this->dq_, this->nspin_, 
+		    this->domag_, this->domag_z_, this->pseudo_mesh_,
+		    *this->p_sf);
+    this->psi_init = atomic_initer;
+    this->psi_init->initialize(this->p_pw_wfc,
+                               this->p_ucell,
+                               this->ik2iktot_,
+                               this->nkstot_,
                                this->random_seed,
-                               this->p_pspot_vnl,
-                               GlobalV::MY_RANK);
-    this->psi_init->tabulate(); // always: new, initialize, tabulate, allocate, proj_ao_onkG
+                               this->lmaxkb,
+                               GlobalV::MY_RANK,
+                               this->npol_,
+                               this->nbands_);
+    this->psi_init->tabulate();
     const int nbands_start = this->psi_init->nbands_start();
-    const int nbasis = this->p_pw_wfc->npwk_max * PARAM.globalv.npol;
+    const int nbasis = this->p_pw_wfc->npwk_max * this->npol_;
     psi::Psi<std::complex<double>>* psi = new psi::Psi<std::complex<double>>(1, nbands_start, nbasis, nbasis, true);
     this->psi_init->init_psig(psi->get_pointer(), 0);
     EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12);
     delete psi;
 }
 
-TEST_F(PsiIntializerUnitTest, CalPsigAtomicSoc) {
-    PARAM.input.init_wfc = "atomic";
-    PARAM.input.nspin = 4;
-    PARAM.sys.npol = 2;
+TEST_F(PsiIntializerUnitTest, CalPsigAtomicSoc)
+{
+    int nspin_save = this->nspin_;
+    int npol_save = this->npol_;
+    this->nspin_ = 4;
+    this->npol_ = 2;
     this->p_ucell->atoms[0].ncpp.has_so = false;
     this->p_ucell->natomwfc *= 2;
-    this->psi_init = new psi_init_atomic<std::complex<double>>();
-    this->psi_init->initialize(this->p_sf, 
-                               this->p_pw_wfc, 
-                               this->p_ucell, 
-                               this->p_kv, 
+    psi_init_atomic<std::complex<double>>* atomic_initer = new psi_init_atomic<std::complex<double>>();
+    atomic_initer->prepare_params(this->nqx_, this->dq_, this->nspin_, this->domag_, this->domag_z_, this->pseudo_mesh_,
+                                  *this->p_sf);
+    this->psi_init = atomic_initer;
+    this->psi_init->initialize(this->p_pw_wfc,
+                               this->p_ucell,
+                               this->ik2iktot_,
+                               this->nkstot_,
                                this->random_seed,
-                               this->p_pspot_vnl,
-                               GlobalV::MY_RANK);
-    this->psi_init->tabulate(); // always: new, initialize, tabulate, allocate, proj_ao_onkG
+                               this->lmaxkb,
+                               GlobalV::MY_RANK,
+                               this->npol_,
+                               this->nbands_);
+    this->psi_init->tabulate();
     const int nbands_start = this->psi_init->nbands_start();
-    const int nbasis = this->p_pw_wfc->npwk_max * PARAM.globalv.npol;
+    const int nbasis = this->p_pw_wfc->npwk_max * this->npol_;
     psi::Psi<std::complex<double>>* psi = new psi::Psi<std::complex<double>>(1, nbands_start, nbasis, nbasis, true);
     this->psi_init->init_psig(psi->get_pointer(), 0);
     EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12);
-    PARAM.input.nspin = 1;
-    PARAM.sys.npol = 1;
+    this->nspin_ = nspin_save;
+    this->npol_ = npol_save;
     this->p_ucell->atoms[0].ncpp.has_so = false;
     this->p_ucell->natomwfc /= 2;
     delete psi;
 }
 
-TEST_F(PsiIntializerUnitTest, CalPsigAtomicSocHasSo) {
-    PARAM.input.init_wfc = "atomic";
-    PARAM.input.nspin = 4;
-    PARAM.sys.npol = 2;
+TEST_F(PsiIntializerUnitTest, CalPsigAtomicSocHasSo)
+{
+    int nspin_save = this->nspin_;
+    int npol_save = this->npol_;
+    this->nspin_ = 4;
+    this->npol_ = 2;
     this->p_ucell->atoms[0].ncpp.has_so = true;
     this->p_ucell->natomwfc *= 2;
-    this->psi_init = new psi_init_atomic<std::complex<double>>();
-    this->psi_init->initialize(this->p_sf, 
-                               this->p_pw_wfc, 
-                               this->p_ucell, 
-                               this->p_kv, 
+    psi_init_atomic<std::complex<double>>* atomic_initer = new psi_init_atomic<std::complex<double>>();
+    atomic_initer->prepare_params(this->nqx_, this->dq_, this->nspin_, this->domag_, this->domag_z_, this->pseudo_mesh_,
+                                  *this->p_sf);
+    this->psi_init = atomic_initer;
+    this->psi_init->initialize(this->p_pw_wfc,
+                               this->p_ucell,
+                               this->ik2iktot_,
+                               this->nkstot_,
                                this->random_seed,
-                               this->p_pspot_vnl,
-                               GlobalV::MY_RANK);
-    this->psi_init->tabulate(); // always: new, initialize, tabulate, allocate, proj_ao_onkG
+                               this->lmaxkb,
+                               GlobalV::MY_RANK,
+                               this->npol_,
+                               this->nbands_);
+    this->psi_init->tabulate();
     const int nbands_start = this->psi_init->nbands_start();
-    const int nbasis = this->p_pw_wfc->npwk_max * PARAM.globalv.npol;
+    const int nbasis = this->p_pw_wfc->npwk_max * this->npol_;
     psi::Psi<std::complex<double>>* psi = new psi::Psi<std::complex<double>>(1, nbands_start, nbasis, nbasis, true);
     this->psi_init->init_psig(psi->get_pointer(), 0);
     EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12);
-    PARAM.input.nspin = 1;
-    PARAM.sys.npol = 1;
+    this->nspin_ = nspin_save;
+    this->npol_ = npol_save;
     this->p_ucell->atoms[0].ncpp.has_so = false;
     this->p_ucell->natomwfc /= 2;
     delete psi;
 }
 
-TEST_F(PsiIntializerUnitTest, CalPsigAtomicRandom) {
-    PARAM.input.init_wfc = "atomic+random";
-    this->psi_init = new psi_init_atomic_random<std::complex<double>>();
-    this->psi_init->initialize(this->p_sf, 
-                               this->p_pw_wfc, 
-                               this->p_ucell, 
-                               this->p_kv, 
+TEST_F(PsiIntializerUnitTest, CalPsigAtomicRandom)
+{
+    psi_init_atomic_random<std::complex<double>>* atomic_rand_initer = new psi_init_atomic_random<std::complex<double>>();
+    atomic_rand_initer->prepare_params(this->nqx_, this->dq_, this->nspin_, this->domag_, this->domag_z_, this->pseudo_mesh_,
+                                       *this->p_sf);
+    this->psi_init = atomic_rand_initer;
+    this->psi_init->initialize(this->p_pw_wfc,
+                               this->p_ucell,
+                               this->ik2iktot_,
+                               this->nkstot_,
                                this->random_seed,
-                               this->p_pspot_vnl,
-                               GlobalV::MY_RANK);
-    this->psi_init->tabulate(); // always: new, initialize, tabulate, allocate, proj_ao_onkG
+                               this->lmaxkb,
+                               GlobalV::MY_RANK,
+                               this->npol_,
+                               this->nbands_);
+    this->psi_init->tabulate();
     const int nbands_start = this->psi_init->nbands_start();
-    const int nbasis = this->p_pw_wfc->npwk_max * PARAM.globalv.npol;
+    const int nbasis = this->p_pw_wfc->npwk_max * this->npol_;
     psi::Psi<std::complex<double>>* psi = new psi::Psi<std::complex<double>>(1, nbands_start, nbasis, nbasis, true);
     this->psi_init->init_psig(psi->get_pointer(), 0);
     EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12);
     delete psi;
 }
 
-TEST_F(PsiIntializerUnitTest, CalPsigNao) {
-    PARAM.input.init_wfc = "nao";
-    this->psi_init = new psi_init_nao<std::complex<double>>();
-    this->psi_init->initialize(this->p_sf, 
-                               this->p_pw_wfc, 
-                               this->p_ucell, 
-                               this->p_kv, 
+TEST_F(PsiIntializerUnitTest, CalPsigNao)
+{
+    psi_init_nao<std::complex<double>>* nao_initer = new psi_init_nao<std::complex<double>>();
+    nao_initer->prepare_params(this->nqx_, this->dq_, this->nspin_, this->orbital_dir_, *this->p_sf);
+    this->psi_init = nao_initer;
+    this->psi_init->initialize(this->p_pw_wfc,
+                               this->p_ucell,
+                               this->ik2iktot_,
+                               this->nkstot_,
                                this->random_seed,
-                               this->p_pspot_vnl,
-                               GlobalV::MY_RANK);
-    this->psi_init->tabulate(); // always: new, initialize, tabulate, allocate, proj_ao_onkG
+                               this->lmaxkb,
+                               GlobalV::MY_RANK,
+                               this->npol_,
+                               this->nbands_);
+    this->psi_init->tabulate();
     const int nbands_start = this->psi_init->nbands_start();
-    const int nbasis = this->p_pw_wfc->npwk_max * PARAM.globalv.npol;
+    const int nbasis = this->p_pw_wfc->npwk_max * this->npol_;
     psi::Psi<std::complex<double>>* psi = new psi::Psi<std::complex<double>>(1, nbands_start, nbasis, nbasis, true);
     this->psi_init->init_psig(psi->get_pointer(), 0);
     EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12);
     delete psi;
 }
 
-TEST_F(PsiIntializerUnitTest, CalPsigNaoRandom) {
-    PARAM.input.init_wfc = "nao+random";
-    this->psi_init = new psi_init_nao_random<std::complex<double>>();
-    this->psi_init->initialize(this->p_sf, 
-                               this->p_pw_wfc, 
-                               this->p_ucell, 
-                               this->p_kv, 
+TEST_F(PsiIntializerUnitTest, CalPsigNaoRandom)
+{
+    psi_init_nao_random<std::complex<double>>* nao_rand_initer = new psi_init_nao_random<std::complex<double>>();
+    nao_rand_initer->prepare_params(this->nqx_, this->dq_, this->nspin_, this->orbital_dir_, *this->p_sf);
+    this->psi_init = nao_rand_initer;
+    this->psi_init->initialize(this->p_pw_wfc,
+                               this->p_ucell,
+                               this->ik2iktot_,
+                               this->nkstot_,
                                this->random_seed,
-                               this->p_pspot_vnl,
-                               GlobalV::MY_RANK);
-    this->psi_init->tabulate(); // always: new, initialize, tabulate, allocate, proj_ao_onkG
+                               this->lmaxkb,
+                               GlobalV::MY_RANK,
+                               this->npol_,
+                               this->nbands_);
+    this->psi_init->tabulate();
     const int nbands_start = this->psi_init->nbands_start();
-    const int nbasis = this->p_pw_wfc->npwk_max * PARAM.globalv.npol;
+    const int nbasis = this->p_pw_wfc->npwk_max * this->npol_;
     psi::Psi<std::complex<double>>* psi = new psi::Psi<std::complex<double>>(1, nbands_start, nbasis, nbasis, true);
     this->psi_init->init_psig(psi->get_pointer(), 0);
     EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12);
     delete psi;
 }
 
-TEST_F(PsiIntializerUnitTest, CalPsigNaoSoc) {
-    PARAM.input.init_wfc = "nao";
-    PARAM.input.nspin = 4;
-    PARAM.sys.npol = 2;
+TEST_F(PsiIntializerUnitTest, CalPsigNaoSoc)
+{
+    int nspin_save = this->nspin_;
+    int npol_save = this->npol_;
+    this->nspin_ = 4;
+    this->npol_ = 2;
     this->p_ucell->atoms[0].ncpp.has_so = false;
-    PARAM.sys.domag = false;
-    PARAM.sys.domag_z = false;
-    this->psi_init = new psi_init_nao<std::complex<double>>();
-    this->psi_init->initialize(this->p_sf, 
-                               this->p_pw_wfc, 
-                               this->p_ucell, 
-                               this->p_kv, 
+    psi_init_nao<std::complex<double>>* nao_initer = new psi_init_nao<std::complex<double>>();
+    nao_initer->prepare_params(this->nqx_, this->dq_, this->nspin_, this->orbital_dir_, *this->p_sf);
+    this->psi_init = nao_initer;
+    this->psi_init->initialize(this->p_pw_wfc,
+                               this->p_ucell,
+                               this->ik2iktot_,
+                               this->nkstot_,
                                this->random_seed,
-                               this->p_pspot_vnl,
-                               GlobalV::MY_RANK);
-    this->psi_init->tabulate(); // always: new, initialize, tabulate, allocate, proj_ao_onkG
+                               this->lmaxkb,
+                               GlobalV::MY_RANK,
+                               this->npol_,
+                               this->nbands_);
+    this->psi_init->tabulate();
     const int nbands_start = this->psi_init->nbands_start();
-    const int nbasis = this->p_pw_wfc->npwk_max * PARAM.globalv.npol;
+    const int nbasis = this->p_pw_wfc->npwk_max * this->npol_;
     psi::Psi<std::complex<double>>* psi = new psi::Psi<std::complex<double>>(1, nbands_start, nbasis, nbasis, true);
     this->psi_init->init_psig(psi->get_pointer(), 0);
     EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12);
+    this->nspin_ = nspin_save;
+    this->npol_ = npol_save;
     delete psi;
 }
 
-TEST_F(PsiIntializerUnitTest, CalPsigNaoSocHasSo) {
-    PARAM.input.init_wfc = "nao";
-    PARAM.input.nspin = 4;
-    PARAM.sys.npol = 2;
+TEST_F(PsiIntializerUnitTest, CalPsigNaoSocHasSo)
+{
+    int nspin_save = this->nspin_;
+    int npol_save = this->npol_;
+    this->nspin_ = 4;
+    this->npol_ = 2;
     this->p_ucell->atoms[0].ncpp.has_so = true;
-    PARAM.sys.domag = false;
-    PARAM.sys.domag_z = false;
-    this->psi_init = new psi_init_nao<std::complex<double>>();
-    this->psi_init->initialize(this->p_sf, 
-                               this->p_pw_wfc, 
-                               this->p_ucell, 
-                               this->p_kv, 
+    psi_init_nao<std::complex<double>>* nao_initer = new psi_init_nao<std::complex<double>>();
+    nao_initer->prepare_params(this->nqx_, this->dq_, this->nspin_, this->orbital_dir_, *this->p_sf);
+    this->psi_init = nao_initer;
+    this->psi_init->initialize(this->p_pw_wfc,
+                               this->p_ucell,
+                               this->ik2iktot_,
+                               this->nkstot_,
                                this->random_seed,
-                               this->p_pspot_vnl,
-                               GlobalV::MY_RANK);
-    this->psi_init->tabulate(); // always: new, initialize, tabulate, allocate, proj_ao_onkG
+                               this->lmaxkb,
+                               GlobalV::MY_RANK,
+                               this->npol_,
+                               this->nbands_);
+    this->psi_init->tabulate();
     const int nbands_start = this->psi_init->nbands_start();
-    const int nbasis = this->p_pw_wfc->npwk_max * PARAM.globalv.npol;
+    const int nbasis = this->p_pw_wfc->npwk_max * this->npol_;
     psi::Psi<std::complex<double>>* psi = new psi::Psi<std::complex<double>>(1, nbands_start, nbasis, nbasis, true);
     this->psi_init->init_psig(psi->get_pointer(), 0);
     EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12);
+    this->nspin_ = nspin_save;
+    this->npol_ = npol_save;
     delete psi;
 }
 
-TEST_F(PsiIntializerUnitTest, CalPsigNaoSocHasSoDOMAG) {
-    PARAM.input.init_wfc = "nao";
-    PARAM.input.nspin = 4;
-    PARAM.sys.npol = 2;
+TEST_F(PsiIntializerUnitTest, CalPsigNaoSocHasSoDOMAG)
+{
+    int nspin_save = this->nspin_;
+    int npol_save = this->npol_;
+    this->nspin_ = 4;
+    this->npol_ = 2;
     this->p_ucell->atoms[0].ncpp.has_so = true;
-    PARAM.sys.domag = true;
-    PARAM.sys.domag_z = false;
-    this->psi_init = new psi_init_nao<std::complex<double>>();
-    this->psi_init->initialize(this->p_sf, 
-                               this->p_pw_wfc, 
-                               this->p_ucell, 
-                               this->p_kv, 
+    psi_init_nao<std::complex<double>>* nao_initer = new psi_init_nao<std::complex<double>>();
+    nao_initer->prepare_params(this->nqx_, this->dq_, this->nspin_, this->orbital_dir_, *this->p_sf);
+    this->psi_init = nao_initer;
+    this->psi_init->initialize(this->p_pw_wfc,
+                               this->p_ucell,
+                               this->ik2iktot_,
+                               this->nkstot_,
                                this->random_seed,
-                               this->p_pspot_vnl,
-                               GlobalV::MY_RANK);
-    this->psi_init->tabulate(); // always: new, initialize, tabulate, allocate, proj_ao_onkG
+                               this->lmaxkb,
+                               GlobalV::MY_RANK,
+                               this->npol_,
+                               this->nbands_);
+    this->psi_init->tabulate();
     const int nbands_start = this->psi_init->nbands_start();
-    const int nbasis = this->p_pw_wfc->npwk_max * PARAM.globalv.npol;
+    const int nbasis = this->p_pw_wfc->npwk_max * this->npol_;
     psi::Psi<std::complex<double>>* psi = new psi::Psi<std::complex<double>>(1, nbands_start, nbasis, nbasis, true);
     this->psi_init->init_psig(psi->get_pointer(), 0);
     EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12);
+    this->nspin_ = nspin_save;
+    this->npol_ = npol_save;
     delete psi;
 }
 
@@ -547,4 +644,4 @@ int main(int argc, char** argv)
 #endif
 
     return result;
-}
\ No newline at end of file
+}
diff --git a/source/source_pw/module_pwdft/test/soc_test.cpp b/source/source_psi/test/soc_base_test.cpp
similarity index 92%
rename from source/source_pw/module_pwdft/test/soc_test.cpp
rename to source/source_psi/test/soc_base_test.cpp
index 9ea9bb53a92..98b13819384 100644
--- a/source/source_pw/module_pwdft/test/soc_test.cpp
+++ b/source/source_psi/test/soc_base_test.cpp
@@ -5,16 +5,16 @@
 #include <complex>
 
 /************************************************
- *  unit test of class Soc and Fcoef
+ *  unit test of class SocBase and Fcoef
  ***********************************************/
 
 /**
  * - Tested Functions:
  *   - Fcoef::create to create a 5 dimensional array of complex numbers
- *   - Soc::set_fcoef to set the fcoef array
- *   - Soc::spinor to calculate the spinor
- *   - Soc::rot_ylm to calculate the rotation matrix
- *   - Soc::sph_ind to calculate the m index of the spherical harmonics
+ *   - SocBase::set_fcoef to set the fcoef array
+ *   - SocBase::spinor to calculate the spinor
+ *   - SocBase::rot_ylm to calculate the rotation matrix
+ *   - SocBase::sph_ind to calculate the m index of the spherical harmonics
 */
 
 //compare two complex by using EXPECT_DOUBLE_EQ()
@@ -25,7 +25,7 @@ void EXPECT_COMPLEX_DOUBLE_EQ(const std::complex<double>& a,const std::complex<d
 }
 
 #define private public
-#include "source_pw/module_pwdft/soc.h"
+#include "../soc_base.h"
 #undef private
 class FcoefTest : public testing::Test
 {
@@ -57,15 +57,15 @@ TEST_F(FcoefTest, CreateWarning)
                   EXPECT_EQ(fcoef.p, nullptr);
 }
 
-class SocTest : public testing::Test
+class SocBaseTest : public testing::Test
 {
 protected:
-                  Soc soc;
+                  SocBase soc;
                   Fcoef fcoef;
                   std::string output;
 };
 
-TEST_F(SocTest, Spinor)
+TEST_F(SocBaseTest, Spinor)
 {
                   EXPECT_DOUBLE_EQ(soc.spinor(0, 0.5, 0, 0), 1.0);
                   EXPECT_DOUBLE_EQ(soc.spinor(0, 0.5, 0, 1), 0.0);
@@ -91,7 +91,7 @@ TEST_F(SocTest, Spinor)
                   EXPECT_THAT(output, testing::HasSubstr("m not allowed"));
 }
 
-TEST_F(SocTest, SphInd)
+TEST_F(SocBaseTest, SphInd)
 {
                   EXPECT_EQ(soc.sph_ind(0, 0.5, 0, 0), 0);
                   EXPECT_EQ(soc.sph_ind(0, 0.5, 0, 1), 0);
@@ -115,7 +115,7 @@ TEST_F(SocTest, SphInd)
                   EXPECT_THAT(output, testing::HasSubstr("m not allowed"));
 }
 
-TEST_F(SocTest, RotYlm)
+TEST_F(SocBaseTest, RotYlm)
 {
                   soc.rot_ylm(0);
                   EXPECT_NE(soc.p_rot, nullptr);
@@ -138,7 +138,7 @@ TEST_F(SocTest, RotYlm)
                   EXPECT_DOUBLE_EQ(rot.imag(),1.0/sqrt(2.0));
 }
 
-TEST_F(SocTest, SetFcoef)
+TEST_F(SocBaseTest, SetFcoef)
 {
                   soc.fcoef.create(1,1,1);
                   soc.rot_ylm(0);
diff --git a/source/source_psi/test/support/atomic_new b/source/source_psi/test/support/atomic_new
deleted file mode 100644
index af7a7d973a4..00000000000
--- a/source/source_psi/test/support/atomic_new
+++ /dev/null
@@ -1,11 +0,0 @@
-INPUT_PARAMETERS
-pseudo_dir .
-orbital_dir .
-basis_type pw
-ecutwfc 60
-scf_thr 1e-8
-scf_nmax 1
-
-init_wfc atomic
-psi_initializer 1
-ks_solver dav
\ No newline at end of file
diff --git a/source/source_psi/test/support/nao_new b/source/source_psi/test/support/nao_new
deleted file mode 100644
index b0b7789a2d5..00000000000
--- a/source/source_psi/test/support/nao_new
+++ /dev/null
@@ -1,11 +0,0 @@
-INPUT_PARAMETERS
-pseudo_dir .
-orbital_dir .
-basis_type pw
-ecutwfc 60
-scf_thr 1e-8
-scf_nmax 1
-
-init_wfc nao
-psi_initializer 1
-ks_solver dav
\ No newline at end of file
diff --git a/source/source_psi/test/support/random_new b/source/source_psi/test/support/random_new
deleted file mode 100644
index 87a70829dd9..00000000000
--- a/source/source_psi/test/support/random_new
+++ /dev/null
@@ -1,11 +0,0 @@
-INPUT_PARAMETERS
-pseudo_dir .
-orbital_dir .
-basis_type pw
-ecutwfc 60
-scf_thr 1e-8
-scf_nmax 1
-
-init_wfc random
-psi_initializer 1
-ks_solver dav
\ No newline at end of file
diff --git a/source/source_pw/module_pwdft/CMakeLists.txt b/source/source_pw/module_pwdft/CMakeLists.txt
index dec6cd86d8e..aa7b00a881e 100644
--- a/source/source_pw/module_pwdft/CMakeLists.txt
+++ b/source/source_pw/module_pwdft/CMakeLists.txt
@@ -41,7 +41,6 @@ list(APPEND objects
     vnl_pw_grad.cpp
     structure_factor.cpp
     structure_factor_k.cpp
-    soc.cpp
     parallel_grid.cpp
     elecond.cpp
     fs_nonlocal_tools.cpp
diff --git a/source/source_pw/module_pwdft/test/CMakeLists.txt b/source/source_pw/module_pwdft/test/CMakeLists.txt
index 08a5c9dd8e4..b2b690d82b1 100644
--- a/source/source_pw/module_pwdft/test/CMakeLists.txt
+++ b/source/source_pw/module_pwdft/test/CMakeLists.txt
@@ -6,7 +6,7 @@ abacus_disable_feature_definitions(__EXX)
 AddTest(
   TARGET MODULE_PW_pwdft_soc
   LIBS parameter
-  SOURCES soc_test.cpp ../soc.cpp
+  SOURCES ../../../source_psi/test/soc_base_test.cpp ../../../source_psi/soc_base.cpp
   	../../../source_base/global_variable.cpp
   	../../../source_base/global_function.cpp
 	../../../source_base/global_file.cpp
diff --git a/source/source_pw/module_pwdft/vnl_pw.h b/source/source_pw/module_pwdft/vnl_pw.h
index 6282b138657..7547424d0fd 100644
--- a/source/source_pw/module_pwdft/vnl_pw.h
+++ b/source/source_pw/module_pwdft/vnl_pw.h
@@ -6,7 +6,7 @@
 #include "source_base/intarray.h"
 #include "source_base/realarray.h"
 #include "source_cell/unitcell.h"
-#include "source_pw/module_pwdft/soc.h"
+#include "source_psi/soc_base.h"
 #include "source_pw/module_pwdft/structure_factor.h"
 #include "source_psi/psi.h"
 #ifdef __LCAO
@@ -202,7 +202,7 @@ class pseudopot_cell_vnl
 
     const ModulePW::PW_Basis_K* wfcpw = nullptr;
 
-    Soc soc;
+    SocBase soc;
 
     double omega_old = 0;
     bool use_gpu_ = false;
diff --git a/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/threshold b/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/threshold
new file mode 100644
index 00000000000..dc467b56d6f
--- /dev/null
+++ b/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/threshold
@@ -0,0 +1 @@
+threshold 1e-6