mohanchen · Oct 10, 2026 · Oct 7, 2026 · Oct 8, 2026 · Oct 8, 2026 · Oct 9, 2026
diff --git a/source/source_cell/klist.cpp b/source/source_cell/klist.cpp
@@ -7,6 +7,7 @@
 #include "source_base/parallel_reduce.h"
 #include "source_cell/module_symmetry/symmetry.h"
 
+
 void K_Vectors::set(const UnitCell& ucell,
                     const ModuleSymmetry::Symmetry& symm,
                     const std::string& k_file_name,
@@ -546,7 +547,6 @@ void K_Vectors::reduce_by_symmetry(const UnitCell& ucell,
                               nrotkm,
                               kvec_d_ibz,
                               symm.epsilon,
-                              [&symm](double a, double b) { return symm.equal(a, b); },
                               this->kstars);
     }
 
@@ -565,10 +565,64 @@ void K_Vectors::reduce_by_symmetry(const UnitCell& ucell,
     return;
 }
 
-void K_Vectors::set_after_vc(const ModuleBase::Matrix3& G, std::ofstream& ofs_running)
+void K_Vectors::set_after_vc(const ModuleBase::Matrix3& G,
+                             const ModuleSymmetry::Symmetry& symm,
+                             const bool update_symmetry,
+                             std::ofstream& ofs_running)
 {
     ofs_running << "\n SETUP K-POINTS" << std::endl;
 
+    // Every pool holds the complete stars. Their identity entries retain the
+    // global IBZ representatives even when kvec_d contains only local points.
+    // Reassign operation indices without changing the density-matrix ordering.
+    if (update_symmetry && !this->kstars.empty())
+    {
+        const int operation_capacity = 2 * symm.nrotk + symm.nrotk_anti;
+        std::vector<ModuleBase::Matrix3> operations(operation_capacity);
+        for (int operation = 0; operation < symm.nrotk; ++operation)
+        {
+            operations[operation] = symm.kgmatrix[operation];
+        }
+        const int operation_count = KListIO::append_time_reversal_ops(symm, operations, symm.nrotk);
+        typedef std::map<int, ModuleBase::Vector3<double>> KStar;
+        std::vector<KStar> updated_stars(this->kstars.size());
+        for (std::size_t ik = 0; ik < this->kstars.size(); ++ik)
+        {
+            const KStar& star = this->kstars[ik];
+            const KStar::const_iterator identity = star.find(0);
+            if (identity == star.end())
+            {
+                ModuleBase::WARNING_QUIT("K_Vectors::set_after_vc", "Missing k-star IBZ representative.");
+            }
+            const ModuleBase::Vector3<double>& representative = identity->second;
+            for (const KStar::value_type& member : star)
+            {
+                bool matched = false;
+                for (int operation = 0; operation < operation_count; ++operation)
+                {
+                    const ModuleBase::Vector3<double> rotated = member.second * operations[operation];
+                    // Representatives may retain the boundary convention of an older tolerance.
+                    if (ModuleCell::kpoints_equivalent(rotated, representative, symm.epsilon))
+                    {
+                        updated_stars[ik].emplace(operation, member.second);
+                        matched = true;
+                        break;
+                    }
+                }
+                if (!matched)
+                {
+                    ModuleBase::WARNING_QUIT("K_Vectors::set_after_vc",
+                                            "New cell symmetry cannot preserve the existing k-star.");
+                }
+            }
+            if (updated_stars[ik].size() != star.size())
+            {
+                ModuleBase::WARNING_QUIT("K_Vectors::set_after_vc", "New cell symmetry changed the k-star size.");
+            }
+        }
+        this->kstars.swap(updated_stars);
+    }
+
     // set cartesian k vectors.
     this->kvec_d2c(G);
 

diff --git a/source/source_cell/klist.h b/source/source_cell/klist.h
@@ -167,9 +167,17 @@ class K_Vectors : public ModuleCell::ReciprocalGrid
      * as up to date. The spin multiplicity is not touched: it was fixed by
      * set() and never changes during a run.
      *
+     * Reindex existing k-stars using the newly analyzed symmetry operations,
+     * preserving IBZ representatives, weights, and pool/spin ownership.
+     *
      * @param G The new reciprocal lattice matrix.
+     * @param symm The symmetry analyzed for the new cell.
+     * @param update_symmetry Whether crystal symmetry is enabled.
      */
-    void set_after_vc(const ModuleBase::Matrix3& G, std::ofstream& ofs_running);
+    void set_after_vc(const ModuleBase::Matrix3& G,
+                      const ModuleSymmetry::Symmetry& symm,
+                      const bool update_symmetry,
+                      std::ofstream& ofs_running);
 
     //====================================================================
     // Test seam.

diff --git a/source/source_cell/klist_io.cpp b/source/source_cell/klist_io.cpp
@@ -176,7 +176,6 @@ void build_kstars(const std::vector<ModuleBase::Vector3<double>>& kvec_d,
                   const int nrotkm,
                   const std::vector<ModuleBase::Vector3<double>>& kvec_d_ibz,
                   const double epsilon,
-                  const std::function<bool(double, double)>& equal,
                   std::vector<std::map<int, ModuleBase::Vector3<double>>>& kstars)
 {
     const int nkstot = static_cast<int>(kvec_d.size());
@@ -194,8 +193,7 @@ void build_kstars(const std::vector<ModuleBase::Vector3<double>>& kvec_d,
             ModuleCell::restrict_kpt(kvec_rot, epsilon);
             for (int k = 0; k < nkstot_ibz; ++k)
             {
-                if (equal(kvec_rot.x, kvec_d_ibz[k].x) && equal(kvec_rot.y, kvec_d_ibz[k].y)
-                    && equal(kvec_rot.z, kvec_d_ibz[k].z))
+                if (ModuleCell::kpoints_equivalent(kvec_rot, kvec_d_ibz[k], epsilon))
                 {
                     isym = j;
                     exist_number = k;

diff --git a/source/source_cell/klist_io.h b/source/source_cell/klist_io.h
@@ -5,7 +5,6 @@
 #include "source_base/vector3.h"
 
 #include <fstream>
-#include <functional>
 #include <map>
 #include <string>
 #include <vector>
@@ -61,15 +60,14 @@ LineK interp_line(std::ifstream& ifk, int nks_special);
 
 /// Build the EXX k-stars: for every k point, find the symmetry operation
 /// (index into `kgmatrix`) that rotates it onto an irreducible k point, and
-/// group k points by that IBZ representative. `equal` compares two doubles
-/// with the symmetry precision; `epsilon` is the k-restriction tolerance.
+/// group k points by that IBZ representative using periodic equivalence.
+/// `epsilon` is the symmetry and k-restriction tolerance.
 /// this-free so the heavy triple loop is isolated and testable.
 void build_kstars(const std::vector<ModuleBase::Vector3<double>>& kvec_d,
                   const std::vector<ModuleBase::Matrix3>& kgmatrix,
                   int nrotkm,
                   const std::vector<ModuleBase::Vector3<double>>& kvec_d_ibz,
                   double epsilon,
-                  const std::function<bool(double, double)>& equal,
                   std::vector<std::map<int, ModuleBase::Vector3<double>>>& kstars);
 
 /// Append the time-reversal-related k-point symmetry operations into

diff --git a/source/source_cell/module_symmetry/irreducible_sector.cpp b/source/source_cell/module_symmetry/irreducible_sector.cpp
@@ -68,33 +68,14 @@ namespace ModuleSymmetry
         const ModuleBase::Matrix3& gmatd, const TCdouble gtransd,
         const TCdouble& posd_a1, const TCdouble& posd_a2)const
     {
-        // auto restrict_center = [&symm](const TCdouble& v) -> TCdouble {
-        //     // in [-0.5, 0.5)
-        //     TCdouble vr;
-        //     vr.x = fmod(v.x + 100.5 + 0.5 * symm.epsilon, 1) - 0.5 - 0.5 * symm.epsilon;
-        //     vr.y = fmod(v.y + 100.5 + 0.5 * symm.epsilon, 1) - 0.5 - 0.5 * symm.epsilon;
-        //     vr.z = fmod(v.z + 100.5 + 0.5 * symm.epsilon, 1) - 0.5 - 0.5 * symm.epsilon;
-        //     if (std::abs(vr.x) < symm.epsilon) vr.x = 0.0;
-        //     if (std::abs(vr.y) < symm.epsilon) vr.y = 0.0;
-        //     if (std::abs(vr.z) < symm.epsilon) vr.z = 0.0;
-        //     return vr;
-        //     };
-        auto restrict_center = [&symm](const TCdouble& v) -> TCdouble {
-            // in [0,1)
-            TCdouble vr;
-            vr.x = fmod(v.x + 100 + symm.epsilon, 1) - symm.epsilon;
-            vr.y = fmod(v.y + 100 + symm.epsilon, 1) - symm.epsilon;
-            vr.z = fmod(v.z + 100 + symm.epsilon, 1) - symm.epsilon;
-            if (std::abs(vr.x) < symm.epsilon) vr.x = 0.0;
-            if (std::abs(vr.y) < symm.epsilon) vr.y = 0.0;
-            if (std::abs(vr.z) < symm.epsilon) vr.z = 0.0;
-            return vr;
-            };
+        // Keep the representatives used by the orbital integrals. Wrapping an
+        // atom near a cell boundary without retaining its integer shift changes
+        // the Bloch phase and the real-space atom-pair mapping.
         auto check_integer = [&symm](const double x) -> void {
             assert(symm.equal(x, std::round(x)));
             };
-        TCdouble rotpos1 = restrict_center(posd_a1) * gmatd + restrict_center(gtransd);  // row vector
-        TCdouble return_lattice_double = rotpos1 - restrict_center(posd_a2);
+        TCdouble rotpos1 = posd_a1 * gmatd + gtransd;  // row vector
+        TCdouble return_lattice_double = rotpos1 - posd_a2;
 #ifdef __DEBUG
         check_integer(return_lattice_double.x);
         check_integer(return_lattice_double.y);

diff --git a/source/source_cell/module_symmetry/irreducible_sector.h b/source/source_cell/module_symmetry/irreducible_sector.h
@@ -69,7 +69,7 @@ namespace ModuleSymmetry
 
         /// Perfoming {R|t} to atom position r in the R=0 lattice, we get Rr+t, which may get out of R=0 lattice,
         /// whose image in R=0 lattice is r'=Rr+t-O. This function is to get O for each atom and each symmetry operation.
-        /// the range of direct position is [-0.5, 0.5).
+        /// Preserve the input coordinate representatives; do not wrap away their integer lattice shifts.
         TCdouble get_return_lattice(const Symmetry& symm,
             const ModuleBase::Matrix3& gmatd, const TCdouble gtransd,
             const TCdouble& posd_a1, const TCdouble& posd_a2)const;

diff --git a/source/source_cell/module_symmetry/symm_rotation_k.cpp b/source/source_cell/module_symmetry/symm_rotation_k.cpp
@@ -10,6 +10,17 @@
 
 namespace ModuleSymmetry
 {
+    void Symmetry_rotation_k::reset_symmetry()
+    {
+        this->irs_ = Irreducible_Sector();
+        this->Ms_.clear();
+        this->little_groups_.clear();
+        this->spin_U_.clear();
+        this->rotmat_Slm_.clear();
+        // Invalidate the derived RI::Tensor rotation cache as well.
+        ++this->rotmat_Slm_version_;
+    }
+
     std::vector<TC> Symmetry_rotation_k::get_bvk_cells(const TC& period)
     {
         std::vector<TC> cells;

diff --git a/source/source_cell/module_symmetry/symm_rotation_k.h b/source/source_cell/module_symmetry/symm_rotation_k.h
@@ -24,6 +24,9 @@ namespace ModuleSymmetry
         Symmetry_rotation_k() {};
         virtual ~Symmetry_rotation_k() {};
 
+        /// Discard geometry-dependent caches before rebuilding symmetry for a new cell.
+        void reset_symmetry();
+
         //--------------------------------------------------------------------------------
         // getters
         const std::map<Tap, std::set<TC>>& get_irreducible_sector()const { return this->irs_.get_irreducible_sector(); }

diff --git a/source/source_cell/module_symmetry/symmetry.cpp b/source/source_cell/module_symmetry/symmetry.cpp
@@ -11,12 +11,11 @@ bool Symmetry::pricell_loop = true;
 void Symmetry::set_atom_map(const Atom* atoms)
 {
     ModuleBase::TITLE("Symmetry", "set_atom_map");
-    if (this->isym_rotiat_.size() == this->nrotk) {
-        return;
-    }
+    // A new analysis can reorder operations without changing their number.
+    // Rebuild their atom mappings from the current geometry every time.
     this->isym_rotiat_.resize(this->nrotk);
     for (int i = 0; i < this->nrotk; ++i) {
-        this->isym_rotiat_[i].resize(this->nat, -1);
+        this->isym_rotiat_[i].assign(this->nat, -1);
     }
 
     double* pos = this->newpos;

diff --git a/source/source_cell/module_symmetry/test/symm_test_analysis.cpp b/source/source_cell/module_symmetry/test/symm_test_analysis.cpp
@@ -39,6 +39,64 @@ SepPot::~SepPot(){}
 Sep_Cell::Sep_Cell() noexcept {}
 Sep_Cell::~Sep_Cell() noexcept {}
 
+TEST_F(SymmetryTest, AtomMapRefreshWithUnchangedOperationCount)
+{
+    stru_ diamond;
+    diamond.cell = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
+    diamond.coordtype = "D";
+    atomtype_ carbon;
+    carbon.atomname = "C";
+    carbon.coordinate = {{0.0, 0.0, 0.0}, {0.0, 0.5, 0.5},
+                         {0.5, 0.0, 0.5}, {0.5, 0.5, 0.0},
+                         {0.25, 0.25, 0.25}, {0.25, 0.75, 0.75},
+                         {0.75, 0.25, 0.75}, {0.75, 0.75, 0.25}};
+    diamond.all_type.push_back(carbon);
+    construct_ucell(diamond);
+    ModuleSymmetry::Symmetry symm;
+    const int representation[2] = {0, 0};
+    const std::string calculation = "cell-relax";
+    symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running,
+                  1e-6, 1, calculation, representation);
+    const int original_count = symm.nrotk;
+    EXPECT_EQ(original_count, 48);
+
+    // Keep the same crystal symmetry while changing atom representatives.
+    std::swap(ucell.atoms[0].taud[0], ucell.atoms[0].taud[1]);
+    ucell.latvec *= 0.98;
+    ucell.a1 *= 0.98;
+    ucell.a2 *= 0.98;
+    ucell.a3 *= 0.98;
+    for (int atom = 0; atom < ucell.nat; ++atom)
+    {
+        ucell.atoms[0].tau[atom] = ucell.atoms[0].taud[atom] * ucell.latvec;
+    }
+    symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running,
+                  1e-6, 1, calculation, representation);
+    EXPECT_EQ(symm.nrotk, original_count);
+    for (int operation = 0; operation < symm.nrotk; ++operation)
+    {
+        for (int atom = 0; atom < ucell.nat; ++atom)
+        {
+            const int mapped = symm.get_rotated_atom(operation, atom);
+            EXPECT_GE(mapped, 0);
+            EXPECT_LT(mapped, ucell.nat);
+            if (mapped < 0 || mapped >= ucell.nat)
+            {
+                continue;
+            }
+            const auto rotated = ucell.atoms[0].taud[atom] * symm.gmatrix[operation];
+            const auto residual = rotated + symm.gtrans[operation] - ucell.atoms[0].taud[mapped];
+            const double lattice_x = std::round(residual.x);
+            const double lattice_y = std::round(residual.y);
+            const double lattice_z = std::round(residual.z);
+            EXPECT_NEAR(residual.x, lattice_x, 1e-8);
+            EXPECT_NEAR(residual.y, lattice_y, 1e-8);
+            EXPECT_NEAR(residual.z, lattice_z, 1e-8);
+        }
+    }
+    ClearUcell();
+}
+
 TEST_F(SymmetryTest, AnalySys)
 {
     for (int stru = 0; stru < stru_lib.size(); stru++)

diff --git a/source/source_cell/reciprocal_grid.cpp b/source/source_cell/reciprocal_grid.cpp
@@ -6,6 +6,8 @@
  */
 #include "reciprocal_grid.h"
 
+#include <cmath>
+
 #include "source_cell/unitcell.h"
 #include "source_cell/module_symmetry/symmetry.h"
 #include "source_base/global_function.h"
@@ -18,6 +20,19 @@
 namespace ModuleCell
 {
 
+bool kpoints_equivalent(const ModuleBase::Vector3<double>& first,
+                        const ModuleBase::Vector3<double>& second,
+                        const double epsilon)
+{
+    ModuleBase::Vector3<double> difference = first - second;
+    difference.x -= std::round(difference.x);
+    difference.y -= std::round(difference.y);
+    difference.z -= std::round(difference.z);
+    return std::abs(difference.x) < epsilon
+           && std::abs(difference.y) < epsilon
+           && std::abs(difference.z) < epsilon;
+}
+
 void restrict_kpt(ModuleBase::Vector3<double>& kvec, double epsilon)
 {
     // fold into (-0.5, 0.5]; the epsilon shift keeps points sitting on the
@@ -329,7 +344,6 @@ void ReciprocalGrid::reduce_ibz(const ModuleBase::Matrix3* rot_ops,
                                 std::vector<int>& ibz_index,
                                 std::vector<int>& ibz2bz)
 {
-    auto equal = [epsilon](double m, double n) { return fabs(m - n) < epsilon; };
 
     // direct coordinates of points in the k-lattice
     std::vector<ModuleBase::Vector3<double>> kvec_d_k(this->nkstot);
@@ -382,15 +396,14 @@ void ReciprocalGrid::reduce_ibz(const ModuleBase::Matrix3* rot_ops,
                     kvec_rot_k = kvec_rot_k * k_lattice * G.Inverse();   // convert to recip lattice
                     restrict_kpt(kvec_rot_k, epsilon);
 
-                    assert(equal(kvec_rot.x, kvec_rot_k.x));
-                    assert(equal(kvec_rot.y, kvec_rot_k.y));
-                    assert(equal(kvec_rot.z, kvec_rot_k.z));
+                    const bool consistent = kpoints_equivalent(kvec_rot, kvec_rot_k, epsilon);
+                    assert(consistent);
+                    (void)consistent;
                     kvec_rot_k = kvec_rot_k * G * k_lattice.Inverse(); // convert back to k-lattice
                 }
                 for (int k = 0; k < nkstot_ibz; ++k)
                 {
-                    if (equal(kvec_rot.x, kvec_d_ibz[k].x) && equal(kvec_rot.y, kvec_d_ibz[k].y)
-                        && equal(kvec_rot.z, kvec_d_ibz[k].z))
+                    if (kpoints_equivalent(kvec_rot, kvec_d_ibz[k], epsilon))
                     {
                         already_exist = true;
                         // find another ibz point,

diff --git a/source/source_cell/reciprocal_grid.h b/source/source_cell/reciprocal_grid.h
@@ -36,6 +36,11 @@ namespace ModuleCell
  */
 void restrict_kpt(ModuleBase::Vector3<double>& kvec, double epsilon);
 
+/// Compare reciprocal fractional coordinates modulo integer reciprocal vectors.
+bool kpoints_equivalent(const ModuleBase::Vector3<double>& first,
+                        const ModuleBase::Vector3<double>& second,
+                        double epsilon);
+
 /**
  * @brief Abstract base class shared by K_Vectors (electrons) and QList (phonons).
  *