Refactor: remove GlobalC::ucell in esolver

source/driver.cpp

@@ -183,7 +183,7 @@ void Driver::atomic_world()
     //--------------------------------------------------
 
     // where the actual stuff is done
-    this->driver_run();
+    this->driver_run(GlobalC::ucell);
 
     ModuleBase::timer::finish(GlobalV::ofs_running);
     ModuleBase::Memory::print_all(GlobalV::ofs_running);

source/driver.h

@@ -1,6 +1,8 @@
 #ifndef DRIVER_H
 #define DRIVER_H
 
+#include "module_cell/unitcell.h"
+
 class Driver
 {
   public:
@@ -34,7 +36,7 @@ class Driver
     void atomic_world();
 
     // the actual calculations
-    void driver_run();
+    void driver_run(UnitCell& ucell);
 };
 
 #endif

source/driver_run.cpp

@@ -24,7 +24,8 @@
  * the configuration-changing subroutine takes force and stress and updates the
  * configuration
  */
-void Driver::driver_run() {
+void Driver::driver_run(UnitCell& ucell)
+{
     ModuleBase::TITLE("Driver", "driver_line");
     ModuleBase::timer::tick("Driver", "driver_line");
 
@@ -39,37 +40,35 @@ void Driver::driver_run() {
 #endif
 
     // the life of ucell should begin here, mohan 2024-05-12
-    // delete ucell as a GlobalC in near future
-    GlobalC::ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running);
-    Check_Atomic_Stru::check_atomic_stru(GlobalC::ucell,
-                                         PARAM.inp.min_dist_coef);
+    ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running);
+    Check_Atomic_Stru::check_atomic_stru(ucell, PARAM.inp.min_dist_coef);
 
     //! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`)
-    ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, GlobalC::ucell);
+    ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, ucell);
 
     //! 3: initialize Esolver and fill json-structure
-    p_esolver->before_all_runners(PARAM.inp, GlobalC::ucell);
+    p_esolver->before_all_runners(PARAM.inp, ucell);
 
     // this Json part should be moved to before_all_runners, mohan 2024-05-12
 #ifdef __RAPIDJSON
-    Json::gen_stru_wrapper(&GlobalC::ucell);
+    Json::gen_stru_wrapper(&ucell);
 #endif
 
     const std::string cal_type = PARAM.inp.calculation;
 
     //! 4: different types of calculations
     if (cal_type == "md")
     {
-        Run_MD::md_line(GlobalC::ucell, p_esolver, PARAM);
+        Run_MD::md_line(ucell, p_esolver, PARAM);
     }
     else if (cal_type == "scf" || cal_type == "relax" || cal_type == "cell-relax" || cal_type == "nscf")
     {
         Relax_Driver rl_driver;
-        rl_driver.relax_driver(p_esolver);
+        rl_driver.relax_driver(p_esolver, ucell);
     }
     else if (cal_type == "get_S")
     {
-        p_esolver->runner(0, GlobalC::ucell);
+        p_esolver->runner(0, ucell);
     }
     else
     {
@@ -79,11 +78,11 @@ void Driver::driver_run() {
         //! test_neighbour(LCAO),
         //! gen_bessel(PW), et al.
         const int istep = 0;
-        p_esolver->others(istep);
+        p_esolver->others(istep, ucell);
     }
 
     //! 5: clean up esolver
-    p_esolver->after_all_runners();
+    p_esolver->after_all_runners(ucell);
 
     ModuleESolver::clean_esolver(p_esolver);
 

source/module_base/opt_TN.hpp

@@ -1,9 +1,10 @@
 #ifndef OPT_TN_H
 #define OPT_TN_H
 
-#include <limits>
-
 #include "./opt_CG.h"
+#include "module_cell/unitcell.h"
+
+#include <limits>
 
 namespace ModuleBase
 {
@@ -53,8 +54,9 @@ class Opt_TN
     void refresh(int nx_new = 0)
     {
         this->iter_ = 0;
-        if (nx_new != 0)
+        if (nx_new != 0) {
             this->nx_ = nx_new;
+}
         this->cg_.refresh(nx_new);
     }
 
@@ -64,10 +66,11 @@ class Opt_TN
         double* pgradient, // df(x)/dx
         int& flag,       // record which truncated condition was triggered, 0 for cond.1, 1 for cond.2, and 2 for cond.3
         double* rdirect, // next optimization direction
-        T* t,            // point of class T, which contains the gradient function
-        void (T::*p_calGradient)(
-            double* ptemp_x,
-            double* rtemp_gradient) // a function point, which calculates the gradient at provided x
+        UnitCell& ucell,
+        T* t, // point of class T, which contains the gradient function
+        void (T::*p_calGradient)(double* ptemp_x,
+                                 double* rtemp_gradient,
+                                 UnitCell& ucell) // a function point, which calculates the gradient at provided x
     );
 
     int get_iter()
@@ -128,8 +131,9 @@ void Opt_TN::next_direct(double* px,
                          double* pgradient,
                          int& flag,
                          double* rdirect,
+                         UnitCell& ucell,
                          T* t,
-                         void (T::*p_calGradient)(double* px, double* rgradient))
+                         void (T::*p_calGradient)(double* px, double* rgradient, UnitCell& ucell))
 {
     // initialize arrays and parameters
     ModuleBase::GlobalFunc::ZEROS(rdirect, this->nx_); // very important
@@ -166,18 +170,21 @@ void Opt_TN::next_direct(double* px,
         // Hcgd = (df(temp_x)/dx - df(x)/x) / epsilon, where temp_x = x + step * cg_direct
         epsilon = this->get_epsilon(px, cg_direct);
         // epsilon = 1e-9;
-        for (int i = 0; i < this->nx_; ++i)
+        for (int i = 0; i < this->nx_; ++i) {
             temp_x[i] = px[i] + epsilon * cg_direct[i];
-        (t->*p_calGradient)(temp_x, temp_gradient);
-        for (int i = 0; i < this->nx_; ++i)
+}
+        (t->*p_calGradient)(temp_x, temp_gradient, ucell);
+        for (int i = 0; i < this->nx_; ++i) {
             temp_Hcgd[i] = (temp_gradient[i] - pgradient[i]) / epsilon;
+}
 
         // get CG step length and update rdirect
         cg_alpha = cg_.step_length(temp_Hcgd, cg_direct, cg_ifPD);
         if (cg_ifPD == -1) // Hessian is not positive definite, and cgiter = 1.
         {
-            for (int i = 0; i < this->nx_; ++i)
+            for (int i = 0; i < this->nx_; ++i) {
                 rdirect[i] += cg_alpha * cg_direct[i];
+}
             flag = -1;
             break;
         }
@@ -187,15 +194,17 @@ void Opt_TN::next_direct(double* px,
             break;
         }
 
-        for (int i = 0; i < this->nx_; ++i)
+        for (int i = 0; i < this->nx_; ++i) {
             rdirect[i] += cg_alpha * cg_direct[i];
+}
 
         // store residuals used in truncated conditions
         last_residual = curr_residual;
         curr_residual = cg_.get_residual();
         cg_iter = cg_.get_iter();
-        if (cg_iter == 1)
+        if (cg_iter == 1) {
             init_residual = curr_residual;
+}
 
         // check truncated conditions
         // if (curr_residual < 1e-12)

source/module_cell/unitcell.cpp

@@ -551,6 +551,9 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) {
     this->atoms = new Atom[this->ntype]; // atom species.
     this->set_atom_flag = true;
 
+    this->symm.epsilon = PARAM.inp.symmetry_prec;
+    this->symm.epsilon_input = PARAM.inp.symmetry_prec;
+
     bool ok = true;
     bool ok2 = true;
 

source/module_esolver/esolver.h

@@ -20,26 +20,26 @@ class ESolver
     }
 
     //! initialize the energy solver by using input parameters and cell modules
-    virtual void before_all_runners(const Input_para& inp, UnitCell& cell) = 0;
+    virtual void before_all_runners(const Input_para& inp, UnitCell& ucell) = 0;
 
     //! run energy solver
     virtual void runner(const int istep, UnitCell& cell) = 0;
 
     //! perform post processing calculations
-    virtual void after_all_runners(){};
+    virtual void after_all_runners(UnitCell& ucell){};
 
     //! deal with exx and other calculation than scf/md/relax/cell-relax:
     //! such as nscf, get_wf and get_pchg
-    virtual void others(const int istep){};
+    virtual void others(const int istep, UnitCell& ucell){};
 
     //! calculate total energy of a given system
     virtual double cal_energy() = 0;
 
     //! calcualte forces for the atoms in the given cell
-    virtual void cal_force(ModuleBase::matrix& force) = 0;
+    virtual void cal_force(ModuleBase::matrix& force, UnitCell& ucell) = 0;
 
     //! calcualte stress of given cell
-    virtual void cal_stress(ModuleBase::matrix& stress) = 0;
+    virtual void cal_stress(ModuleBase::matrix& stress, UnitCell& ucell) = 0;
 
     bool conv_esolver = true; // whether esolver is converged
 

source/module_esolver/esolver_dp.cpp

@@ -127,13 +127,13 @@ double ESolver_DP::cal_energy()
     return dp_potential;
 }
 
