diff --git a/source/module_elecstate/elecstate_energy_terms.cpp b/source/module_elecstate/elecstate_energy_terms.cpp
index d820ba064e..535e9c83a3 100644
--- a/source/module_elecstate/elecstate_energy_terms.cpp
+++ b/source/module_elecstate/elecstate_energy_terms.cpp
@@ -26,12 +26,12 @@ double ElecState::get_etot_gatefield()
 
 double ElecState::get_solvent_model_Ael()
 {
-    return GlobalC::solvent_model.Ael;
+    return surchem::Ael;
 }
 
 double ElecState::get_solvent_model_Acav()
 {
-    return GlobalC::solvent_model.Acav;
+    return surchem::Acav;
 }
 
 double ElecState::get_dftu_energy()
diff --git a/source/module_elecstate/potentials/efield.cpp b/source/module_elecstate/potentials/efield.cpp
index 0297015095..ed602d5405 100644
--- a/source/module_elecstate/potentials/efield.cpp
+++ b/source/module_elecstate/potentials/efield.cpp
@@ -34,7 +34,7 @@ ModuleBase::matrix Efield::add_efield(const UnitCell& cell,
                                       const ModulePW::PW_Basis* rho_basis,
                                       const int& nspin,
                                       const double* const* const rho,
-                                      surchem& solvent)
+                                      const surchem& solvent)
 {
     ModuleBase::TITLE("Efield", "add_efield");
     ModuleBase::timer::tick("Efield", "add_efield");
@@ -202,7 +202,7 @@ double Efield::cal_elec_dipole(const UnitCell& cell,
 
 double Efield::cal_induced_dipole(const UnitCell& cell,
                                   const ModulePW::PW_Basis* rho_basis,
-                                  surchem& solvent,
+                                  const surchem& solvent,
                                   const double& bmod)
 {
     double induced_dipole = 0;
diff --git a/source/module_elecstate/potentials/efield.h b/source/module_elecstate/potentials/efield.h
index 79a4860b03..14db473155 100644
--- a/source/module_elecstate/potentials/efield.h
+++ b/source/module_elecstate/potentials/efield.h
@@ -1,9 +1,10 @@
 #ifndef EFIELD_H
 #define EFIELD_H
 
-#include "module_cell/unitcell.h"
 #include "module_basis/module_pw/pw_basis.h"
+#include "module_cell/unitcell.h"
 #include "module_hamilt_general/module_surchem/surchem.h"
+#include "module_parameter/parameter.h"
 
 namespace elecstate
 {
@@ -17,7 +18,7 @@ class Efield
                                          const ModulePW::PW_Basis* rho_basis,
                                          const int& nspin,
                                          const double* const* const rho,
-                                         surchem& solvent);
+                                         const surchem& solvent);
 
     static double cal_elec_dipole(const UnitCell& cell,
                                   const ModulePW::PW_Basis* rho_basis,
@@ -29,7 +30,7 @@ class Efield
 
     static double cal_induced_dipole(const UnitCell& cell,
                                      const ModulePW::PW_Basis* rho_basis,
-                                     surchem& solvent,
+                                     const surchem& solvent,
                                      const double& bmod);
 
     static double saw_function(const double &a, const double &b, const double &x);
@@ -59,7 +60,8 @@ namespace elecstate
 class PotEfield : public PotBase
 {
   public:
-    PotEfield(const ModulePW::PW_Basis* rho_basis_in, const UnitCell* ucell_in, bool dipole) : ucell_(ucell_in)
+    PotEfield(const ModulePW::PW_Basis* rho_basis_in, const UnitCell* ucell_in, const surchem* solvent_in, bool dipole)
+        : ucell_(ucell_in), solvent_(solvent_in)
     {
         this->rho_basis_ = rho_basis_in;
         if (!dipole)
@@ -81,7 +83,7 @@ class PotEfield : public PotBase
                                       const_cast<const ModulePW::PW_Basis*>(rho_basis_),
                                       PARAM.inp.nspin,
                                       nullptr,
-                                      GlobalC::solvent_model);
+                                      *solvent_);
         for (int ir = 0; ir < rho_basis_->nrxx; ++ir)
         {
             vl_pseudo[ir] += v_efield(0, ir);
@@ -94,11 +96,12 @@ class PotEfield : public PotBase
                                     const_cast<const ModulePW::PW_Basis*>(rho_basis_),
                                     v_eff.nr,
                                     chg->rho,
-                                    GlobalC::solvent_model);
+                                    *solvent_);
     }
 
   private:
-    const UnitCell *ucell_ = nullptr;
+    const UnitCell* ucell_ = nullptr;
+    const surchem* solvent_ = nullptr;
 };
 
 } // namespace elecstate
diff --git a/source/module_elecstate/potentials/potential_new.cpp b/source/module_elecstate/potentials/potential_new.cpp
index 48b90f65b4..a4443c46d8 100644
--- a/source/module_elecstate/potentials/potential_new.cpp
+++ b/source/module_elecstate/potentials/potential_new.cpp
@@ -23,9 +23,11 @@ Potential::Potential(const ModulePW::PW_Basis* rho_basis_in,
                      const UnitCell* ucell_in,
                      const ModuleBase::matrix* vloc_in,
                      Structure_Factor* structure_factors_in,
+                     surchem* solvent_in,
                      double* etxc_in,
                      double* vtxc_in)
