From d9db0d41cfd55c7ce84a56a164e01d9e607ffc94 Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Wed, 22 Nov 2023 14:33:52 +0800
Subject: [PATCH 1/4] put Symmetry into UnitCell

---
 source/module_cell/klist.cpp                  |   6 +-
 source/module_cell/klist.h                    |   2 +-
 .../module_symmetry/run_symmetry.cpp          |   2 +-
 .../module_cell/module_symmetry/symmetry.cpp  | 117 ++++++++----------
 source/module_cell/module_symmetry/symmetry.h |  33 ++---
 .../module_symmetry/symmetry_basic.cpp        |  10 --
 .../module_symmetry/symmetry_basic.h          |   4 +-
 .../module_symmetry/test/symmetry_test.h      |   2 +-
 .../test/symmetry_test_analysis.cpp           |   4 +-
 .../test/symmetry_test_symtrz.cpp             |   6 +-
 source/module_cell/test/klist_test.cpp        |   4 +-
 source/module_cell/test/klist_test_para.cpp   |   6 +-
 source/module_cell/unitcell.h                 |  76 +++++-------
 source/module_cell/unitcell_data.h            |  49 ++++++++
 source/module_esolver/esolver_fp.cpp          |   6 +-
 source/module_esolver/esolver_fp.h            |   1 -
 source/module_esolver/esolver_ks.cpp          |   4 +-
 source/module_esolver/esolver_ks_lcao.cpp     |   8 +-
 .../module_esolver/esolver_ks_lcao_elec.cpp   |   2 +-
 .../module_esolver/esolver_ks_lcao_tddft.cpp  |   2 +-
 source/module_esolver/esolver_ks_pw.cpp       |   8 +-
 source/module_esolver/esolver_of.cpp          |  14 +--
 source/module_esolver/esolver_sdft_pw.cpp     |   6 +-
 .../hamilt_lcaodft/FORCE_STRESS.cpp           |   4 +-
 .../hamilt_ofdft/of_stress_pw.cpp             |   2 +-
 .../hamilt_pwdft/stress_func_kin.cpp          |   2 +-
 .../hamilt_pwdft/stress_func_nl.cpp           |   2 +-
 .../hamilt_pwdft/stress_pw.cpp                |   2 +-
 .../hamilt_stodft/sto_stress_pw.cpp           |   6 +-
 .../module_io/test/for_testing_input_conv.h   |  12 --
 30 files changed, 201 insertions(+), 201 deletions(-)
 create mode 100644 source/module_cell/unitcell_data.h

diff --git a/source/module_cell/klist.cpp b/source/module_cell/klist.cpp
index 4931017828..a526b735b7 100644
--- a/source/module_cell/klist.cpp
+++ b/source/module_cell/klist.cpp
@@ -621,7 +621,7 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s
         ModuleBase::Vector3<double> gb01(b_optlat.e11, b_optlat.e12, b_optlat.e13);
         ModuleBase::Vector3<double> gb02(b_optlat.e21, b_optlat.e22, b_optlat.e23);
         ModuleBase::Vector3<double> gb03(b_optlat.e31, b_optlat.e32, b_optlat.e33);
-        symm.lattice_type(gb1, gb2, gb3, gb01, gb02, gb03, b_const, b0_const, bbrav, bbrav_name, ucell, false, nullptr);
+        symm.lattice_type(gb1, gb2, gb3, gb01, gb02, gb03, b_const, b0_const, bbrav, bbrav_name, ucell.atoms, false, nullptr);
         GlobalV::ofs_running<<"(for reciprocal lattice: )"<<std::endl;
         ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BRAVAIS TYPE", bbrav);
         ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BRAVAIS LATTICE NAME", bbrav_name);
@@ -646,7 +646,7 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s
             return;
         }
 
-        symm.lattice_type(gk1, gk2, gk3, gk01, gk02, gk03, bk_const, bk0_const, bkbrav, bkbrav_name, ucell, false, nullptr);
+        symm.lattice_type(gk1, gk2, gk3, gk01, gk02, gk03, bk_const, bk0_const, bkbrav, bkbrav_name, ucell.atoms, false, nullptr);
         GlobalV::ofs_running<<"(for k-lattice: )"<<std::endl;
         ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BRAVAIS TYPE", bkbrav);
         ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"BRAVAIS LATTICE NAME", bkbrav_name);
@@ -656,7 +656,7 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s
         ModuleBase::Matrix3 bsymop[48];
         int bnop = 0;
         // search again
-        symm.lattice_type(gb1, gb2, gb3, gb1, gb2, gb3, b_const, b0_const, bbrav, bbrav_name, ucell, false, nullptr);
+        symm.lattice_type(gb1, gb2, gb3, gb1, gb2, gb3, b_const, b0_const, bbrav, bbrav_name, ucell.atoms, false, nullptr);
         ModuleBase::Matrix3 b_optlat_new(gb1.x, gb1.y, gb1.z, gb2.x, gb2.y, gb2.z, gb3.x, gb3.y, gb3.z);
         symm.setgroup(bsymop, bnop, bbrav);
         symm.gmatrix_convert(bsymop, bsymop, bnop, b_optlat_new, ucell.G);
diff --git a/source/module_cell/klist.h b/source/module_cell/klist.h
index b7b0c3785d..a9e06f8614 100644
--- a/source/module_cell/klist.h
+++ b/source/module_cell/klist.h
@@ -6,7 +6,7 @@
 #include "module_base/global_function.h"
 #include "module_base/global_variable.h"
 #include "module_base/matrix3.h"
-#include "module_cell/module_symmetry/symmetry.h"
+#include "module_cell/unitcell.h"
 
 class K_Vectors
 {
diff --git a/source/module_cell/module_symmetry/run_symmetry.cpp b/source/module_cell/module_symmetry/run_symmetry.cpp
index f615bb21a2..f6604ac72e 100644
--- a/source/module_cell/module_symmetry/run_symmetry.cpp
+++ b/source/module_cell/module_symmetry/run_symmetry.cpp
@@ -40,7 +40,7 @@ void calculate()
 	"STRU", 
 	ofs_running);
 	std::cout << "set up cell classic done." << std::endl;
-	symm.analy_sys(ucell, out, ofs_running);
+    symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running);
 	ofs_running.close();
 //	ooo.set_orb_tables();
 