-void ESolver_DP::cal_force(ModuleBase::matrix& force)
+void ESolver_DP::cal_force(ModuleBase::matrix& force, UnitCell& ucell)
 {
     force = dp_force;
     ModuleIO::print_force(GlobalV::ofs_running, *ucell_, "TOTAL-FORCE (eV/Angstrom)", force, false);
 }
 
-void ESolver_DP::cal_stress(ModuleBase::matrix& stress)
+void ESolver_DP::cal_stress(ModuleBase::matrix& stress, UnitCell& ucell)
 {
     stress = dp_virial;
 
@@ -148,7 +148,7 @@ void ESolver_DP::cal_stress(ModuleBase::matrix& stress)
     ModuleIO::print_stress("TOTAL-STRESS", stress, true, false);
 }
 
-void ESolver_DP::after_all_runners()
+void ESolver_DP::after_all_runners(UnitCell& ucell)
 {
     GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl;
     GlobalV::ofs_running << std::setprecision(16);

source/module_esolver/esolver_dp.h

@@ -36,7 +36,7 @@ class ESolver_DP : public ESolver
      * @param inp input parameters
      * @param cell unitcell information
      */
-    void before_all_runners(const Input_para& inp, UnitCell& cell) override;
+    void before_all_runners(const Input_para& inp, UnitCell& ucell) override;
 
     /**
      * @brief Run the DP solver for a given ion/md step and unit cell
@@ -59,21 +59,21 @@ class ESolver_DP : public ESolver
      *
      * @param force the computed atomic forces
      */
-    void cal_force(ModuleBase::matrix& force) override;
+    void cal_force(ModuleBase::matrix& force, UnitCell& ucell) override;
 
     /**
      * @brief get the computed lattice virials
      *
      * @param stress the computed lattice virials
      */
-    void cal_stress(ModuleBase::matrix& stress) override;
+    void cal_stress(ModuleBase::matrix& stress, UnitCell& ucell) override;
 
     /**
      * @brief Prints the final total energy of the DP model to the output file
      *
      * This function prints the final total energy of the DP model in eV to the output file along with some formatting.
      */
-    void after_all_runners() override;
+    void after_all_runners(UnitCell& ucell) override;
 
   private:
     /**