-    : ucell_(ucell_in), vloc_(vloc_in), structure_factors_(structure_factors_in), etxc_(etxc_in), vtxc_(vtxc_in)
+    : ucell_(ucell_in), vloc_(vloc_in), structure_factors_(structure_factors_in), solvent_(solvent_in), etxc_(etxc_in),
+      vtxc_(vtxc_in)
 {
     this->rho_basis_ = rho_basis_in;
     this->rho_basis_smooth_ = rho_basis_smooth_in;
diff --git a/source/module_elecstate/potentials/potential_new.h b/source/module_elecstate/potentials/potential_new.h
index d1f8e41f78..fd6e087534 100644
--- a/source/module_elecstate/potentials/potential_new.h
+++ b/source/module_elecstate/potentials/potential_new.h
@@ -1,13 +1,14 @@
 #ifndef POTENTIALNEW_H
 #define POTENTIALNEW_H
 
-#include <vector>
-
 #include "module_base/complexmatrix.h"
+#include "module_hamilt_general/module_surchem/surchem.h"
 #include "module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h"
 #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h"
 #include "pot_base.h"
 
+#include <vector>
+
 namespace elecstate
 {
 /**
@@ -58,6 +59,7 @@ class Potential : public PotBase
               const UnitCell* ucell_in,
               const ModuleBase::matrix* vloc_in,
               Structure_Factor* structure_factors_in,
+              surchem* solvent_in,
               double* etxc_in,
               double* vtxc_in);
     ~Potential();
@@ -211,6 +213,7 @@ class Potential : public PotBase
     const UnitCell* ucell_ = nullptr;
     const ModuleBase::matrix* vloc_ = nullptr;
     Structure_Factor* structure_factors_ = nullptr;
+    surchem* solvent_ = nullptr;
 };
 
 } // namespace elecstate
diff --git a/source/module_elecstate/potentials/potential_types.cpp b/source/module_elecstate/potentials/potential_types.cpp
index bc15110798..7ca1bacc8a 100644
--- a/source/module_elecstate/potentials/potential_types.cpp
+++ b/source/module_elecstate/potentials/potential_types.cpp
@@ -47,11 +47,11 @@ PotBase* Potential::get_pot_type(const std::string& pot_type)
         return new PotSurChem(this->rho_basis_,
                               this->structure_factors_,
                               this->v_effective_fixed.data(),
-                              &GlobalC::solvent_model);
+                              this->solvent_);
     }
     else if (pot_type == "efield")
     {
-        return new PotEfield(this->rho_basis_, this->ucell_, PARAM.inp.dip_cor_flag);
+        return new PotEfield(this->rho_basis_, this->ucell_, this->solvent_, PARAM.inp.dip_cor_flag);
     }
     else if (pot_type == "gatefield")
     {
diff --git a/source/module_elecstate/test/potential_new_test.cpp b/source/module_elecstate/test/potential_new_test.cpp
index b38361d76f..7115886bff 100644
--- a/source/module_elecstate/test/potential_new_test.cpp
+++ b/source/module_elecstate/test/potential_new_test.cpp
@@ -38,6 +38,12 @@ Charge::Charge()
 Charge::~Charge()
 {
 }
+surchem::surchem()
+{
+}
+surchem::~surchem()
+{
+}
 
 namespace elecstate
 {
@@ -101,6 +107,7 @@ class PotentialNewTest : public ::testing::Test
     UnitCell* ucell = nullptr;
     ModuleBase::matrix* vloc = nullptr;
     Structure_Factor* structure_factors = nullptr;
+    surchem* solvent = nullptr;
     double* etxc = nullptr;
     double* vtxc = nullptr;
     elecstate::Potential* pot = nullptr;
@@ -111,43 +118,56 @@ class PotentialNewTest : public ::testing::Test
         ucell = new UnitCell;
         vloc = new ModuleBase::matrix;
         structure_factors = new Structure_Factor();
+        solvent = new surchem();
         etxc = new double;
         vtxc = new double;
         elecstate::Set_GlobalV_Default();
     }
     virtual void TearDown()
     {
-        if (rhopw != nullptr) {
+        if (rhopw != nullptr)
+        {
             delete rhopw;
-}
-        if (rhodpw != nullptr) {
+        }
+        if (rhodpw != nullptr)
+        {
             delete rhodpw;
-}
-        if (ucell != nullptr) {
+        }
+        if (ucell != nullptr)
+        {
             delete ucell;
-}
-        if (vloc != nullptr) {
+        }
+        if (vloc != nullptr)
+        {
             delete vloc;
-}
-        if (structure_factors != nullptr) {
+        }
+        if (structure_factors != nullptr)
+        {
             delete structure_factors;
-}
-        if (etxc != nullptr) {
+        }
+        if (solvent != nullptr)
+        {
+            delete solvent;
+        }
+        if (etxc != nullptr)
+        {
             delete etxc;
-}
-        if (vtxc != nullptr) {
+        }
+        if (vtxc != nullptr)
+        {
             delete vtxc;
-}
-        if (pot != nullptr) {
+        }
+        if (pot != nullptr)
+        {
             delete pot;
-}
+        }
     }
 };
 
 TEST_F(PotentialNewTest, ConstructorCPUDouble)
 {
     rhopw->nrxx = 100;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     EXPECT_TRUE(pot->fixed_mode);
     EXPECT_TRUE(pot->dynamic_mode);
     EXPECT_EQ(pot->v_effective_fixed.size(), 100);
@@ -159,7 +179,7 @@ TEST_F(PotentialNewTest, ConstructorCPUSingle)
 {
     rhopw->nrxx = 100;
     PARAM.input.precision = "single";
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     EXPECT_TRUE(pot->fixed_mode);
     EXPECT_TRUE(pot->dynamic_mode);
     EXPECT_EQ(pot->v_effective_fixed.size(), 100);
@@ -170,7 +190,7 @@ TEST_F(PotentialNewTest, ConstructorCPUSingle)
 TEST_F(PotentialNewTest, ConstructorNRXX0)
 {
     rhopw->nrxx = 0;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     EXPECT_TRUE(pot->fixed_mode);
     EXPECT_TRUE(pot->dynamic_mode);
 }
@@ -179,7 +199,7 @@ TEST_F(PotentialNewTest, ConstructorXC3)
 {
     elecstate::tmp_xc_func_type = 3;
     rhopw->nrxx = 100;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     EXPECT_TRUE(pot->fixed_mode);
     EXPECT_TRUE(pot->dynamic_mode);
     EXPECT_EQ(pot->v_effective_fixed.size(), 100);
@@ -194,7 +214,7 @@ TEST_F(PotentialNewTest, ConstructorGPUDouble)
     // this is just a trivial call to the GPU code
     rhopw->nrxx = 100;
     PARAM.input.device = "gpu";
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     EXPECT_TRUE(pot->fixed_mode);
     EXPECT_TRUE(pot->dynamic_mode);
     EXPECT_EQ(pot->v_effective_fixed.size(), 100);
@@ -208,7 +228,7 @@ TEST_F(PotentialNewTest, ConstructorGPUSingle)
     rhopw->nrxx = 100;
     PARAM.input.device = "gpu";
     PARAM.input.precision = "single";
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     EXPECT_TRUE(pot->fixed_mode);
     EXPECT_TRUE(pot->dynamic_mode);
     EXPECT_EQ(pot->v_effective_fixed.size(), 100);
@@ -251,7 +271,7 @@ TEST_F(PotentialNewTest, CalFixedV)
 {
     // construct potential
     rhopw->nrxx = 100;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
     std::vector<std::string> compnents_list = {
         "local",
@@ -279,7 +299,7 @@ TEST_F(PotentialNewTest, CalVeff)
 {
     // construct potential
     rhopw->nrxx = 100;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
     std::vector<std::string> compnents_list = {
         "local",
@@ -309,7 +329,7 @@ TEST_F(PotentialNewTest, UpdateFromCharge)
 {
     // construct potential
     rhopw->nrxx = 100;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
     std::vector<std::string> compnents_list = {
         "local",
@@ -337,7 +357,7 @@ TEST_F(PotentialNewTest, InitPot)
 {
     // construct potential
     rhopw->nrxx = 100;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
     std::vector<std::string> compnents_list = {
         "local",
@@ -365,7 +385,7 @@ TEST_F(PotentialNewTest, GetVnew)
 {
     // construct potential
     rhopw->nrxx = 100;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
     std::vector<std::string> compnents_list = {
         "local",
@@ -394,7 +414,7 @@ TEST_F(PotentialNewTest, GetEffectiveVmatrix)
 {
     // construct potential
     rhopw->nrxx = 100;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
     ModuleBase::matrix v_eff_tmp = pot->get_effective_v();
     const ModuleBase::matrix v_eff_tmp_const = pot->get_effective_v();
@@ -416,7 +436,7 @@ TEST_F(PotentialNewTest, GetEffectiveVarray)
 {
     // construct potential
     rhopw->nrxx = 100;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
     double* v_eff_tmp = pot->get_effective_v(0);
     const double* v_eff_tmp_const = pot->get_effective_v(0);
@@ -445,7 +465,7 @@ TEST_F(PotentialNewTest, GetEffectiveVofkmatrix)
     // construct potential
     elecstate::tmp_xc_func_type = 3;
     rhopw->nrxx = 100;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
     ModuleBase::matrix vofk_eff_tmp = pot->get_effective_vofk();
     const ModuleBase::matrix vofk_eff_tmp_const = pot->get_effective_vofk();
@@ -467,7 +487,7 @@ TEST_F(PotentialNewTest, GetEffectiveVofkarray)
 {
     // construct potential
     rhopw->nrxx = 100;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
     double* vofk_eff_tmp = pot->get_effective_vofk(0);
     const double* vofk_eff_tmp_const = pot->get_effective_vofk(0);
@@ -494,7 +514,7 @@ TEST_F(PotentialNewTest, GetEffectiveVofkarrayNullptr)
 TEST_F(PotentialNewTest, GetFixedV)
 {
     rhopw->nrxx = 100;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     EXPECT_TRUE(pot->fixed_mode);
     EXPECT_TRUE(pot->dynamic_mode);
     EXPECT_EQ(pot->v_effective_fixed.size(), 100);
@@ -513,7 +533,7 @@ TEST_F(PotentialNewTest, GetVeffSmooth)
     // construct potential
     rhopw->nrxx = 100;
     elecstate::tmp_xc_func_type = 3;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
     ModuleBase::matrix veff_smooth_tmp = pot->get_veff_smooth();
     const ModuleBase::matrix veff_smooth_const_tmp = pot->get_veff_smooth();
@@ -535,7 +555,7 @@ TEST_F(PotentialNewTest, GetVofkSmooth)
 {
     // construct potential
     rhopw->nrxx = 100;
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
     //
     ModuleBase::matrix vofk_smooth_tmp = pot->get_veff_smooth();
     const ModuleBase::matrix vofk_smooth_const_tmp = pot->get_veff_smooth();
@@ -568,7 +588,7 @@ TEST_F(PotentialNewTest, InterpolateVrsDoubleGrids)
     static_cast<ModulePW::PW_Basis_Sup*>(rhodpw)->setuptransform(rhopw);
     rhodpw->collect_local_pw();
 
-    pot = new elecstate::Potential(rhodpw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhodpw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
 
     for (int ir = 0; ir < pot->v_effective.nr; ir++)
     {
@@ -613,7 +633,7 @@ TEST_F(PotentialNewTest, InterpolateVrsWarningQuit)
     rhodpw->collect_local_pw();
     rhodpw->gamma_only = true;
 
-    pot = new elecstate::Potential(rhodpw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhodpw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
 
     EXPECT_EXIT(pot->interpolate_vrs(), ::testing::ExitedWithCode(1), "");
 }
@@ -628,7 +648,7 @@ TEST_F(PotentialNewTest, InterpolateVrsSingleGrids)
     rhopw->setuptransform();
     rhopw->collect_local_pw();
 
-    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc);
+    pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc);
 
     for (int ir = 0; ir < pot->v_effective.nr; ir++)
     {
diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp
index 94daacd37d..d2c820a45a 100644
--- a/source/module_esolver/esolver_fp.cpp
+++ b/source/module_esolver/esolver_fp.cpp
@@ -234,7 +234,8 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep)
                 this->pw_rhod,
                 this->pelec->charge,
                 &(ucell),
-                this->pelec->pot->get_fixed_v());
+                this->pelec->pot->get_fixed_v(),
+                this->solvent);
         }
 
         // 5) write ELF
diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h
index fc66cb42a4..1f2e6cf1ec 100644
--- a/source/module_esolver/esolver_fp.h
+++ b/source/module_esolver/esolver_fp.h
@@ -6,6 +6,7 @@
 #include "module_cell/module_symmetry/symmetry.h"
 #include "module_elecstate/elecstate.h"
 #include "module_elecstate/module_charge/charge_extra.h"
+#include "module_hamilt_general/module_surchem/surchem.h"
 #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h"
 
 #include <fstream>
@@ -67,6 +68,9 @@ namespace ModuleESolver
 
         //! Charge extrapolation
         Charge_Extra CE;
+
+        // solvent model
+        surchem solvent;
     };
 }
 
diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
index 65c4c774b1..712617948a 100644
--- a/source/module_esolver/esolver_ks_lcao.cpp
+++ b/source/module_esolver/esolver_ks_lcao.cpp
@@ -214,6 +214,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
                                                     &ucell,
                                                     &(this->ppcell.vloc),
                                                     &(this->sf),
+                                                    &(this->solvent),
                                                     &(this->pelec->f_en.etxc),
                                                     &(this->pelec->f_en.vtxc));
     }
@@ -320,6 +321,7 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
                        this->sf,
                        this->kv,
                        this->pw_rho,
+                       this->solvent,
 #ifdef __EXX
                        *this->exx_lri_double,
                        *this->exx_lri_complex,
@@ -460,6 +462,7 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
                                     *this->psi,
                                     ucell,
                                     this->sf,
+                                    this->solvent,
                                     *this->pw_rho,
                                     *this->pw_rhod,
                                     this->ppcell.vloc,
@@ -487,6 +490,7 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
                                             *this->psi,
                                             ucell,
                                             this->sf,
+                                            this->solvent,
                                             *this->pw_rho,
                                             *this->pw_rhod,
                                             this->ppcell.vloc,
diff --git a/source/module_esolver/esolver_ks_lcaopw.cpp b/source/module_esolver/esolver_ks_lcaopw.cpp
index e79346cd6d..6abf05e9d7 100644
--- a/source/module_esolver/esolver_ks_lcaopw.cpp
+++ b/source/module_esolver/esolver_ks_lcaopw.cpp
@@ -247,6 +247,7 @@ namespace ModuleESolver
                                 *this->kspw_psi,
                                 ucell,
                                 this->sf,
+                                this->solvent,
                                 *this->pw_wfc,
                                 *this->pw_rho,
                                 *this->pw_rhod,
diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
index 5248bcce0b..2171948312 100644
--- a/source/module_esolver/esolver_ks_pw.cpp
+++ b/source/module_esolver/esolver_ks_pw.cpp
@@ -181,6 +181,7 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
                                                     &ucell,
                                                     &this->ppcell.vloc,
                                                     &(this->sf),
+                                                    &(this->solvent),
                                                     &(this->pelec->f_en.etxc),
                                                     &(this->pelec->f_en.vtxc));
     }
@@ -778,6 +779,7 @@ void ESolver_KS_PW<T, Device>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
                  this->pw_rhod,
                  &ucell.symm,
                  &this->sf,
+                 this->solvent,
                  &this->ppcell,
                  &this->ppcell,
                  &this->kv,
diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp
index 7abb63553b..ac7d01b558 100644
--- a/source/module_esolver/esolver_of.cpp
+++ b/source/module_esolver/esolver_of.cpp
@@ -548,7 +548,7 @@ double ESolver_OF::cal_energy()
 void ESolver_OF::cal_force(UnitCell& ucell, ModuleBase::matrix& force)
 {
     Forces<double> ff(ucell.nat);
-    ff.cal_force(ucell,force, *pelec, this->pw_rho, &ucell.symm, &sf, &this->locpp);
+    ff.cal_force(ucell, force, *pelec, this->pw_rho, &ucell.symm, &sf, this->solvent, &this->locpp);
 }
 
 /**
diff --git a/source/module_esolver/esolver_of_tool.cpp b/source/module_esolver/esolver_of_tool.cpp
index 30ff2f5175..750598ab2e 100644
--- a/source/module_esolver/esolver_of_tool.cpp
+++ b/source/module_esolver/esolver_of_tool.cpp
@@ -30,6 +30,7 @@ void ESolver_OF::init_elecstate(UnitCell& ucell)
                                                 &ucell,
                                                 &(this->locpp.vloc),
                                                 &(this->sf),
+                                                &(this->solvent),
                                                 &(this->pelec->f_en.etxc),
                                                 &(this->pelec->f_en.vtxc));
     // There is no Operator in ESolver_OF, register Potentials here!
diff --git a/source/module_hamilt_general/module_surchem/cal_totn.cpp b/source/module_hamilt_general/module_surchem/cal_totn.cpp
index 05fe619761..3650a67dfa 100644
--- a/source/module_hamilt_general/module_surchem/cal_totn.cpp
+++ b/source/module_hamilt_general/module_surchem/cal_totn.cpp
@@ -33,7 +33,7 @@ void surchem::cal_totn(const UnitCell& cell,
     return;
 }
 
-void surchem::induced_charge(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, double* induced_rho)
+void surchem::induced_charge(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, double* induced_rho) const
 {
     std::complex<double> *delta_phig = new complex<double>[rho_basis->npw];
     std::complex<double> *induced_rhog = new complex<double>[rho_basis->npw];
diff --git a/source/module_hamilt_general/module_surchem/cal_vcav.cpp b/source/module_hamilt_general/module_surchem/cal_vcav.cpp
index 458034e558..3f7e768f4b 100644
--- a/source/module_hamilt_general/module_surchem/cal_vcav.cpp
+++ b/source/module_hamilt_general/module_surchem/cal_vcav.cpp
@@ -1,5 +1,6 @@
 #include "module_base/timer.h"
 #include "module_hamilt_general/module_xc/xc_functional.h"
+#include "module_parameter/parameter.h"
 #include "surchem.h"
 
 void lapl_rho(const double& tpiba2,
@@ -25,8 +26,9 @@ void lapl_rho(const double& tpiba2,
         // bring the gdr from G --> R
         rho_basis->recip2real(aux, aux);
         // remember to multily 2pi/a0, which belongs to G vectors.
-        for (int ir = 0; ir < rho_basis->nrxx; ir++)
+        for (int ir = 0; ir < rho_basis->nrxx; ir++) {
             lapn[ir] -= aux[ir].real() * tpiba2;
+}
 
     }
 
diff --git a/source/module_hamilt_general/module_surchem/sol_force.cpp b/source/module_hamilt_general/module_surchem/sol_force.cpp
index 1769f8886d..00c4da4404 100644
--- a/source/module_hamilt_general/module_surchem/sol_force.cpp
+++ b/source/module_hamilt_general/module_surchem/sol_force.cpp
@@ -2,10 +2,10 @@
 #include "module_base/timer.h"
 #include "module_parameter/parameter.h"
 
-void force_cor_one(const UnitCell& cell,
-                   const ModulePW::PW_Basis* rho_basis,
-                   const ModuleBase::matrix& vloc,
-                   ModuleBase::matrix& forcesol)
+void surchem::force_cor_one(const UnitCell& cell,
+                            const ModulePW::PW_Basis* rho_basis,
+                            const ModuleBase::matrix& vloc,
+                            ModuleBase::matrix& forcesol)
 {
    
    
@@ -15,10 +15,11 @@ void force_cor_one(const UnitCell& cell,
     std::complex<double> *delta_phi_g = new complex<double>[rho_basis->npw];
     //ModuleBase::GlobalFunc::ZEROS(delta_phi_g, rho_basis->npw);
 
-    rho_basis->real2recip(GlobalC::solvent_model.delta_phi, delta_phi_g);
-    //GlobalC::UFFT.ToReciSpace(GlobalC::solvent_model.delta_phi, delta_phi_g,rho_basis);
-    double Ael=0;double Ael1 = 0;
-    //ModuleBase::GlobalFunc::ZEROS(vg, ngmc);
+    rho_basis->real2recip(this->delta_phi, delta_phi_g);
+    // GlobalC::UFFT.ToReciSpace(this->delta_phi, delta_phi_g,rho_basis);
+    // double Ael = 0;
+    // double Ael1 = 0;
+    // ModuleBase::GlobalFunc::ZEROS(vg, ngmc);
     int iat = 0;
 
     for (int it = 0;it < cell.ntype;it++)
@@ -68,14 +69,14 @@ void force_cor_one(const UnitCell& cell,
 
 }
 
-void force_cor_two(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, ModuleBase::matrix& forcesol)
+void surchem::force_cor_two(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, ModuleBase::matrix& forcesol)
 {
    
     complex<double> *n_pseudo = new complex<double>[rho_basis->npw];
     ModuleBase::GlobalFunc::ZEROS(n_pseudo,rho_basis->npw);
 
-    //GlobalC::solvent_model.gauss_charge(cell, pwb, n_pseudo);
-    
+    // this->gauss_charge(cell, pwb, n_pseudo);
+
     double *Vcav_sum =  new double[rho_basis->nrxx];
     ModuleBase::GlobalFunc::ZEROS(Vcav_sum, rho_basis->nrxx);
     std::complex<double> *Vcav_g = new complex<double>[rho_basis->npw];
@@ -86,18 +87,18 @@ void force_cor_two(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, Mo
 	{
 		for (int ir=0; ir<rho_basis->nrxx; ir++)
 		{
-        	Vcav_sum[ir] += GlobalC::solvent_model.Vcav(is, ir);
-		}
-	}
+            Vcav_sum[ir] += this->Vcav(is, ir);
+        }
+    }
 
     rho_basis->real2recip(Vcav_sum, Vcav_g);
-    rho_basis->real2recip(GlobalC::solvent_model.epspot, Vel_g);
+    rho_basis->real2recip(this->epspot, Vel_g);
 
     int iat = 0;
-    double Ael1 = 0;
+    // double Ael1 = 0;
     for (int it = 0;it < cell.ntype;it++)
     {
-        double RCS = GlobalC::solvent_model.GetAtom.atom_RCS[cell.atoms[it].ncpp.psd];
+        double RCS = this->GetAtom.atom_RCS[cell.atoms[it].ncpp.psd];
         double sigma_rc_k = RCS / 2.5;
         for (int ia = 0;ia < cell.atoms[it].na;ia++)
         {
@@ -111,10 +112,10 @@ void force_cor_two(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, Mo
                 gg = gg * cell.tpiba2;
                 complex<double> phase = exp( ModuleBase::NEG_IMAG_UNIT *ModuleBase::TWO_PI * ( rho_basis->gcar[ig] * cell.atoms[it].tau[ia]));
 
-                n_pseudo[ig].real((GlobalC::solvent_model.GetAtom.atom_Z[cell.atoms[it].ncpp.psd] - cell.atoms[it].ncpp.zv) * phase.real()
-                             * exp(-0.5 * gg * (sigma_rc_k * sigma_rc_k)));
-                n_pseudo[ig].imag((GlobalC::solvent_model.GetAtom.atom_Z[cell.atoms[it].ncpp.psd] - cell.atoms[it].ncpp.zv) * phase.imag()
-                             * exp(-0.5 * gg * (sigma_rc_k * sigma_rc_k)));
+                n_pseudo[ig].real((this->GetAtom.atom_Z[cell.atoms[it].ncpp.psd] - cell.atoms[it].ncpp.zv)
+                                  * phase.real() * exp(-0.5 * gg * (sigma_rc_k * sigma_rc_k)));
+                n_pseudo[ig].imag((this->GetAtom.atom_Z[cell.atoms[it].ncpp.psd] - cell.atoms[it].ncpp.zv)
+                                  * phase.imag() * exp(-0.5 * gg * (sigma_rc_k * sigma_rc_k)));
             }
             
             for (int ig = 0; ig < rho_basis->npw; ig++)
diff --git a/source/module_hamilt_general/module_surchem/surchem.cpp b/source/module_hamilt_general/module_surchem/surchem.cpp
index 97acd8de54..9fb3541d9a 100644
--- a/source/module_hamilt_general/module_surchem/surchem.cpp
+++ b/source/module_hamilt_general/module_surchem/surchem.cpp
@@ -1,9 +1,7 @@
 #include "surchem.h"
 
-namespace GlobalC
-{
-surchem solvent_model;
-}
+double surchem::Acav = 0;
+double surchem::Ael = 0;
 
 surchem::surchem()
 {
diff --git a/source/module_hamilt_general/module_surchem/surchem.h b/source/module_hamilt_general/module_surchem/surchem.h
index 60d6b7e62f..091a14ddec 100644
--- a/source/module_hamilt_general/module_surchem/surchem.h
+++ b/source/module_hamilt_general/module_surchem/surchem.h
@@ -8,7 +8,6 @@
 #include "module_base/parallel_reduce.h"
 #include "module_basis/module_pw/pw_basis.h"
 #include "module_cell/unitcell.h"
-#include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h"
 
 class surchem
@@ -24,8 +23,8 @@ class surchem
     ModuleBase::matrix Vel;
     double qs;
 
-    double Acav;
-    double Ael;
+    static double Acav;
+    static double Ael;
 
     // get atom info
     atom_in GetAtom;
@@ -114,16 +113,18 @@ class surchem
                        const ModuleBase::matrix& vloc,
                        ModuleBase::matrix& forcesol);
 
+    void force_cor_one(const UnitCell& cell,
+                       const ModulePW::PW_Basis* rho_basis,
+                       const ModuleBase::matrix& vloc,
+                       ModuleBase::matrix& forcesol);
+
+    void force_cor_two(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, ModuleBase::matrix& forcesol);
+
     void get_totn_reci(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, complex<double>* totn_reci);
 
-    void induced_charge(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, double* induced_rho);
+    void induced_charge(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, double* induced_rho) const;
 
   private:
 };
 
-namespace GlobalC
-{
-extern surchem solvent_model;
-}
-
 #endif
diff --git a/source/module_hamilt_general/module_surchem/test/setcell.h b/source/module_hamilt_general/module_surchem/test/setcell.h
index c6a4545438..5d3ccd8ede 100644
--- a/source/module_hamilt_general/module_surchem/test/setcell.h
+++ b/source/module_hamilt_general/module_surchem/test/setcell.h
@@ -7,6 +7,7 @@
 #include "module_cell/module_neighbor/sltk_atom_arrange.h"
 #include "module_cell/module_neighbor/sltk_grid_driver.h"
 #include "module_cell/unitcell.h"
+#include "module_hamilt_pw/hamilt_pwdft/parallel_grid.h"
 #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h"
 
 namespace GlobalC
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
index edd8def7a8..9620a5f33c 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
@@ -50,6 +50,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                                           const Structure_Factor& sf,
                                           const K_Vectors& kv,
                                           ModulePW::PW_Basis* rhopw,
+                                          surchem& solvent,
 #ifdef __EXX
                                           Exx_LRI<double>& exx_lri_double,
                                           Exx_LRI<std::complex<double>>& exx_lri_complex,
@@ -279,7 +280,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
     if (PARAM.inp.imp_sol && isforce)
     {
         fsol.create(nat, 3);
-        GlobalC::solvent_model.cal_force_sol(ucell, rhopw, nlpp.vloc, fsol);
+        solvent.cal_force_sol(ucell, rhopw, nlpp.vloc, fsol);
     }
 
     //! atomic forces from DFT+U (Quxin version)
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h
index af6446ba7b..c995d57b09 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h
@@ -48,6 +48,7 @@ class Force_Stress_LCAO
                         const Structure_Factor& sf,
                         const K_Vectors& kv,
                         ModulePW::PW_Basis* rhopw,
+                        surchem& solvent,
 #ifdef __EXX
                         Exx_LRI<double>& exx_lri_double,
                         Exx_LRI<std::complex<double>>& exx_lri_complex,
diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp
index dd24ce82c7..372836eaa9 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp
@@ -30,6 +30,7 @@ void Forces<FPTYPE, Device>::cal_force(const UnitCell& ucell,
                                        ModulePW::PW_Basis* rho_basis,
                                        ModuleSymmetry::Symmetry* p_symm,
                                        Structure_Factor* p_sf,
+                                       surchem& solvent,
                                        const pseudopot_cell_vl* locpp,
                                        const pseudopot_cell_vnl* p_nlpp,
                                        K_Vectors* pkv,
@@ -231,7 +232,7 @@ void Forces<FPTYPE, Device>::cal_force(const UnitCell& ucell,
     if (PARAM.inp.imp_sol)
     {
         forcesol.create(this->nat, 3);
-        GlobalC::solvent_model.cal_force_sol(ucell, rho_basis, locpp->vloc, forcesol);
+        solvent.cal_force_sol(ucell, rho_basis, locpp->vloc, forcesol);
         if (PARAM.inp.test_force)
         {
             ModuleIO::print_force(GlobalV::ofs_running, ucell, "IMP_SOL      FORCE (Ry/Bohr)", forcesol);
diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.h b/source/module_hamilt_pw/hamilt_pwdft/forces.h
index 90b419199d..a01aaa7f02 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/forces.h
+++ b/source/module_hamilt_pw/hamilt_pwdft/forces.h
@@ -39,6 +39,7 @@ class Forces
                    ModulePW::PW_Basis* rho_basis,
                    ModuleSymmetry::Symmetry* p_symm,
                    Structure_Factor* p_sf,
+                   surchem& solvent,
                    const pseudopot_cell_vl* locpp,
                    const pseudopot_cell_vnl* nlpp = nullptr,
                    K_Vectors* pkv = nullptr,
diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp
index 12162d88e4..3346724deb 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp
@@ -9,12 +9,13 @@
 #include "module_base/mathzone.h"
 #include "module_base/timer.h"
 #include "module_base/tool_threading.h"
+#include "module_elecstate/cal_ux.h"
 #include "module_elecstate/potentials/efield.h"
 #include "module_elecstate/potentials/gatefield.h"
 #include "module_hamilt_general/module_ewald/H_Ewald_pw.h"
 #include "module_hamilt_general/module_surchem/surchem.h"
 #include "module_hamilt_general/module_vdw/vdw.h"
-#include "module_elecstate/cal_ux.h"
+#include "module_hamilt_general/module_xc/xc_functional.h"
 
 #ifdef _OPENMP
 #include <omp.h>
diff --git a/source/module_io/write_eband_terms.hpp b/source/module_io/write_eband_terms.hpp
index 41f0a42112..eff46a20b9 100644
--- a/source/module_io/write_eband_terms.hpp
+++ b/source/module_io/write_eband_terms.hpp
@@ -6,30 +6,32 @@
 #include "module_basis/module_nao/two_center_bundle.h"
 namespace ModuleIO
 {
-    template <typename TK, typename TR>
-    void write_eband_terms(const int nspin,
-        const int nbasis,
-        const int drank,
-        const Parallel_Orbitals* pv,
-        const psi::Psi<TK>& psi,
-        const UnitCell& ucell,
-        Structure_Factor& sf,
-        const ModulePW::PW_Basis& rho_basis,
-        const ModulePW::PW_Basis& rhod_basis,
-        const ModuleBase::matrix& vloc,
-        const Charge& chg,
-        Gint_Gamma& gint_gamma, // mohan add 2024-04-01
-        Gint_k& gint_k,         // mohan add 2024-04-01
-        const K_Vectors& kv,
-        const ModuleBase::matrix& wg,
-        Grid_Driver& gd,
-        const std::vector<double>& orb_cutoff,
-        const TwoCenterBundle& two_center_bundle
+template <typename TK, typename TR>
+void write_eband_terms(const int nspin,
+                       const int nbasis,
+                       const int drank,
+                       const Parallel_Orbitals* pv,
+                       const psi::Psi<TK>& psi,
+                       const UnitCell& ucell,
+                       Structure_Factor& sf,
+                       surchem& solvent,
+                       const ModulePW::PW_Basis& rho_basis,
+                       const ModulePW::PW_Basis& rhod_basis,
+                       const ModuleBase::matrix& vloc,
+                       const Charge& chg,
+                       Gint_Gamma& gint_gamma, // mohan add 2024-04-01
+                       Gint_k& gint_k,         // mohan add 2024-04-01
+                       const K_Vectors& kv,
+                       const ModuleBase::matrix& wg,
+                       Grid_Driver& gd,
+                       const std::vector<double>& orb_cutoff,
+                       const TwoCenterBundle& two_center_bundle
 #ifdef __EXX
-        , std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd = nullptr
-        , std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc = nullptr
+                       ,
+                       std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd = nullptr,
+                       std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc = nullptr
 #endif
-    )
+)
     {
         // 0. prepare
         const int& nbands = wg.nc;
@@ -80,7 +82,7 @@ namespace ModuleIO
         // 2. pp: local
         if (PARAM.inp.vl_in_h)
         {
-            elecstate::Potential pot_local(&rhod_basis, &rho_basis, &ucell, &vloc, &sf, &etxc, &vtxc);
+            elecstate::Potential pot_local(&rhod_basis, &rho_basis, &ucell, &vloc, &sf, &solvent, &etxc, &vtxc);
             pot_local.pot_register({ "local" });
             pot_local.update_from_charge(&chg, &ucell);
             hamilt::HS_Matrix_K<TK> v_pp_local_k_ao(pv, 1);
@@ -128,7 +130,7 @@ namespace ModuleIO
         // 4. hartree
         if (PARAM.inp.vh_in_h)
         {
-            elecstate::Potential pot_hartree(&rhod_basis, &rho_basis, &ucell, &vloc, &sf, &etxc, &vtxc);
+            elecstate::Potential pot_hartree(&rhod_basis, &rho_basis, &ucell, &vloc, &sf, &solvent, &etxc, &vtxc);
             pot_hartree.pot_register({ "hartree" });
             pot_hartree.update_from_charge(&chg, &ucell);
             std::vector<hamilt::HContainer<TR>> v_hartree_R_ao(nspin0, hamilt::HContainer<TR>(pv));
@@ -164,9 +166,28 @@ namespace ModuleIO
         // 5. xc (including exx)
         if (!PARAM.inp.out_mat_xc)  // avoid duplicate output
         {
-            write_Vxc<TK, TR>(nspin, nbasis, drank, pv, psi, ucell, sf, rho_basis, rhod_basis, vloc, chg, gint_gamma, gint_k, kv, orb_cutoff, wg, gd
+            write_Vxc<TK, TR>(nspin,
+                              nbasis,
+                              drank,
+                              pv,
+                              psi,
+                              ucell,
+                              sf,
+                              solvent,
+                              rho_basis,
+                              rhod_basis,
+                              vloc,
+                              chg,
+                              gint_gamma,
+                              gint_k,
+                              kv,
+                              orb_cutoff,
+                              wg,
+                              gd
 #ifdef __EXX
-                , Hexxd, Hexxc
+                              ,
+                              Hexxd,
+                              Hexxc
 #endif
             );
         }
diff --git a/source/module_io/write_elecstat_pot.cpp b/source/module_io/write_elecstat_pot.cpp
index f4561bb606..6f9fb42a57 100644
--- a/source/module_io/write_elecstat_pot.cpp
+++ b/source/module_io/write_elecstat_pot.cpp
@@ -21,7 +21,8 @@ void write_elecstat_pot(
     ModulePW::PW_Basis* rho_basis,
     const Charge* const chr,
     const UnitCell* ucell,
-    const double* v_eff)
+    const double* v_eff,
+    const surchem& solvent)
 {
     ModuleBase::TITLE("ModuleIO", "write_elecstat_pot");
     ModuleBase::timer::tick("ModuleIO", "write_elecstat_pot");
@@ -50,7 +51,7 @@ void write_elecstat_pot(
                                                  const_cast<ModulePW::PW_Basis*>(rho_basis),
                                                  nspin,
                                                  chr->rho,
-                                                 GlobalC::solvent_model);
+                                                 solvent);
     }
 
     //==========================================
@@ -67,7 +68,7 @@ void write_elecstat_pot(
         }
         if(imp_sol == true)
         {
-            v_elecstat[ir] += GlobalC::solvent_model.delta_phi[ir];
+            v_elecstat[ir] += solvent.delta_phi[ir];
         }
     }
 
diff --git a/source/module_io/write_elecstat_pot.h b/source/module_io/write_elecstat_pot.h
index 14640d9e43..19541ee22f 100644
--- a/source/module_io/write_elecstat_pot.h
+++ b/source/module_io/write_elecstat_pot.h
@@ -1,9 +1,11 @@
 #ifndef POTENTIAL_IO_H
 #define POTENTIAL_IO_H
-#include <string>
-#include "module_cell/unitcell.h"
 #include "module_basis/module_pw/pw_basis.h"
+#include "module_cell/unitcell.h"
 #include "module_elecstate/module_charge/charge.h"
+#include "module_hamilt_general/module_surchem/surchem.h"
+
+#include <string>
 
 namespace ModuleIO
 {
@@ -27,7 +29,8 @@ void write_elecstat_pot(
     ModulePW::PW_Basis* rho_basis,
     const Charge* const chr,
     const UnitCell* ucell_,
-    const double* v_effective_fixed);
+    const double* v_effective_fixed,
+    const surchem& solvent);
 
 } // namespace ModuleIO
 
diff --git a/source/module_io/write_vxc.hpp b/source/module_io/write_vxc.hpp
index fcd03e07a3..fa9da0ec62 100644
--- a/source/module_io/write_vxc.hpp
+++ b/source/module_io/write_vxc.hpp
@@ -182,25 +182,27 @@ inline void write_orb_energy(const K_Vectors& kv,
 /// including terms: local/semi-local XC, EXX, DFTU
 template <typename TK, typename TR>
 void write_Vxc(const int nspin,
-    const int nbasis,
-    const int drank,
-    const Parallel_Orbitals* pv,
-    const psi::Psi<TK>& psi,
-    const UnitCell& ucell,
-    Structure_Factor& sf,
-    const ModulePW::PW_Basis& rho_basis,
-    const ModulePW::PW_Basis& rhod_basis,
-    const ModuleBase::matrix& vloc,
-    const Charge& chg,
-    Gint_Gamma& gint_gamma, // mohan add 2024-04-01
-    Gint_k& gint_k,         // mohan add 2024-04-01
-    const K_Vectors& kv,
-    const std::vector<double>& orb_cutoff,
-    const ModuleBase::matrix& wg,
-    Grid_Driver& gd
+               const int nbasis,
+               const int drank,
+               const Parallel_Orbitals* pv,
+               const psi::Psi<TK>& psi,
+               const UnitCell& ucell,
+               Structure_Factor& sf,
+               surchem& solvent,
+               const ModulePW::PW_Basis& rho_basis,
+               const ModulePW::PW_Basis& rhod_basis,
+               const ModuleBase::matrix& vloc,
+               const Charge& chg,
+               Gint_Gamma& gint_gamma, // mohan add 2024-04-01
+               Gint_k& gint_k,         // mohan add 2024-04-01
+               const K_Vectors& kv,
+               const std::vector<double>& orb_cutoff,
+               const ModuleBase::matrix& wg,
+               Grid_Driver& gd
 #ifdef __EXX
-    , std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd = nullptr
-    , std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc = nullptr
+               ,
+               std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd = nullptr,
+               std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc = nullptr
 #endif
 )
 {
@@ -212,7 +214,8 @@ void write_Vxc(const int nspin,
     double vtxc = 0.0;
     // elecstate::PotXC* potxc(&rho_basis, &etxc, vtxc, nullptr);
     // potxc.cal_v_eff(&chg, &ucell, vr_xc);
-    elecstate::Potential* potxc = new elecstate::Potential(&rhod_basis, &rho_basis, &ucell, &vloc, &sf, &etxc, &vtxc);
+    elecstate::Potential* potxc
+        = new elecstate::Potential(&rhod_basis, &rho_basis, &ucell, &vloc, &sf, &solvent, &etxc, &vtxc);
     std::vector<std::string> compnents_list = {"xc"};
     potxc->pot_register(compnents_list);
     potxc->update_from_charge(&chg, &ucell);
diff --git a/source/module_io/write_vxc_lip.hpp b/source/module_io/write_vxc_lip.hpp
index f3c30bd5e2..205fdbb057 100644
--- a/source/module_io/write_vxc_lip.hpp
+++ b/source/module_io/write_vxc_lip.hpp
@@ -84,21 +84,23 @@ namespace ModuleIO
     /// including terms: local/semi-local XC and EXX
     template <typename FPTYPE>
     void write_Vxc(int nspin,
-        int naos,
-        int drank,
-        const psi::Psi<std::complex<FPTYPE>>& psi_pw,
-        // const psi::Psi<T>& psi_lcao,
-        const UnitCell& ucell,
-        Structure_Factor& sf,
-        const ModulePW::PW_Basis_K& wfc_basis,
-        const ModulePW::PW_Basis& rho_basis,
-        const ModulePW::PW_Basis& rhod_basis,
-        const ModuleBase::matrix& vloc,
-        const Charge& chg,
-        const K_Vectors& kv,
-        const ModuleBase::matrix& wg
+                   int naos,
+                   int drank,
+                   const psi::Psi<std::complex<FPTYPE>>& psi_pw,
+                   // const psi::Psi<T>& psi_lcao,
+                   const UnitCell& ucell,
+                   Structure_Factor& sf,
+                   surchem& solvent,
+                   const ModulePW::PW_Basis_K& wfc_basis,
+                   const ModulePW::PW_Basis& rho_basis,
+                   const ModulePW::PW_Basis& rhod_basis,
+                   const ModuleBase::matrix& vloc,
+                   const Charge& chg,
+                   const K_Vectors& kv,
+                   const ModuleBase::matrix& wg
 #ifdef __EXX
-        , const Exx_Lip<std::complex<FPTYPE>>& exx_lip
+                   ,
+                   const Exx_Lip<std::complex<FPTYPE>>& exx_lip
 #endif
     )
     {
@@ -111,7 +113,8 @@ namespace ModuleIO
         double vtxc = 0.0;
         // elecstate::PotXC* potxc(&rho_basis, &etxc, vtxc, nullptr);
         // potxc.cal_v_eff(&chg, &ucell, vr_xc);
-        elecstate::Potential* potxc = new elecstate::Potential(&rhod_basis, &rho_basis, &ucell, &vloc, &sf, &etxc, &vtxc);
+        elecstate::Potential* potxc
+            = new elecstate::Potential(&rhod_basis, &rho_basis, &ucell, &vloc, &sf, &solvent, &etxc, &vtxc);
         std::vector<std::string> compnents_list = { "xc" };
 
         potxc->pot_register(compnents_list);