diff --git a/source/module_cell/module_symmetry/symmetry.cpp b/source/module_cell/module_symmetry/symmetry.cpp
index 2244e81787..b25146e080 100644
--- a/source/module_cell/module_symmetry/symmetry.cpp
+++ b/source/module_cell/module_symmetry/symmetry.cpp
@@ -8,23 +8,11 @@
 
 namespace ModuleSymmetry
 {
-Symmetry::Symmetry()
-{
-    this->epsilon = 1e-6;
-    this->tab = 12;
-    this->available = true;
-}
-
-Symmetry::~Symmetry()
-{
-}
-
-
 int Symmetry::symm_flag = 0;
 bool Symmetry::symm_autoclose = false;
 bool Symmetry::pricell_loop = true;
 
-void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
+void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, std::ofstream& ofs_running)
 {
     const double MAX_EPS = std::max(1e-3, epsilon_input * 1.001);
     const double MULT_EPS = 2.0;
@@ -47,9 +35,9 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
 
 
     // number of total atoms
-    this->nat = ucell.nat;
+    this->nat = st.nat;
     // number of atom species
-    this->ntype = ucell.ntype;
+    this->ntype = st.ntype;
     this->na = new int[ntype];
     this->istart = new int[ntype];
     this->index = new int [nat + 2];
@@ -64,9 +52,9 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
 	ModuleBase::GlobalFunc::ZEROS(newpos, 3*nat);
     ModuleBase::GlobalFunc::ZEROS(rotpos, 3*nat);
 
-    this->a1 = ucell.a1;
-    this->a2 = ucell.a2;
-    this->a3 = ucell.a3;
+    this->a1 = lat.a1;
+    this->a2 = lat.a2;
+    this->a3 = lat.a3;
 
 	ModuleBase::Matrix3 latvec1;
 	latvec1.e11 = a1.x; latvec1.e12 = a1.y; latvec1.e13 = a1.z;
@@ -83,7 +71,7 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
     this->itmin_start = 0;
     for (int it = 0; it < ntype; ++it)
     {
-		Atom* atom = &ucell.atoms[it];
+        Atom* atom = &atoms[it];
         this->na[it] = atom->na;
         if (it > 0) {
             istart[it] = istart[it-1] + na[it-1];
@@ -104,7 +92,7 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
         //a: optimized config
         // find the lattice type accordiing to lattice vectors.
         this->lattice_type(this->a1, this->a2, this->a3, this->s1, this->s2, this->s3,
-            this->cel_const, this->pre_const, this->real_brav, ilattname, ucell, true, this->newpos);
+            this->cel_const, this->pre_const, this->real_brav, ilattname, atoms, true, this->newpos);
 
         ofs_running << "(for optimal symmetric configuration:)" << std::endl;
         ModuleBase::GlobalFunc::OUT(ofs_running, "BRAVAIS TYPE", real_brav);
@@ -118,9 +106,8 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
         optlat.e21 = a2.x; optlat.e22 = a2.y; optlat.e23 = a2.z;
         optlat.e31 = a3.x; optlat.e32 = a3.y; optlat.e33 = a3.z;
 
-        this->pricell(this->newpos);         // pengfei Li 2018-05-14 
-        //for( iat =0 ; iat < ucell.nat ; iat++)   
-//         std::cout << " newpos_now = " << newpos[3*iat] << " " << newpos[3*iat+1] << " " << newpos[3*iat+2] << std::endl;
+        this->pricell(this->newpos, atoms);         // pengfei Li 2018-05-14 
+
         test_brav = true; // output the real ibrav and point group
         this->setgroup(this->symop, this->nop, this->real_brav);
         this->getgroup(nrot_out, nrotk_out, ofs_running);
@@ -235,16 +222,16 @@ void Symmetry::analy_sys(const UnitCell &ucell, std::ofstream &ofs_running)
     }
 
     //convert gmatrix to reciprocal space
-    this->gmatrix_convert_int(gmatrix, kgmatrix, nrotk, optlat, ucell.G);
+    this->gmatrix_convert_int(gmatrix, kgmatrix, nrotk, optlat, lat.G);
     
 // convert the symmetry operations from the basis of optimal symmetric configuration 
 // to the basis of input configuration
     this->gmatrix_convert_int(gmatrix, gmatrix, nrotk, optlat, latvec1);
     this->gtrans_convert(gtrans, gtrans, nrotk, optlat, latvec1);
 
-    this->set_atom_map(ucell);
+    this->set_atom_map(atoms);
 
-    if (GlobalV::NSPIN > 1) pricell_loop = this->magmom_same_check(ucell);
+    if (GlobalV::NSPIN > 1) pricell_loop = this->magmom_same_check(atoms);
 
 	delete[] newpos;
     delete[] na;
@@ -546,9 +533,9 @@ void Symmetry::lattice_type(
     double *cel_const,
     double *pre_const,
     int& real_brav,
-    std::string &bravname,
-    const UnitCell &ucell, 
-    bool convert_atoms, 
+    std::string& bravname,
+    const Atom* atoms,
+    bool convert_atoms,
     double* newpos)const
 {
     ModuleBase::TITLE("Symmetry","lattice_type");
@@ -657,13 +644,13 @@ void Symmetry::lattice_type(
             GlobalV::ofs_running <<" The lattice vectors have been changed (STRU_SIMPLE.cif)"<<std::endl;
             GlobalV::ofs_running <<std::endl;
             int at=0;
-            for(int it=0; it<ucell.ntype; ++it)
+            for (int it = 0; it < this->ntype; ++it)
             {
-                    for(int ia=0; ia<ucell.atoms[it].na; ++ia)
+                for (int ia = 0; ia < this->na[it]; ++ia)
                     {
-                            ModuleBase::Mathzone::Cartesian_to_Direct(ucell.atoms[it].tau[ia].x,
-                                            ucell.atoms[it].tau[ia].y,
-                                            ucell.atoms[it].tau[ia].z,
+                    ModuleBase::Mathzone::Cartesian_to_Direct(atoms[it].tau[ia].x,
+                        atoms[it].tau[ia].y,
+                        atoms[it].tau[ia].z,
                                             q1.x, q1.y, q1.z,
                                             q2.x, q2.y, q2.z,
                                             q3.x, q3.y, q3.z,
@@ -695,11 +682,11 @@ void Symmetry::lattice_type(
             ofs << std::endl;
             at =0;
             
-            for (int it=0; it<ucell.ntype; it++)
+            for (int it = 0; it < this->ntype; it++)
             {
-                for (int ia=0; ia<ucell.atoms[it].na; ia++)
+                for (int ia = 0; ia < this->na[it]; ia++)
                 {
-                    ofs << ucell.atoms[it].label
+                    ofs << atoms[it].label
                     << " " << newpos[3*at]
                     << " " << newpos[3*at+1]
                     << " " << newpos[3*at+2] << std::endl;
@@ -724,13 +711,13 @@ void Symmetry::lattice_type(
         if(convert_atoms)
         {
             int at=0;
-            for(int it=0; it<ucell.ntype; ++it)
+            for (int it = 0; it < this->ntype; ++it)
             {
-                for(int ia=0; ia<ucell.atoms[it].na; ++ia)
+                for (int ia = 0; ia < this->na[it]; ++ia)
                 {
-                    ModuleBase::Mathzone::Cartesian_to_Direct(ucell.atoms[it].tau[ia].x,
-                                    ucell.atoms[it].tau[ia].y,
-                                    ucell.atoms[it].tau[ia].z,
+                    ModuleBase::Mathzone::Cartesian_to_Direct(atoms[it].tau[ia].x,
+                        atoms[it].tau[ia].y,
+                        atoms[it].tau[ia].z,
                                     v1.x, v1.y, v1.z,
                                     v2.x, v2.y, v2.z,
                                     v3.x, v3.y, v3.z,
@@ -903,13 +890,8 @@ void Symmetry::checksym(ModuleBase::Matrix3 &s, ModuleBase::Vector3<double> &gtr
             this->check_boundary(pos[j*3+1]);
             this->check_boundary(pos[j*3+2]);
         }
-         //for( int iat =0 ; iat < ucell.nat ; iat++)
-         //std::cout << " newpos_now1 = " << newpos[3*iat] << " " << newpos[3*iat+1] << " " << newpos[3*iat+2] << std::endl;
-
         //order original atomic positions for current species
         this->atom_ordering_new(pos + istart[it] * 3, na[it], index + istart[it]);
-         //for( int iat =0 ; iat < ucell.nat ; iat++)
-         //std::cout << " newpos_now2 = " << newpos[3*iat] << " " << newpos[3*iat+1] << " " << newpos[3*iat+2] << std::endl;
 
         //Rotate atoms of current species
         for (int j = istart[it]; j < istart[it] + na[it]; ++j)
@@ -1063,7 +1045,7 @@ void Symmetry::checksym(ModuleBase::Matrix3 &s, ModuleBase::Vector3<double> &gtr
     return;
 }
 
-void Symmetry::pricell(double* pos)
+void Symmetry::pricell(double* pos, const Atom* atoms)
 {
     bool no_diff = 0;
     s_flag = 0;
@@ -1290,12 +1272,11 @@ void Symmetry::pricell(double* pos)
 #endif
 
     // get the optimized primitive cell
-    UnitCell tmp_ucell;
     std::string pbravname;
     ModuleBase::Vector3<double> p01=p1, p02=p2, p03=p3;
     double pcel_pre_const[6];
     for(int i=0;i<6;++i) pcel_pre_const[i]=pcel_const[i];
-    this->lattice_type(p1, p2, p3, p01, p02, p03, pcel_const, pcel_pre_const, pbrav, pbravname,tmp_ucell, false, nullptr);
+    this->lattice_type(p1, p2, p3, p01, p02, p03, pcel_const, pcel_pre_const, pbrav, pbravname, atoms, false, nullptr);
 
     this->plat.e11=p1.x;
     this->plat.e12=p1.y;
@@ -1636,7 +1617,7 @@ for (int g_index = 0; g_index < group_index; g_index++)
     ModuleBase::timer::tick("Symmetry","rhog_symmetry");
 }
 
-void Symmetry::set_atom_map(const UnitCell& ucell)
+void Symmetry::set_atom_map(const Atom* atoms)
 {
     ModuleBase::TITLE("Symmetry", "set_atom_map");
     if (this->isym_rotiat_.size() > 0) return;
@@ -1645,15 +1626,15 @@ void Symmetry::set_atom_map(const UnitCell& ucell)
 
     double* pos = this->newpos;
     double* rotpos = this->rotpos;
-    ModuleBase::GlobalFunc::ZEROS(pos, ucell.nat * 3);
+    ModuleBase::GlobalFunc::ZEROS(pos, this->nat * 3);
     int iat = 0;
-    for (int it = 0; it < ucell.ntype; it++)
+    for (int it = 0; it < this->ntype; it++)
     {
-        for (int ia = 0; ia < ucell.atoms[it].na; ia++)
+        for (int ia = 0; ia < this->na[it]; ia++)
         {
-            pos[3 * iat] = ucell.atoms[it].taud[ia].x;
-            pos[3 * iat + 1] = ucell.atoms[it].taud[ia].y;
-            pos[3 * iat + 2] = ucell.atoms[it].taud[ia].z;
+            pos[3 * iat] = atoms[it].taud[ia].x;
+            pos[3 * iat + 1] = atoms[it].taud[ia].y;
+            pos[3 * iat + 2] = atoms[it].taud[ia].z;
             for (int k = 0; k < 3; ++k)
             {
                 this->check_translation(pos[iat * 3 + k], -floor(pos[iat * 3 + k]));
@@ -1662,7 +1643,7 @@ void Symmetry::set_atom_map(const UnitCell& ucell)
             iat++;
         }
     }
-    for (int it = 0; it < ntype; it++)
+    for (int it = 0; it < this->ntype; it++)
     {
         for (int ia = istart[it]; ia < istart[it] + na[it]; ++ia)
         {
@@ -1729,12 +1710,12 @@ void Symmetry::symmetrize_vec3_nat(double* v)const   // pengfei 2016-12-20
 	return;
 }
 
-void Symmetry::symmetrize_mat3(ModuleBase::matrix& sigma, const UnitCell& ucell)const   //zhengdy added 2017
+void Symmetry::symmetrize_mat3(ModuleBase::matrix& sigma, const Lattice& lat)const   //zhengdy added 2017
 {
-    ModuleBase::matrix A = ucell.latvec.to_matrix();
-    ModuleBase::matrix AT = ucell.latvec.Transpose().to_matrix();
-    ModuleBase::matrix invA = ucell.GT.to_matrix();
-    ModuleBase::matrix invAT = ucell.G.to_matrix();
+    ModuleBase::matrix A = lat.latvec.to_matrix();
+    ModuleBase::matrix AT = lat.latvec.Transpose().to_matrix();
+    ModuleBase::matrix invA = lat.GT.to_matrix();
+    ModuleBase::matrix invAT = lat.G.to_matrix();
     ModuleBase::matrix tot_sigma(3, 3, true);
     sigma = A * sigma * AT;
     for (int k = 0; k < nrotk; ++k)
@@ -2094,18 +2075,18 @@ void Symmetry::hermite_normal_form(const ModuleBase::Matrix3 &s3, ModuleBase::Ma
     return;
 }
 
-bool Symmetry::magmom_same_check(const UnitCell& ucell) const
+bool Symmetry::magmom_same_check(const Atom* atoms)const
 {
     ModuleBase::TITLE("Symmetry", "magmom_same_check");
     bool pricell_loop = true;
     for (int it = 0;it < ntype;++it)
     {
         if (pricell_loop)
-            for (int ia = 1;ia < ucell.atoms[it].na;++ia)
+            for (int ia = 1;ia < atoms[it].na;++ia)
             {
-                if (!equal(ucell.atoms[it].m_loc_[ia].x, ucell.atoms[it].m_loc_[0].x) ||
-                    !equal(ucell.atoms[it].m_loc_[ia].y, ucell.atoms[it].m_loc_[0].y) ||
-                    !equal(ucell.atoms[it].m_loc_[ia].z, ucell.atoms[it].m_loc_[0].z))
+                if (!equal(atoms[it].m_loc_[ia].x, atoms[it].m_loc_[0].x) ||
+                    !equal(atoms[it].m_loc_[ia].y, atoms[it].m_loc_[0].y) ||
+                    !equal(atoms[it].m_loc_[ia].z, atoms[it].m_loc_[0].z))
                 {
                     pricell_loop = false;
                     break;
diff --git a/source/module_cell/module_symmetry/symmetry.h b/source/module_cell/module_symmetry/symmetry.h
index 669d0d4faf..d72c9d29b8 100644
--- a/source/module_cell/module_symmetry/symmetry.h
+++ b/source/module_cell/module_symmetry/symmetry.h
@@ -1,7 +1,8 @@
 #ifndef SYMMETRY_H
 #define SYMMETRY_H
 
-#include "module_cell/unitcell.h"
+#include "module_cell/unitcell_data.h"
+#include "module_cell/atom_spec.h"
 #include "symmetry_basic.h"
 
 namespace ModuleSymmetry
@@ -9,9 +10,12 @@ namespace ModuleSymmetry
 class Symmetry : public Symmetry_Basic
 {
 public:
-
-	 Symmetry();
-	~Symmetry();
+    Symmetry() {
+        this->epsilon = 1e-6;
+        this->tab = 12;
+        this->available = true;
+    };
+    ~Symmetry() {};
 
 	//symmetry flag for levels
 	//-1 : no symmetry at all, k points would be total nks in KPT
@@ -21,7 +25,7 @@ class Symmetry : public Symmetry_Basic
     static bool symm_autoclose;
     static bool pricell_loop;   ///< whether to loop primitive cell in rhog_symmetry
 
-	void analy_sys(const UnitCell &ucell, std::ofstream &ofs_running);
+    void analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, std::ofstream& ofs_running);
 	bool available;
 
 	ModuleBase::Vector3<double> s1, s2, s3;
@@ -79,9 +83,9 @@ class Symmetry : public Symmetry_Basic
 	int standard_lat(ModuleBase::Vector3<double> &a,ModuleBase::Vector3<double> &b,ModuleBase::Vector3<double> &c,double *celconst )const;
 
 	void lattice_type(ModuleBase::Vector3<double> &v1,ModuleBase::Vector3<double> &v2,ModuleBase::Vector3<double> &v3, 
-	    	ModuleBase::Vector3<double> &v01, ModuleBase::Vector3<double> &v02, ModuleBase::Vector3<double> &v03,
-			double *cel_const, double* pre_const, int& real_brav, std::string &bravname, const UnitCell &ucell, 
-			bool convert_atoms, double* newpos=nullptr)const;
+        ModuleBase::Vector3<double>& v01, ModuleBase::Vector3<double>& v02, ModuleBase::Vector3<double>& v03,
+        double* cel_const, double* pre_const, int& real_brav, std::string& bravname, const Atom* atoms,
+        bool convert_atoms, double* newpos = nullptr)const;
 
 	void recip(
 			const double a,
@@ -95,11 +99,10 @@ class Symmetry : public Symmetry_Basic
 	
 	void change_lattice(void);
 
-	// check if the input cell is a primitive cell.
-	//void pricell(const UnitCell &ucell);
 	void getgroup(int &nrot, int &nrotk, std::ofstream &ofs_running);
-	void checksym(ModuleBase::Matrix3 &s, ModuleBase::Vector3<double> &gtrans, double *pos);
-	void pricell(double* pos);
+    void checksym(ModuleBase::Matrix3& s, ModuleBase::Vector3<double>& gtrans, double* pos);
+    /// @brief  primitive cell analysis
+    void pricell(double* pos, const Atom* atoms);
 	void rho_symmetry(double *rho, const int &nr1, const int &nr2, const int &nr3);
 	void rhog_symmetry(std::complex<double> *rhogtot, int* ixyz2ipw, const int &nx, 
 			const int &ny, const int &nz, const int & fftnx, const int &fftny, const int &fftnz);
@@ -107,7 +110,7 @@ class Symmetry : public Symmetry_Basic
     /// symmetrize a vector3 with nat elements, which can be forces or variation of atom positions in relax
     void symmetrize_vec3_nat(double* v)const;
     /// symmetrize a 3*3 tensor, which can be stress or variation of unitcell in cell-relax
-    void symmetrize_mat3(ModuleBase::matrix& sigma, const UnitCell& ucell)const;
+    void symmetrize_mat3(ModuleBase::matrix& sigma, const Lattice& lat)const;
 
     void write();
 
@@ -134,7 +137,7 @@ class Symmetry : public Symmetry_Basic
 
 
     /// @brief  set atom map for each symmetry operation
-    void set_atom_map(const UnitCell& ucell);
+    void set_atom_map(const Atom* atoms);
     // to be called in lattice_type
 	void get_shortest_latvec(ModuleBase::Vector3<double> &a1, 
 			ModuleBase::Vector3<double> &a2, ModuleBase::Vector3<double> &a3)const;
@@ -144,7 +147,7 @@ class Symmetry : public Symmetry_Basic
         int& real_brav, double* cel_const, double* tmp_const)const;
 
     /// Loop the magmom of each atoms in its type when NSPIN>1. If not all the same, primitive cells should not be looped in rhog_symmetry.
-    bool magmom_same_check(const UnitCell& ucell)const;
+    bool magmom_same_check(const Atom* atoms)const;
 };
 }
 
diff --git a/source/module_cell/module_symmetry/symmetry_basic.cpp b/source/module_cell/module_symmetry/symmetry_basic.cpp
index 1ad6c1c66c..b9ba19d08b 100644
--- a/source/module_cell/module_symmetry/symmetry_basic.cpp
+++ b/source/module_cell/module_symmetry/symmetry_basic.cpp
@@ -8,16 +8,6 @@ bool ModuleSymmetry::test_brav = 0;
 
 namespace ModuleSymmetry
 {
-Symmetry_Basic::Symmetry_Basic()
-{
-	this->epsilon = 1e-6;
-}
-
-Symmetry_Basic::~Symmetry_Basic()
-{
-}
-
-
 // Find the type of bravais lattice.
 std::string Symmetry_Basic::get_brav_name(const int ibrav) const
 {
diff --git a/source/module_cell/module_symmetry/symmetry_basic.h b/source/module_cell/module_symmetry/symmetry_basic.h
index fa4db8d846..2e815147b9 100644
--- a/source/module_cell/module_symmetry/symmetry_basic.h
+++ b/source/module_cell/module_symmetry/symmetry_basic.h
@@ -14,8 +14,8 @@ class Symmetry_Basic
 {
 	public:
 
-	 Symmetry_Basic();
-	~Symmetry_Basic();
+        Symmetry_Basic() {};
+        ~Symmetry_Basic() {};
 
     double epsilon;
     double epsilon_input;   ///< the input value of symmetry_prec, should not  be changed
diff --git a/source/module_cell/module_symmetry/test/symmetry_test.h b/source/module_cell/module_symmetry/test/symmetry_test.h
index 816bbeeafe..f50e045ff3 100644
--- a/source/module_cell/module_symmetry/test/symmetry_test.h
+++ b/source/module_cell/module_symmetry/test/symmetry_test.h
@@ -1,6 +1,6 @@
 #pragma once
 #include "module_base/mathzone.h"
-#include "../symmetry.h"
+#include "module_cell/unitcell.h"
 #include "gtest/gtest.h"
 
 #define DOUBLETHRESHOLD 1e-8
diff --git a/source/module_cell/module_symmetry/test/symmetry_test_analysis.cpp b/source/module_cell/module_symmetry/test/symmetry_test_analysis.cpp
index aa4c2af3a3..5d76f74957 100644
--- a/source/module_cell/module_symmetry/test/symmetry_test_analysis.cpp
+++ b/source/module_cell/module_symmetry/test/symmetry_test_analysis.cpp
@@ -41,7 +41,7 @@ TEST_F(SymmetryTest, AnalySys)
     {
         ModuleSymmetry::Symmetry symm;
         construct_ucell(stru_lib[stru]);
-        symm.analy_sys(ucell, ofs_running);
+        symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running);
 
         //1. ibrav
         std::string ref_point_group = stru_lib[stru].point_group;
@@ -244,7 +244,7 @@ TEST_F(SymmetryTest, SG_Pricell)
         ModuleSymmetry::Symmetry symm;
         symm.epsilon = 1e-5;
         construct_ucell(supercell_lib[stru]);
-        symm.analy_sys(ucell, ofs_running);
+        symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running);
 
         std::string ref_point_group = supercell_lib[stru].point_group;
         std::string cal_point_group = symm.pgname;
diff --git a/source/module_cell/module_symmetry/test/symmetry_test_symtrz.cpp b/source/module_cell/module_symmetry/test/symmetry_test_symtrz.cpp
index e8981ec613..eee2cdaa85 100644
--- a/source/module_cell/module_symmetry/test/symmetry_test_symtrz.cpp
+++ b/source/module_cell/module_symmetry/test/symmetry_test_symtrz.cpp
@@ -66,7 +66,7 @@ TEST_F(SymmetryTest, ForceSymmetry)
     {
         ModuleSymmetry::Symmetry symm;
         construct_ucell(supercell_lib[stru]);
-        symm.analy_sys(ucell, ofs_running);
+        symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running);
 
         ModuleBase::matrix force(ucell.nat, 3, true);
         //generate random number for force and restrict to [-100,100)
@@ -97,7 +97,7 @@ TEST_F(SymmetryTest, StressSymmetry)
     {
         ModuleSymmetry::Symmetry symm;
         construct_ucell(supercell_lib[stru]);
-        symm.analy_sys(ucell, ofs_running);
+        symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running);
 
         ModuleBase::matrix stress(3, 3, true);
         //generate random number for stress and restrict to [-1e5,1e5)
@@ -105,7 +105,7 @@ TEST_F(SymmetryTest, StressSymmetry)
             for (int j = 0;j < 3;++j)
                 stress(i, j) = double(rand()) / double(RAND_MAX) * 2e5 - 1e5;
 
-        symm.symmetrize_mat3(stress, ucell);
+        symm.symmetrize_mat3(stress, ucell.lat);
         check_stress(supercell_lib[stru], stress);
 
     }
diff --git a/source/module_cell/test/klist_test.cpp b/source/module_cell/test/klist_test.cpp
index 2c17a59534..54bfa5e358 100644
--- a/source/module_cell/test/klist_test.cpp
+++ b/source/module_cell/test/klist_test.cpp
@@ -732,7 +732,7 @@ TEST_F(KlistTest, IbzKpoint)
 	ModuleSymmetry::Symmetry symm;
 	construct_ucell(stru_lib[0]);
 	GlobalV::ofs_running.open("tmp_klist_3");
-	symm.analy_sys(ucell,GlobalV::ofs_running);
+    symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
 	//read KPT
 	std::string k_file = "./support/KPT1";
 	kv->nspin = 1;
@@ -756,7 +756,7 @@ TEST_F(KlistTest, IbzKpointIsMP)
 	ModuleSymmetry::Symmetry symm;
 	construct_ucell(stru_lib[0]);
 	GlobalV::ofs_running.open("tmp_klist_4");
-	symm.analy_sys(ucell,GlobalV::ofs_running);
+    symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
 	//read KPT
 	std::string k_file = "./support/KPT1";
 	kv->nspin = 1;
diff --git a/source/module_cell/test/klist_test_para.cpp b/source/module_cell/test/klist_test_para.cpp
index 36ea7eb574..9eaeac86f4 100644
--- a/source/module_cell/test/klist_test_para.cpp
+++ b/source/module_cell/test/klist_test_para.cpp
@@ -174,8 +174,8 @@ TEST_F(KlistParaTest,Set)
 	//construct cell and symmetry
 	ModuleSymmetry::Symmetry symm;
 	construct_ucell(stru_lib[0]);
-	if(GlobalV::MY_RANK == 0) GlobalV::ofs_running.open("tmp_klist_5");
-	symm.analy_sys(GlobalC::ucell,GlobalV::ofs_running);
+    if (GlobalV::MY_RANK == 0) GlobalV::ofs_running.open("tmp_klist_5");
+    symm.analy_sys(GlobalC::ucell.lat, GlobalC::ucell.st, GlobalC::ucell.atoms, GlobalV::ofs_running);
 	//read KPT
 	std::string k_file = "./support/KPT1";
 	//set klist
@@ -210,7 +210,7 @@ TEST_F(KlistParaTest,SetAfterVC)
 	ModuleSymmetry::Symmetry symm;
 	construct_ucell(stru_lib[0]);
 	if(GlobalV::MY_RANK == 0) GlobalV::ofs_running.open("tmp_klist_6");
-	symm.analy_sys(GlobalC::ucell,GlobalV::ofs_running);
+    symm.analy_sys(GlobalC::ucell.lat, GlobalC::ucell.st, GlobalC::ucell.atoms, GlobalV::ofs_running);
 	//read KPT
 	std::string k_file = "./support/KPT1";
 	//set klist
diff --git a/source/module_cell/unitcell.h b/source/module_cell/unitcell.h
index f47903da8b..6223e3b9f6 100644
--- a/source/module_cell/unitcell.h
+++ b/source/module_cell/unitcell.h
@@ -3,17 +3,16 @@
 
 #include "module_base/global_function.h"
 #include "module_base/global_variable.h"
-#include "module_base/matrix3.h"
-#include "module_base/intarray.h"
 #include "module_io/output.h"
 #include "module_elecstate/magnetism.h"
-#include "atom_spec.h"
+#include "module_symmetry/symmetry.h"
 
 #ifdef __LCAO
 #include "module_basis/module_ao/ORB_read.h"
 #include "setup_nonlocal.h"
 #endif
 
+
 // provide the basic information about unitcell.
 class UnitCell
 {
@@ -27,46 +26,37 @@ class UnitCell
 	double *atom_mag;
 	int n_mag_at;
 
-    int ntype;// number of atom species in UnitCell
-    int nat; // total number of atoms of all species in unitcell
-    std::string Coordinate; // "Direct" or "Cartesian" or "Cartesian_angstrom"
-    std::string latName; // Lattice name
-    double lat0; // Lattice constant(bohr)(a.u.)
-    double lat0_angstrom;// Lattice constant(angstrom)
-    double tpiba;// 2*pi / lat0;
-    double tpiba2; // tpiba ^ 2
-    double omega;// the volume of the unit cell
-
-    ModuleBase::Matrix3 latvec; // Unitcell lattice vectors
-	int *lc;  // Change the lattice vectors or not
-	ModuleBase::Vector3<double> a1,a2,a3; // Same as latvec, just at another form.
-	ModuleBase::Vector3<double> latcenter; // (a1+a2+a3)/2 the center of vector
-    ModuleBase::Matrix3 latvec_supercell; // Supercell lattice vectors
-    ModuleBase::Matrix3 G; // reciprocal lattice vector (2pi*inv(R) )
-    ModuleBase::Matrix3 GT; // traspose of G
-    ModuleBase::Matrix3 GGT; // GGT = G*GT
-    ModuleBase::Matrix3 invGGT; // inverse G
-
-    //========================================================
-    // relationship between:
-    // ntype, it
-    // nat, iat
-    // atoms[it].na, ia,
-    // atoms[it].nw, iw
-    //
-    // if know it ==> atoms[it].na; atoms[it].nw
-    // if know iat ==> it; ia;
-    // if know ia, mush have known it ==> iat
-    // if know iwt, must have known it, ia ==> iwt
-    //========================================================
-    int namax;// the max na among all atom species
-    int nwmax;// the max nw among all atom species
-
-    int *iat2it; //iat==>it, distinguish a atom belong to which type
-    int *iat2ia; //iat==>ia
-	int *iwt2iat; // iwt ==> iat.
-	int *iwt2iw; // iwt ==> iw, Peize Lin add 2018-07-02
-    ModuleBase::IntArray itia2iat;//(it, ia)==>iat, the index in nat, add 2009-3-2 by mohan
+    std::string& Coordinate = lat.Coordinate;
+    std::string& latName = lat.latName;
+    double& lat0 = lat.lat0;
+    double& lat0_angstrom = lat.lat0_angstrom;
+    double& tpiba = lat.tpiba;
+    double& tpiba2 = lat.tpiba2;
+    double& omega = lat.omega;
+    int*& lc = lat.lc;
+
+    Lattice lat;
+    ModuleBase::Matrix3& latvec = lat.latvec;
+    ModuleBase::Vector3<double>& a1 = lat.a1, & a2 = lat.a2, & a3 = lat.a3;
+    ModuleBase::Vector3<double>& latcenter = lat.latcenter;
+    ModuleBase::Matrix3& latvec_supercell = lat.latvec_supercell;
+    ModuleBase::Matrix3& G = lat.G;
+    ModuleBase::Matrix3& GT = lat.GT;
+    ModuleBase::Matrix3& GGT = lat.GGT;
+    ModuleBase::Matrix3& invGGT = lat.invGGT;
+
+    Statistics st;
+    int& ntype = st.ntype;
+    int& nat = st.nat;
+    int*& iat2it = st.iat2it;
+    int*& iat2ia = st.iat2ia;
+    int*& iwt2iat = st.iwt2iat;
+    int*& iwt2iw = st.iwt2iw;
+    ModuleBase::IntArray& itia2iat = st.itia2iat;
+    int& namax = st.namax;
+    int& nwmax = st.nwmax;
+
+    ModuleSymmetry::Symmetry symm;
 
     // ========================================================
     // iat2iwt is the atom index iat to the first global index for orbital of this atom
diff --git a/source/module_cell/unitcell_data.h b/source/module_cell/unitcell_data.h
new file mode 100644
index 0000000000..4ed8debd28
--- /dev/null
+++ b/source/module_cell/unitcell_data.h
@@ -0,0 +1,49 @@
+#include "module_base/matrix3.h"
+#include "module_base/intarray.h"
+/// @brief info of lattice 
+struct Lattice
+{
+    std::string Coordinate; // "Direct" or "Cartesian" or "Cartesian_angstrom"
+    std::string latName; // Lattice name
+    double lat0; // Lattice constant(bohr)(a.u.)
+    double lat0_angstrom;// Lattice constant(angstrom)
+    double tpiba;// 2*pi / lat0;
+    double tpiba2; // tpiba ^ 2
+    double omega;// the volume of the unit cell
+    int* lc;  // Change the lattice vectors or not
+
+    ModuleBase::Matrix3 latvec; // Unitcell lattice vectors
+    ModuleBase::Vector3<double> a1, a2, a3; // Same as latvec, just at another form.
+    ModuleBase::Vector3<double> latcenter; // (a1+a2+a3)/2 the center of vector
+    ModuleBase::Matrix3 latvec_supercell; // Supercell lattice vectors
+    ModuleBase::Matrix3 G; // reciprocal lattice vector (2pi*inv(R) )
+    ModuleBase::Matrix3 GT; // traspose of G
+    ModuleBase::Matrix3 GGT; // GGT = G*GT
+    ModuleBase::Matrix3 invGGT; // inverse G
+};
+
+//========================================================
+// relationship between:
+// ntype, it
+// nat, iat
+// atoms[it].na, ia,
+// atoms[it].nw, iw
+//
+// if know it ==> atoms[it].na; atoms[it].nw
+// if know iat ==> it; ia;
+// if know ia, mush have known it ==> iat
+// if know iwt, must have known it, ia ==> iwt
+//========================================================
+/// @brief usefull data and index maps
+struct Statistics
+{
+    int ntype;// number of atom species in UnitCell
+    int nat; // total number of atoms of all species in unitcell
+    int* iat2it; //iat==>it, distinguish a atom belong to which type
+    int* iat2ia; //iat==>ia
+    int* iwt2iat; // iwt ==> iat.
+    int* iwt2iw; // iwt ==> iw, Peize Lin add 2018-07-02
+    ModuleBase::IntArray itia2iat;//(it, ia)==>iat, the index in nat, add 2009-3-2 by mohan
+    int namax;// the max na among all atom species
+    int nwmax;// the max nw among all atom species
+};
\ No newline at end of file
diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp
index 4b6d2f4e80..9c231acf94 100644
--- a/source/module_esolver/esolver_fp.cpp
+++ b/source/module_esolver/esolver_fp.cpp
@@ -24,7 +24,7 @@ namespace ModuleESolver
         pw_big->setbxyz(INPUT.bx, INPUT.by, INPUT.bz);
         sf.set(pw_rhod, INPUT.nbspline);
 
-        this->symm.epsilon = this->symm.epsilon_input = INPUT.symmetry_prec;
+        GlobalC::ucell.symm.epsilon = GlobalC::ucell.symm.epsilon_input = INPUT.symmetry_prec;
 }
     ESolver_FP::~ESolver_FP()
     {
@@ -131,11 +131,11 @@ namespace ModuleESolver
 
         if(ModuleSymmetry::Symmetry::symm_flag == 1)
         {
-            symm.analy_sys(cell, GlobalV::ofs_running);
+            cell.symm.analy_sys(cell.lat, cell.st, cell.atoms, GlobalV::ofs_running);
             ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
         }
 
-        kv.set_after_vc(symm, GlobalV::global_kpoint_card, GlobalV::NSPIN, cell.G, cell.latvec);
+        kv.set_after_vc(cell.symm, GlobalV::global_kpoint_card, GlobalV::NSPIN, cell.G, cell.latvec);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
     }
 
diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h
index 79cdbafd52..21b0588ff2 100644
--- a/source/module_esolver/esolver_fp.h
+++ b/source/module_esolver/esolver_fp.h
@@ -30,7 +30,6 @@ namespace ModuleESolver
 
         elecstate::ElecState* pelec = nullptr;
         Charge chr;
-        ModuleSymmetry::Symmetry symm;
         //--------------temporary----------------------------
         // this is the interface of non-self-consistant calculation
         virtual void nscf(){};
diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp
index e532492cd0..a97ebfd9dd 100644
--- a/source/module_esolver/esolver_ks.cpp
+++ b/source/module_esolver/esolver_ks.cpp
@@ -182,12 +182,12 @@ namespace ModuleESolver
         // symmetry analysis should be performed every time the cell is changed
         if (ModuleSymmetry::Symmetry::symm_flag == 1)
         {
-            this->symm.analy_sys(ucell, GlobalV::ofs_running);
+            ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
             ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
         }
 
         // Setup the k points according to symmetry.
-        this->kv.set(this->symm, GlobalV::global_kpoint_card, GlobalV::NSPIN, ucell.G, ucell.latvec);
+        this->kv.set(ucell.symm, GlobalV::global_kpoint_card, GlobalV::NSPIN, ucell.G, ucell.latvec);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
 
         // print information
diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
index 190f4215d3..e192ec3af6 100644
--- a/source/module_esolver/esolver_ks_lcao.cpp
+++ b/source/module_esolver/esolver_ks_lcao.cpp
@@ -82,12 +82,12 @@ namespace ModuleESolver
 
         if (ModuleSymmetry::Symmetry::symm_flag == 1)
         {
-            this->symm.analy_sys(ucell, GlobalV::ofs_running);
+            ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
             ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
         }
 
         // Setup the k points according to symmetry.
-        this->kv.set(this->symm, GlobalV::global_kpoint_card, GlobalV::NSPIN, ucell.G, ucell.latvec);
+        this->kv.set(ucell.symm, GlobalV::global_kpoint_card, GlobalV::NSPIN, ucell.G, ucell.latvec);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
 
         Print_Info::setup_parameters(ucell, this->kv);
@@ -280,7 +280,7 @@ namespace ModuleESolver
         * this->exx_lri_double,
         * this->exx_lri_complex,
 #endif  
-        & this->symm);
+        & GlobalC::ucell.symm);
     // delete RA after cal_Force
     this->RA.delete_grid();
     this->have_force = true;
@@ -650,7 +650,7 @@ namespace ModuleESolver
     Symmetry_rho srho;
     for (int is = 0; is < GlobalV::NSPIN; is++)
     {
-        srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::Pgrid, this->symm);
+        srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm);
     }
 
     // (6) compute magnetization, only for spin==2
diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp
index b85dfd31db..1fb831dbb7 100644
--- a/source/module_esolver/esolver_ks_lcao_elec.cpp
+++ b/source/module_esolver/esolver_ks_lcao_elec.cpp
@@ -299,7 +299,7 @@ void ESolver_KS_LCAO<TK, TR>::beforescf(int istep)
     Symmetry_rho srho;
     for (int is = 0; is < GlobalV::NSPIN; is++)
     {
-        srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::Pgrid, this->symm);
+        srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm);
     }
 // Peize Lin add 2016-12-03
 #ifdef __EXX
diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp
index 78fa01e76e..a26e97f75e 100644
--- a/source/module_esolver/esolver_ks_lcao_tddft.cpp
+++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp
@@ -223,7 +223,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(int istep, int iter, double ethr)
         Symmetry_rho srho;
         for (int is = 0; is < GlobalV::NSPIN; is++)
         {
-            srho.begin(is, *(pelec->charge), pw_rho, GlobalC::Pgrid, this->symm);
+            srho.begin(is, *(pelec->charge), pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm);
         }
     }
 
diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
index 82ca5ada93..d814cf8406 100644
--- a/source/module_esolver/esolver_ks_pw.cpp
+++ b/source/module_esolver/esolver_ks_pw.cpp
@@ -407,7 +407,7 @@ void ESolver_KS_PW<T, Device>::beforescf(int istep)
     Symmetry_rho srho;
     for (int is = 0; is < GlobalV::NSPIN; is++)
     {
-        srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::Pgrid, this->symm);
+        srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::Pgrid, GlobalC::ucell.symm);
     }
 
     // liuyu move here 2023-10-09
@@ -706,7 +706,7 @@ void ESolver_KS_PW<T, Device>::hamilt2density(const int istep, const int iter, c
     Symmetry_rho srho;
     for (int is = 0; is < GlobalV::NSPIN; is++)
     {
-        srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::Pgrid, this->symm);
+        srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::Pgrid, GlobalC::ucell.symm);
     }
 
     // compute magnetization, only for LSDA(spin==2)
@@ -857,7 +857,7 @@ void ESolver_KS_PW<T, Device>::cal_Force(ModuleBase::matrix& force)
                                ? new psi::Psi<std::complex<double>, Device>(this->kspw_psi[0])
                                : reinterpret_cast<psi::Psi<std::complex<double>, Device>*>(this->kspw_psi);
     }
-    ff.cal_force(force, *this->pelec, this->pw_rhod, &this->symm, &this->sf, &this->kv, this->pw_wfc, this->__kspw_psi);
+    ff.cal_force(force, *this->pelec, this->pw_rhod, &GlobalC::ucell.symm, &this->sf, &this->kv, this->pw_wfc, this->__kspw_psi);
 }
 
 template <typename T, typename Device>
@@ -875,7 +875,7 @@ void ESolver_KS_PW<T, Device>::cal_Stress(ModuleBase::matrix& stress)
     ss.cal_stress(stress,
                   GlobalC::ucell,
                   this->pw_rhod,
-                  &this->symm,
+        &GlobalC::ucell.symm,
                   &this->sf,
                   &this->kv,
                   this->pw_wfc,
diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp
index 915a884c48..53e1224e44 100644
--- a/source/module_esolver/esolver_of.cpp
+++ b/source/module_esolver/esolver_of.cpp
@@ -83,13 +83,13 @@ void ESolver_OF::Init(Input& inp, UnitCell& ucell)
     // symmetry analysis should be performed every time the cell is changed
     if (ModuleSymmetry::Symmetry::symm_flag == 1)
     {
-        this->symm.analy_sys(ucell, GlobalV::ofs_running);
+        ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
     }
 
     // Setup the k points according to symmetry.
-    kv.set(this->symm, GlobalV::global_kpoint_card, GlobalV::NSPIN, ucell.G, ucell.latvec);
-    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
+    kv.set(ucell.symm, GlobalV::global_kpoint_card, GlobalV::NSPIN, ucell.G, ucell.latvec);
+    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"INIT K-POINTS");
 
     // print information
     // mohan add 2021-01-30
@@ -286,7 +286,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
     Symmetry_rho srho;
     for (int is = 0; is < GlobalV::NSPIN; is++)
     {
-        srho.begin(is, *(pelec->charge), this->pw_rho, GlobalC::Pgrid, this->symm);
+        srho.begin(is, *(pelec->charge), this->pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm);
     }
 
     for (int is = 0; is < GlobalV::NSPIN; ++is)
@@ -433,7 +433,7 @@ void ESolver_OF::update_rho()
     //     Symmetry_rho srho;
     //     for (int is = 0; is < GlobalV::NSPIN; is++)
     //     {
-    //         srho.begin(is, *(pelec->charge), this->pw_rho, GlobalC::Pgrid, this->symm);
+    //         srho.begin(is, *(pelec->charge), this->pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm);
     //         for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs)
     //         {
     //             this->pphi_[is][ibs] = sqrt(pelec->charge->rho[is][ibs]);
@@ -624,7 +624,7 @@ double ESolver_OF::cal_Energy()
 void ESolver_OF::cal_Force(ModuleBase::matrix& force)
 {
     Forces<double> ff(GlobalC::ucell.nat);
-    ff.cal_force(force, *pelec, this->pw_rho, &this->symm, &sf);
+    ff.cal_force(force, *pelec, this->pw_rho, &GlobalC::ucell.symm, &sf);
 }
 
 /**
@@ -639,6 +639,6 @@ void ESolver_OF::cal_Stress(ModuleBase::matrix& stress)
     this->kinetic_stress(kinetic_stress_);
 
     OF_Stress_PW ss(this->pelec, this->pw_rho);
-    ss.cal_stress(stress, kinetic_stress_, GlobalC::ucell, &this->symm, &sf, &kv);
+    ss.cal_stress(stress, kinetic_stress_, GlobalC::ucell, &GlobalC::ucell.symm, &sf, &kv);
 }
 } // namespace ModuleESolver
\ No newline at end of file
diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp
index c5b033dc59..26390d99b5 100644
--- a/source/module_esolver/esolver_sdft_pw.cpp
+++ b/source/module_esolver/esolver_sdft_pw.cpp
@@ -181,7 +181,7 @@ void ESolver_SDFT_PW::hamilt2density(int istep, int iter, double ethr)
         Symmetry_rho srho;
         for (int is = 0; is < GlobalV::NSPIN; is++)
         {
-            srho.begin(is, *(this->pelec->charge), pw_rho, GlobalC::Pgrid, this->symm);
+            srho.begin(is, *(this->pelec->charge), pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm);
         }
         this->pelec->f_en.deband = this->pelec->cal_delta_eband();
     }
@@ -202,12 +202,12 @@ double ESolver_SDFT_PW::cal_Energy()
 void ESolver_SDFT_PW::cal_Force(ModuleBase::matrix& force)
 {
     Sto_Forces ff(GlobalC::ucell.nat);
-    ff.cal_stoforce(force, *this->pelec, pw_rho, &this->symm, &sf, &kv, pw_wfc, this->psi, this->stowf);
+    ff.cal_stoforce(force, *this->pelec, pw_rho, &GlobalC::ucell.symm, &sf, &kv, pw_wfc, this->psi, this->stowf);
 }
 void ESolver_SDFT_PW::cal_Stress(ModuleBase::matrix& stress)
 {
     Sto_Stress_PW ss;
-    ss.cal_stress(stress, *this->pelec, pw_rho, &this->symm, &sf, &kv, pw_wfc, this->psi, this->stowf, pelec->charge);
+    ss.cal_stress(stress, *this->pelec, pw_rho, &GlobalC::ucell.symm, &sf, &kv, pw_wfc, this->psi, this->stowf, pelec->charge);
 }
 void ESolver_SDFT_PW::postprocess()
 {
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
index fdd8e01a62..eca982cbf1 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
@@ -544,7 +544,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
         }
         if (ModuleSymmetry::Symmetry::symm_flag == 1)
         {
-            symm->symmetrize_mat3(scs, GlobalC::ucell);
+            symm->symmetrize_mat3(scs, GlobalC::ucell.lat);
         } // end symmetry
 
 #ifdef __DEEPKS
@@ -558,7 +558,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
         {
             if (ModuleSymmetry::Symmetry::symm_flag == 1)
             {
-                symm->symmetrize_mat3(svnl_dalpha, GlobalC::ucell);
+                symm->symmetrize_mat3(svnl_dalpha, GlobalC::ucell.lat);
             } // end symmetry
             for (int i = 0; i < 3; i++)
             {
diff --git a/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp b/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp
index 5665217548..10a08feeff 100644
--- a/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp
+++ b/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp
@@ -95,7 +95,7 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot,
 
     if (ModuleSymmetry::Symmetry::symm_flag == 1)
     {
-        p_symm->symmetrize_mat3(sigmatot, ucell);
+        p_symm->symmetrize_mat3(sigmatot, ucell.lat);
     }
 
     bool ry = false;
diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_kin.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_kin.cpp
index c9f767245f..49006099d1 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_kin.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_kin.cpp
@@ -149,7 +149,7 @@ void Stress_Func<FPTYPE, Device>::stress_kin(ModuleBase::matrix& sigma,
 	//do symmetry
     if (ModuleSymmetry::Symmetry::symm_flag == 1)
     {
-        p_symm->symmetrize_mat3(sigma, GlobalC::ucell);
+        p_symm->symmetrize_mat3(sigma, GlobalC::ucell.lat);
     } // end symmetry
 
     delete[] gk[0];
diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_nl.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_nl.cpp
index a36ab13ca7..8a3d21a515 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_nl.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_nl.cpp
@@ -301,7 +301,7 @@ void Stress_Func<FPTYPE, Device>::stress_nl(ModuleBase::matrix& sigma,
 	//do symmetry
     if (ModuleSymmetry::Symmetry::symm_flag == 1)
     {
-        p_symm->symmetrize_mat3(sigma, GlobalC::ucell);
+        p_symm->symmetrize_mat3(sigma, GlobalC::ucell.lat);
     } // end symmetry
 
     delete [] h_atom_nh;
diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp
index bdbcdea57c..5386d11fcd 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp
@@ -123,7 +123,7 @@ void Stress_PW<FPTYPE, Device>::cal_stress(ModuleBase::matrix& sigmatot,
     
 	if(ModuleSymmetry::Symmetry::symm_flag == 1)                          
 	{
-        p_symm->symmetrize_mat3(sigmatot, ucell);
+        p_symm->symmetrize_mat3(sigmatot, ucell.lat);
     }
 
 	bool ry = false;
diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp
index 3ad989be74..e2299807ea 100644
--- a/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp
+++ b/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp
@@ -65,7 +65,7 @@ void Sto_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot,
 
     if (ModuleSymmetry::Symmetry::symm_flag == 1)
     {
-        p_symm->symmetrize_mat3(sigmatot, GlobalC::ucell);
+        p_symm->symmetrize_mat3(sigmatot, GlobalC::ucell.lat);
     }
 
     bool ry = false;
@@ -196,7 +196,7 @@ void Sto_Stress_PW::sto_stress_kin(ModuleBase::matrix& sigma,
     // do symmetry
     if (ModuleSymmetry::Symmetry::symm_flag == 1)
     {
-        p_symm->symmetrize_mat3(sigma, GlobalC::ucell);
+        p_symm->symmetrize_mat3(sigma, GlobalC::ucell.lat);
     }
     delete[] gk[0];
     delete[] gk[1];
@@ -461,7 +461,7 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma,
     // do symmetry
     if (ModuleSymmetry::Symmetry::symm_flag == 1)
     {
-        p_symm->symmetrize_mat3(sigma, GlobalC::ucell);
+        p_symm->symmetrize_mat3(sigma, GlobalC::ucell.lat);
     }
 
     //  this->print(ofs_running, "nonlocal stress", stresnl);
diff --git a/source/module_io/test/for_testing_input_conv.h b/source/module_io/test/for_testing_input_conv.h
index f39953f7a1..b444732606 100644
--- a/source/module_io/test/for_testing_input_conv.h
+++ b/source/module_io/test/for_testing_input_conv.h
@@ -187,18 +187,6 @@ Structure_Factor::Structure_Factor()
 Structure_Factor::~Structure_Factor()
 {
 }
-ModuleSymmetry::Symmetry::Symmetry()
-{
-}
-ModuleSymmetry::Symmetry::~Symmetry()
-{
-}
-ModuleSymmetry::Symmetry_Basic::Symmetry_Basic()
-{
-}
-ModuleSymmetry::Symmetry_Basic::~Symmetry_Basic()
-{
-}
 WF_atomic::WF_atomic()
 {
 }

From bc1cf7cac3d1518eacf7611fccc6fb353a41d031 Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Thu, 23 Nov 2023 01:04:15 +0800
Subject: [PATCH 2/4] symmetrization for cell-relax

---
 source/module_relax/relax_new/relax.cpp       | 31 +++++++++++++------
 .../module_relax/relax_new/test/relax_test.h  |  5 +--
 .../relax_old/lattice_change_basic.cpp        |  9 ++++++
 source/module_relax/relax_old/test/for_test.h |  2 ++
 4 files changed, 35 insertions(+), 12 deletions(-)

diff --git a/source/module_relax/relax_new/relax.cpp b/source/module_relax/relax_new/relax.cpp
index 8972c86f39..95b2d66c7d 100644
--- a/source/module_relax/relax_new/relax.cpp
+++ b/source/module_relax/relax_new/relax.cpp
@@ -462,16 +462,27 @@ void Relax::move_cell_ions(const bool is_new_dir)
     {
         // imo matrix3 class is not a very clever way to store 3*3 matrix ...
         ModuleBase::Matrix3 sr_dr_cell;
-
-        sr_dr_cell.e11 = search_dr_cell(0,0);
-        sr_dr_cell.e12 = search_dr_cell(0,1);
-        sr_dr_cell.e13 = search_dr_cell(0,2);
-        sr_dr_cell.e21 = search_dr_cell(1,0);
-        sr_dr_cell.e22 = search_dr_cell(1,1);
-        sr_dr_cell.e23 = search_dr_cell(1,2);
-        sr_dr_cell.e31 = search_dr_cell(2,0);
-        sr_dr_cell.e32 = search_dr_cell(2,1);
-        sr_dr_cell.e33 = search_dr_cell(2,2);
+        auto cp_mat_to_mat3 = [&sr_dr_cell, this]() -> void
+            {
+                sr_dr_cell.e11 = search_dr_cell(0, 0);
+                sr_dr_cell.e12 = search_dr_cell(0, 1);
+                sr_dr_cell.e13 = search_dr_cell(0, 2);
+                sr_dr_cell.e21 = search_dr_cell(1, 0);
+                sr_dr_cell.e22 = search_dr_cell(1, 1);
+                sr_dr_cell.e23 = search_dr_cell(1, 2);
+                sr_dr_cell.e31 = search_dr_cell(2, 0);
+                sr_dr_cell.e32 = search_dr_cell(2, 1);
+                sr_dr_cell.e33 = search_dr_cell(2, 2);
+            };
+        cp_mat_to_mat3();
+
+        if (ModuleSymmetry::Symmetry::symm_flag)
+        {
+            search_dr_cell = sr_dr_cell.Transpose().to_matrix();
+            GlobalC::ucell.symm.symmetrize_mat3(search_dr_cell, GlobalC::ucell.lat);
+            cp_mat_to_mat3();
+            sr_dr_cell = sr_dr_cell.Transpose();
+        }
 
         // The logic here is as follows: a line search is a continuation
         // in the new direction; but GlobalC::ucell.latvec now is already
diff --git a/source/module_relax/relax_new/test/relax_test.h b/source/module_relax/relax_new/test/relax_test.h
index 2666b1ad77..d6a77d9faa 100644
--- a/source/module_relax/relax_new/test/relax_test.h
+++ b/source/module_relax/relax_new/test/relax_test.h
@@ -43,8 +43,9 @@ Atom_pseudo::Atom_pseudo(){};
 Atom_pseudo::~Atom_pseudo(){};
 pseudo::pseudo(){};
 pseudo::~pseudo(){};
-
-Structure_Factor::Structure_Factor(){};
+int ModuleSymmetry::Symmetry::symm_flag = 0;
+void ModuleSymmetry::Symmetry::symmetrize_mat3(ModuleBase::matrix& sigma, const Lattice& lat)const {};
+Structure_Factor::Structure_Factor() {};
 Structure_Factor::~Structure_Factor(){};
 void Structure_Factor::setup_structure_factor(UnitCell* Ucell, const ModulePW::PW_Basis* rho_basis){};
 
diff --git a/source/module_relax/relax_old/lattice_change_basic.cpp b/source/module_relax/relax_old/lattice_change_basic.cpp
index 1968241812..3cd8543b0f 100644
--- a/source/module_relax/relax_old/lattice_change_basic.cpp
+++ b/source/module_relax/relax_old/lattice_change_basic.cpp
@@ -82,6 +82,15 @@ void Lattice_Change_Basic::change_lattice(UnitCell &ucell, double *move, double
        "<<std::setprecision(12)<<ucell.latvec.e31<<"   "<<ucell.latvec.e32<<"
        "<<ucell.latvec.e33<<std::endl;
     */
+    if (ModuleSymmetry::Symmetry::symm_flag)
+    {
+        ModuleBase::matrix move_mat_t(3, 3);
+        for (int i = 0;i < 3;++i)for (int j = 0;j < 3;++j)move_mat_t(j, i) = move[i * 3 + j] / ucell.lat0;    //transpose
+        ModuleBase::matrix symm_move_mat_t = (move_mat_t * ucell.G.to_matrix());//symmetrize (latvec^{-1} * move_mat)^T
+        ucell.symm.symmetrize_mat3(symm_move_mat_t, ucell.lat);
+        move_mat_t = symm_move_mat_t * ucell.latvec.Transpose().to_matrix();//G^{-1}=latvec^T
+        for (int i = 0;i < 3;++i)for (int j = 0;j < 3;++j)move[i * 3 + j] = move_mat_t(j, i) * ucell.lat0;//transpose back
+    }
 
     if (ucell.lc[0] != 0)
     {
diff --git a/source/module_relax/relax_old/test/for_test.h b/source/module_relax/relax_old/test/for_test.h
index 61ed8112c6..40edd0a640 100644
--- a/source/module_relax/relax_old/test/for_test.h
+++ b/source/module_relax/relax_old/test/for_test.h
@@ -103,5 +103,7 @@ pseudo::pseudo()
 pseudo::~pseudo()
 {
 }
+int ModuleSymmetry::Symmetry::symm_flag = 0;
+void ModuleSymmetry::Symmetry::symmetrize_mat3(ModuleBase::matrix& sigma, const Lattice& lat)const {};
 
 #endif
\ No newline at end of file

From dee73a86cf2a0df7c6aa3b92b95936107a3a8741 Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Thu, 23 Nov 2023 03:30:24 +0800
Subject: [PATCH 3/4] relax

---
 .../module_cell/module_symmetry/symmetry.cpp  | 21 +++++++++++++++++--
 source/module_cell/module_symmetry/symmetry.h |  8 ++++++-
 source/module_relax/relax_new/relax.cpp       |  2 ++
 .../module_relax/relax_new/test/relax_test.h  |  1 +
 .../relax_old/ions_move_basic.cpp             |  3 +++
 source/module_relax/relax_old/test/for_test.h |  1 +
 6 files changed, 33 insertions(+), 3 deletions(-)

diff --git a/source/module_cell/module_symmetry/symmetry.cpp b/source/module_cell/module_symmetry/symmetry.cpp
index b25146e080..5bb91e4b8b 100644
--- a/source/module_cell/module_symmetry/symmetry.cpp
+++ b/source/module_cell/module_symmetry/symmetry.cpp
@@ -233,7 +233,9 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms,
 
     if (GlobalV::NSPIN > 1) pricell_loop = this->magmom_same_check(atoms);
 
-	delete[] newpos;
+    if (GlobalV::CALCULATION == "relax") this->all_mbl = this->is_all_movable(atoms, st);
+
+    delete[] newpos;
     delete[] na;
     delete[] rotpos;
     delete[] index;
@@ -2096,4 +2098,19 @@ bool Symmetry::magmom_same_check(const Atom* atoms)const
     return pricell_loop;
 }
 
-}
\ No newline at end of file
+bool Symmetry::is_all_movable(const Atom* atoms, const Statistics& st)const
+{
+    bool all_mbl = true;
+    for (int iat = 0;iat < st.nat;++iat)
+    {
+        int it = st.iat2it[iat];
+        int ia = st.iat2ia[iat];
+        if (!atoms[it].mbl[ia].x || !atoms[it].mbl[ia].y || !atoms[it].mbl[ia].z)
+        {
+            all_mbl = false;
+            break;
+        }
+    }
+    return all_mbl;
+}
+}
diff --git a/source/module_cell/module_symmetry/symmetry.h b/source/module_cell/module_symmetry/symmetry.h
index d72c9d29b8..bf00c758af 100644
--- a/source/module_cell/module_symmetry/symmetry.h
+++ b/source/module_cell/module_symmetry/symmetry.h
@@ -80,7 +80,9 @@ class Symmetry : public Symmetry_Basic
 
 	int tab;
 
-	int standard_lat(ModuleBase::Vector3<double> &a,ModuleBase::Vector3<double> &b,ModuleBase::Vector3<double> &c,double *celconst )const;
+    bool all_mbl = true;    ///< whether all the atoms are movable in all the directions
+
+    int standard_lat(ModuleBase::Vector3<double>& a, ModuleBase::Vector3<double>& b, ModuleBase::Vector3<double>& c, double* celconst)const;
 
 	void lattice_type(ModuleBase::Vector3<double> &v1,ModuleBase::Vector3<double> &v2,ModuleBase::Vector3<double> &v3, 
         ModuleBase::Vector3<double>& v01, ModuleBase::Vector3<double>& v02, ModuleBase::Vector3<double>& v03,
@@ -138,6 +140,10 @@ class Symmetry : public Symmetry_Basic
 
     /// @brief  set atom map for each symmetry operation
     void set_atom_map(const Atom* atoms);
+    /// @brief check if all the atoms are movable
+    ///  delta_pos symmetrization in relax is only meaningful when all the atoms are movable in all the directions.
+    bool is_all_movable(const Atom* atoms, const Statistics& st)const;
+
     // to be called in lattice_type
 	void get_shortest_latvec(ModuleBase::Vector3<double> &a1, 
 			ModuleBase::Vector3<double> &a2, ModuleBase::Vector3<double> &a3)const;
diff --git a/source/module_relax/relax_new/relax.cpp b/source/module_relax/relax_new/relax.cpp
index 95b2d66c7d..1fe0a6b98c 100644
--- a/source/module_relax/relax_new/relax.cpp
+++ b/source/module_relax/relax_new/relax.cpp
@@ -558,6 +558,8 @@ void Relax::move_cell_ions(const bool is_new_dir)
         }
     }
 
+    if (ModuleSymmetry::Symmetry::symm_flag && GlobalC::ucell.symm.all_mbl)GlobalC::ucell.symm.symmetrize_vec3_nat(move_ion);
+
 	GlobalC::ucell.update_pos_taud(move_ion);
 
 	// Print the structure file.
diff --git a/source/module_relax/relax_new/test/relax_test.h b/source/module_relax/relax_new/test/relax_test.h
index d6a77d9faa..4415e831df 100644
--- a/source/module_relax/relax_new/test/relax_test.h
+++ b/source/module_relax/relax_new/test/relax_test.h
@@ -45,6 +45,7 @@ pseudo::pseudo(){};
 pseudo::~pseudo(){};
 int ModuleSymmetry::Symmetry::symm_flag = 0;
 void ModuleSymmetry::Symmetry::symmetrize_mat3(ModuleBase::matrix& sigma, const Lattice& lat)const {};
+void ModuleSymmetry::Symmetry::symmetrize_vec3_nat(double* v)const {};
 Structure_Factor::Structure_Factor() {};
 Structure_Factor::~Structure_Factor(){};
 void Structure_Factor::setup_structure_factor(UnitCell* Ucell, const ModulePW::PW_Basis* rho_basis){};
diff --git a/source/module_relax/relax_old/ions_move_basic.cpp b/source/module_relax/relax_old/ions_move_basic.cpp
index 2d646cd334..7a56ac8281 100644
--- a/source/module_relax/relax_old/ions_move_basic.cpp
+++ b/source/module_relax/relax_old/ions_move_basic.cpp
@@ -92,6 +92,9 @@ void Ions_Move_Basic::move_atoms(UnitCell &ucell, double *move, double *pos)
 
     const double move_threshold = 1.0e-10;
     const int total_freedom = ucell.nat * 3;
+
+    if (ModuleSymmetry::Symmetry::symm_flag && ucell.symm.all_mbl)ucell.symm.symmetrize_vec3_nat(move);
+
     for (int i = 0; i < total_freedom; i++)
     {
         if (std::abs(move[i]) > move_threshold)
diff --git a/source/module_relax/relax_old/test/for_test.h b/source/module_relax/relax_old/test/for_test.h
index 40edd0a640..7c48614aea 100644
--- a/source/module_relax/relax_old/test/for_test.h
+++ b/source/module_relax/relax_old/test/for_test.h
@@ -105,5 +105,6 @@ pseudo::~pseudo()
 }
 int ModuleSymmetry::Symmetry::symm_flag = 0;
 void ModuleSymmetry::Symmetry::symmetrize_mat3(ModuleBase::matrix& sigma, const Lattice& lat)const {};
+void ModuleSymmetry::Symmetry::symmetrize_vec3_nat(double* v)const {};
 
 #endif
\ No newline at end of file

From 59df6a151dae6667411616e139ecfcefaa77726e Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Thu, 23 Nov 2023 13:44:52 +0800
Subject: [PATCH 4/4] fix a segfalt

---
 source/module_cell/module_symmetry/symmetry.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/source/module_cell/module_symmetry/symmetry.cpp b/source/module_cell/module_symmetry/symmetry.cpp
index 5bb91e4b8b..fa6c9f1eef 100644
--- a/source/module_cell/module_symmetry/symmetry.cpp
+++ b/source/module_cell/module_symmetry/symmetry.cpp
@@ -1622,7 +1622,7 @@ for (int g_index = 0; g_index < group_index; g_index++)
 void Symmetry::set_atom_map(const Atom* atoms)
 {
     ModuleBase::TITLE("Symmetry", "set_atom_map");
-    if (this->isym_rotiat_.size() > 0) return;
+    if (this->isym_rotiat_.size() == this->nrotk) return;
     this->isym_rotiat_.resize(this->nrotk);
     for (int i = 0; i < this->nrotk; ++i)this->isym_rotiat_[i].resize(this->nat, -1);