Fix #8078: Refactor: remove the UnitCell Statistics layer and flattening reference members

AGENTS.md

@@ -115,11 +115,6 @@ rules. Read the complete governance document before making or reviewing changes:
 - Do not relax existing tests or references merely to make a failure pass.
   Update references only when the intended behavior changed and the PR explains
   why.
-- When mocking `UnitCell` in a test fixture, do not `delete[] iat2it` or
-  `iat2ia` in `TearDown`: they are owned by `UnitCell`'s internal `Statistics`
-  member, whose destructor releases them. Deleting them again causes a double
-  free. Mirror the ownership pattern of existing fixtures such as
-  `source/source_lcao/module_dftu/test/dftu_lcao_test.cpp`.
 
 ## Review And Exception Flow
 

source/Makefile.Objects

@@ -632,7 +632,7 @@ OBJS_SYMMETRY=symm_other.o\
     symm_rho.o\
     little_group.o\
     symmetry.o\
-    irreducible_sector.o\
+    irred_sec.o\
     irred_sec_bvk.o\
     symm_rotation_k.o\
 

source/source_cell/cal_wfc.cpp

@@ -53,10 +53,8 @@ namespace unitcell
         // This assertion ensures consistency between the two calculation paths.
         assert(nlocal_tmp > 0);
         assert(nlocal_tmp == nlocal);
-        delete[] ucell.iwt2iat;
-        delete[] ucell.iwt2iw;
-        ucell.iwt2iat = new int[nlocal_tmp];
-        ucell.iwt2iw = new int[nlocal_tmp];
+        ucell.iwt2iat.resize(nlocal_tmp);
+        ucell.iwt2iw.resize(nlocal_tmp);
 
         ucell.itia2iat.create(ntype, ucell.namax);
         ucell.set_iat2iwt(npol);

source/source_cell/module_symmetry/CMakeLists.txt

@@ -17,7 +17,7 @@ add_library(
     symmetry.cpp
     symm_rot_spin.cpp
     little_group.cpp
-    irreducible_sector.cpp
+    irred_sec.cpp
     irred_sec_bvk.cpp
     symm_rotation_k.cpp
 )

source/source_cell/module_symmetry/irreducible_sector.cpp → source/source_cell/module_symmetry/irred_sec.cpp

@@ -1,4 +1,11 @@
-#include "source_cell/module_symmetry/irreducible_sector.h"
+#include "source_cell/module_symmetry/irred_sec.h"
+#include "source_base/global_variable.h"
+
+#include <cassert>
+#include <cmath>
+#include <fstream>
+#include <iostream>
+
 namespace ModuleSymmetry
 {
     // Raw-index dispatch shared by the real-space sector helpers, matching the convention used
@@ -16,7 +23,7 @@ namespace ModuleSymmetry
                                    : symm.get_rotated_atom_anti(isym - symm.nrotk, iat);
     }
 
-    TC Irreducible_Sector::rotate_R(const Symmetry& symm,
+    TC Irred_Sector::rotate_R(const Symmetry& symm,
         const int isym, const int iat1, const int iat2, const TC& R, const char gauge) const
     {
         auto round2int = [symm](const double x) -> int { return x > 0 ? static_cast<int>(x + symm.epsilon) : static_cast<int>(x - symm.epsilon); };
@@ -28,27 +35,14 @@ namespace ModuleSymmetry
             : R_double * gmat + this->return_lattice_[iat2][isym] - this->return_lattice_[iat1][isym];
         return { round2int(Rrot_double.x), round2int(Rrot_double.y), round2int(Rrot_double.z) };
     }
-    TapR Irreducible_Sector::rotate_apR_by_formula(const Symmetry& symm,
+    TapR Irred_Sector::rotate_apR_by_formula(const Symmetry& symm,
         const int isym, const TapR& apR, const char gauge) const
     {
         const Tap& aprot = { sector_rotated_atom(symm, isym, apR.first.first),
                              sector_rotated_atom(symm, isym, apR.first.second) };
         return { aprot, this->rotate_R(symm, isym, apR.first.first, apR.first.second, apR.second, gauge) };
     }
 
-    TCdouble Irreducible_Sector::get_aRb_direct(const Atom* atoms, const Statistics& st,
-        const int iat1, const int iat2, const TCdouble& R, const char gauge)const
-    {
-        return TCdouble(atoms[st.iat2it[iat1]].taud[st.iat2ia[iat1]] - atoms[st.iat2it[iat2]].taud[st.iat2ia[iat2]]) + (gauge == 'L' ? R : TCdouble(-R));
-    }
-
-    TCdouble Irreducible_Sector::get_aRb_direct(const Atom* atoms, const Statistics& st,
-        const int iat1, const int iat2, const TC& R, const char gauge) const
-    {
-        const TCdouble R_double(static_cast<double>(R[0]), static_cast<double>(R[1]), static_cast<double>(R[2]));
-        return get_aRb_direct(atoms, st, iat1, iat2, R_double);
-    }
-
     inline void output_return_lattice(const std::vector<std::vector<TCdouble>>& return_lattice)
     {
         std::cout << "return lattice:" << std::endl;
@@ -64,7 +58,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).
-    TCdouble Irreducible_Sector::get_return_lattice(const Symmetry& symm,
+    TCdouble Irred_Sector::get_return_lattice(const Symmetry& symm,
         const ModuleBase::Matrix3& gmatd, const TCdouble gtransd,
         const TCdouble& posd_a1, const TCdouble& posd_a2)const
     {
@@ -84,35 +78,36 @@ namespace ModuleSymmetry
         return TCdouble(std::round(return_lattice_double.x), std::round(return_lattice_double.y), std::round(return_lattice_double.z));
     }
 
-    void Irreducible_Sector::cal_return_lattice_all(const Symmetry& symm, const Atom* atoms, const Statistics& st)
+    void Irred_Sector::cal_return_lattice_all(const Symmetry& symm, const Atom* atoms,
+        const int nat, const std::vector<int>& iat2it, const std::vector<int>& iat2ia)
     {
-        ModuleBase::TITLE("Symmetry_rotation", "cal_return_lattice_all");
+        ModuleBase::TITLE("Irred_Sector", "cal_return_lattice_all");
         // Columns [0, nrotk) are the unitary operations; columns [nrotk, nrotk+nrotk_anti) are the
         // spatial parts of the antiunitary elements Theta*g of the Shubnikov group (nspin=4 magnetic),
         // so that Symmetry_rotation can address both with one raw index.
-        this->return_lattice_.resize(st.nat, std::vector<TCdouble>(symm.nrotk + symm.nrotk_anti));
-        for (int iat1 = 0;iat1 < st.nat;++iat1)
+        this->return_lattice_.resize(nat, std::vector<TCdouble>(symm.nrotk + symm.nrotk_anti));
+        for (int iat1 = 0;iat1 < nat;++iat1)
         {
-            int it = st.iat2it[iat1];
-            int ia1 = st.iat2ia[iat1];
+            int it = iat2it[iat1];
+            int ia1 = iat2ia[iat1];
             for (int isym = 0;isym < symm.nrotk;++isym)
             {
                 int iat2 = symm.get_rotated_atom(isym, iat1);
-                int ia2 = st.iat2ia[iat2];
+                int ia2 = iat2ia[iat2];
                 this->return_lattice_[iat1][isym] = get_return_lattice(symm, symm.gmatrix[isym], symm.gtrans[isym], atoms[it].taud[ia1], atoms[it].taud[ia2]);
             }
             for (int j = 0;j < symm.nrotk_anti;++j)
             {
                 int iat2 = symm.get_rotated_atom_anti(j, iat1);
-                int ia2 = st.iat2ia[iat2];
+                int ia2 = iat2ia[iat2];
                 this->return_lattice_[iat1][symm.nrotk + j] = get_return_lattice(symm, symm.gmatrix_anti[j], symm.gtrans_anti[j], atoms[it].taud[ia1], atoms[it].taud[ia2]);
             }
         }
         // test: output return_lattice
         // output_return_lattice(this->return_lattice_);
     }
 
-    void Irreducible_Sector::output_full_map_to_irreducible_sector(const int nat)
+    void Irred_Sector::output_full_map_to_irreducible_sector(const int nat)
     {
         std::cout << "Final map to irreducible sector: " << std::endl;
         for (auto& apR_isym_irapR : this->full_map_to_irreducible_sector_)
@@ -127,7 +122,7 @@ namespace ModuleSymmetry
         }
     }
 
-    void Irreducible_Sector::output_sector_star()
+    void Irred_Sector::output_sector_star()
     {
         std::cout << "Found " << this->sector_stars_.size() << " irreducible sector stars:" << std::endl;
         // for (auto& irs_star : this->sector_stars_)
@@ -150,7 +145,7 @@ namespace ModuleSymmetry
             std::cout << std::endl;
         }
     }
-    void Irreducible_Sector::write_irreducible_sector(const std::string& output_dir)
+    void Irred_Sector::write_irreducible_sector(const std::string& output_dir)
     {
         if(GlobalV::MY_RANK == 0 && !output_dir.empty())
         {
@@ -164,20 +159,23 @@ namespace ModuleSymmetry
         }
     }
 
-    void Irreducible_Sector::find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const Statistics& st, const std::vector<TC>& Rs, const TC& period, const Lattice& lat, const std::string& output_dir)
+    void Irred_Sector::find_irred_sector(const Symmetry& symm, const Atom* atoms,
+        const int nat, const int ntype,
+        const std::vector<int>& iat2it, const std::vector<int>& iat2ia,
+        const std::vector<TC>& Rs, const TC& period, const Lattice& lat, const std::string& output_dir)
     {
         this->full_map_to_irreducible_sector_.clear();
         this->irreducible_sector_.clear();
         this->sector_stars_.clear();
 
-        if (this->return_lattice_.empty()) this->cal_return_lattice_all(symm, atoms, st);
+        if (this->return_lattice_.empty()) this->cal_return_lattice_all(symm, atoms, nat, iat2it, iat2ia);
         // if (this->atompair_stars_.empty()) this->find_irreducible_atom_pairs(symm);
 
         // contruct {atom pair, R} set
         // constider different number of Rs for different atom pairs later.
         std::map<Tap, std::set<TC, len_less_func>> apR_all;
-        for (int iat1 = 0;iat1 < st.nat; iat1++)
-            for (int iat2 = 0; iat2 < st.nat; iat2++)
+        for (int iat1 = 0;iat1 < nat; iat1++)
+            for (int iat2 = 0; iat2 < nat; iat2++)
                 for (auto& R : Rs)
                     apR_all[{iat1, iat2}].insert(R);
 
@@ -200,7 +198,7 @@ namespace ModuleSymmetry
 
         // get symmetry of BvK supercell
         if (this->isymbvk_to_isym_.empty())
-            this->gen_symmetry_BvK(symm, atoms, lat, st, period);
+            this->gen_symm_bvk(symm, atoms, lat, nat, ntype, period);
         assert(!this->isymbvk_to_isym_.empty());
         // std::vector<bool> in_2d_plain;
         // const bool judge_2d = (symm.real_brav == 4);
@@ -255,7 +253,7 @@ namespace ModuleSymmetry
         for (auto& sector : this->sector_stars_)
             total_apR_in_star += sector.second.size();
         assert(total_apR_in_star == this->full_map_to_irreducible_sector_.size());
-        // this->output_full_map_to_irreducible_sector(st.nat);
+        // this->output_full_map_to_irreducible_sector(nat);
         // this->output_sector_star();
         this->write_irreducible_sector(output_dir);
     }

source/source_cell/module_symmetry/irreducible_sector.h → source/source_cell/module_symmetry/irred_sec.h

@@ -1,5 +1,5 @@
-#ifndef ABACUS_SOURCE_CELL_MODULE_SYMMETRY_IRREDUCIBLE_SECTOR_H
-#define ABACUS_SOURCE_CELL_MODULE_SYMMETRY_IRREDUCIBLE_SECTOR_H
+#ifndef ABACUS_SOURCE_CELL_MODULE_SYMMETRY_IRRED_SECTOR_H
+#define ABACUS_SOURCE_CELL_MODULE_SYMMETRY_IRRED_SECTOR_H
 
 #include <vector>
 #include <map>
@@ -16,7 +16,7 @@ namespace ModuleSymmetry
     using TapR = std::pair<Tap, TC>;
     using TCdouble = ModuleBase::Vector3<double>;
 
-    class Irreducible_Sector
+    class Irred_Sector
     {
     public:
         struct ap_less_func
@@ -61,9 +61,11 @@ namespace ModuleSymmetry
         ///                    output directory setting); kept as an explicit argument rather than
         ///                    reading the global config directly so this LibRI-free class has no
         ///                    module_parameter link dependency.
-        void find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const Statistics& st,
+        void find_irred_sector(const Symmetry& symm, const Atom* atoms,
+            const int nat, const int ntype,
+            const std::vector<int>& iat2it, const std::vector<int>& iat2ia,
             const std::vector<TC>& Rs, const TC& period, const Lattice& lat, const std::string& output_dir = "");
-        const std::map<Tap, std::set<TC>>& get_irreducible_sector()const { return this->irreducible_sector_; }
+        const std::map<Tap, std::set<TC>>& get_irred_sector()const { return this->irreducible_sector_; }
         // const std::map<int, std::set<std::pair<int, TC>>> convirt_irreducible_sector() {};
         //--------------------------------------------------------------------------------
 
@@ -86,7 +88,8 @@ namespace ModuleSymmetry
             return this->return_lattice_[iat][isym];
         }
     protected:
-        void cal_return_lattice_all(const Symmetry& symm, const Atom* atoms, const Statistics& st);
+        void cal_return_lattice_all(const Symmetry& symm, const Atom* atoms,
+            const int nat, const std::vector<int>& iat2it, const std::vector<int>& iat2ia);
 
         //--------------------------------------------------------------------------------
         /// The sub functions to find irreducible sector: {abR}
@@ -95,9 +98,6 @@ namespace ModuleSymmetry
         /// gauge='L': R'=R+O_1-O_2; gauge='R': R'=R+O_2-O_1
         TC rotate_R(const Symmetry& symm, const int isym, const int iat1, const int iat2, const TC& R, const char gauge = 'R')const;
         TapR rotate_apR_by_formula(const Symmetry& symm, const int isym, const TapR& iapR, const char gauge = 'R')const;
-        /// gauge='L': tau_a + R - tau_b; gauge='R': tau_a - tau_b - R (direct)
-        TCdouble get_aRb_direct(const Atom* atoms, const Statistics& st, const int iat1, const int iat2, const TC& R, const char gauge = 'R')const;
-        TCdouble get_aRb_direct(const Atom* atoms, const Statistics& st, const int iat1, const int iat2, const TCdouble& R, const char gauge = 'R')const;
 
         ModuleBase::Matrix3 direct_to_cartesian(const ModuleBase::Matrix3& d, const ModuleBase::Matrix3& latvec)const;
 
@@ -115,7 +115,8 @@ namespace ModuleSymmetry
 
         //--------------------------------------------------------------------------------
         /// The sub functions judge special symmetry
-        void gen_symmetry_BvK(const Symmetry& symm, const Atom* atoms, const Lattice& lat, const Statistics& st, const TC bvk_period);
+        void gen_symm_bvk(const Symmetry& symm, const Atom* atoms, const Lattice& lat,
+            const int nat, const int ntype, const TC bvk_period);
         /// whether in 2D plain or not for each symmetry operation
         // std::vector<bool> in_plain(const ModuleSymmetry::Symmetry& symm, const ModuleBase::Matrix3& latvec)const;
         //--------------------------------------------------------------------------------

source/source_cell/module_symmetry/irred_sec_bvk.cpp

@@ -1,7 +1,8 @@
-#include "irreducible_sector.h"
+#include "irred_sec.h"
+#include "source_base/global_variable.h"
 namespace ModuleSymmetry
 {
-    ModuleBase::Matrix3 Irreducible_Sector::direct_to_cartesian(const ModuleBase::Matrix3& d, const ModuleBase::Matrix3& latvec)const
+    ModuleBase::Matrix3 Irred_Sector::direct_to_cartesian(const ModuleBase::Matrix3& d, const ModuleBase::Matrix3& latvec)const
     {
         return latvec.Inverse() * d * latvec;
     }
@@ -43,7 +44,7 @@ namespace ModuleSymmetry
             // Unmatched stays -1. That is legitimate for nspin=4 magnetic: the unit-cell set is the
             // Shubnikov group H (union) A, generally a PROPER subset of the crystallographic group
             // (operations that merely tilt the moment belong to neither), so a BvK operation may
-            // have no counterpart. The consumer in find_irreducible_sector skips negative entries;
+            // have no counterpart. The consumer in find_irred_sector skips negative entries;
             // it must never use one as an index.
         }
         return isymbvk2isym;
@@ -55,9 +56,10 @@ namespace ModuleSymmetry
         int c = a % b;
         return (c == 0) ? b : gcd(b, c);
     }
-    void Irreducible_Sector::gen_symmetry_BvK(const ModuleSymmetry::Symmetry& symm, const Atom* atoms, const Lattice& lat, const Statistics& st, const TC bvk_period)
+    void Irred_Sector::gen_symm_bvk(const ModuleSymmetry::Symmetry& symm, const Atom* atoms, const Lattice& lat,
+        const int nat, const int ntype, const TC bvk_period)
     {
-        ModuleBase::TITLE("Irreducible_Sector", "gen_symmetry_BvK");
+        ModuleBase::TITLE("Irred_Sector", "gen_symm_bvk");
         auto set_matrix3 = [](const ModuleBase::Vector3<double>& a1, const ModuleBase::Vector3<double>& a2, const ModuleBase::Vector3<double>& a3)
             -> ModuleBase::Matrix3 {return ModuleBase::Matrix3(a1.x, a1.y, a1.z, a2.x, a2.y, a2.z, a3.x, a3.y, a3.z);};
         auto set_bvk_same_as_ucell = [&symm, this]()->void
@@ -78,11 +80,11 @@ namespace ModuleSymmetry
         // extern lattice to minimal BvK lattice, and set direct coordinates in min BvK lattice
         int bvk_gcd = gcd(bvk_period[0], gcd(bvk_period[1], bvk_period[2]));
         const TC bvk_min_period = TC({ bvk_period[0] / bvk_gcd, bvk_period[1] / bvk_gcd, bvk_period[2] / bvk_gcd });
-        const int bvk_nat = st.nat * bvk_min_period[0] * bvk_min_period[1] * bvk_min_period[2];
-        std::vector<int> bvk_na(st.ntype);
-        std::vector<int> bvk_istart(st.ntype, 0);
+        const int bvk_nat = nat * bvk_min_period[0] * bvk_min_period[1] * bvk_min_period[2];
+        std::vector<int> bvk_na(ntype);
+        std::vector<int> bvk_istart(ntype, 0);
         int bvk_itmin_start = 0, bvk_itmin_type = 0;
-        for (int it = 0;it < st.ntype;++it)
+        for (int it = 0;it < ntype;++it)
         {
             bvk_na[it] = atoms[it].na * bvk_min_period[0] * bvk_min_period[1] * bvk_min_period[2];
             if (it > 0) { bvk_istart[it] = bvk_istart[it - 1] + bvk_na[it - 1];
@@ -103,7 +105,7 @@ namespace ModuleSymmetry
         s3 = a3 = lat.a3 * static_cast<double>(bvk_min_period[2]);
         ModuleBase::Matrix3 bvk_min_lat = set_matrix3(s1, s2, s3);
         int at = 0;
-        for (int it = 0; it < st.ntype; ++it) {
+        for (int it = 0; it < ntype; ++it) {
             for (int c1 = 0;c1 < bvk_min_period[0];++c1) {
                 for (int c2 = 0;c2 < bvk_min_period[1];++c2) {
                     for (int c3 = 0;c3 < bvk_min_period[2];++c3) {
@@ -184,7 +186,7 @@ namespace ModuleSymmetry
         return;
     }
 
-    // std::vector<bool> Irreducible_Sector::in_plain(const ModuleSymmetry::Symmetry& symm, const ModuleBase::Matrix3& latvec)const
+    // std::vector<bool> Irred_Sector::in_plain(const ModuleSymmetry::Symmetry& symm, const ModuleBase::Matrix3& latvec)const
     // {
     //     // get euler angel of the cartesian gmatrix in optimal lattice
     //     std::vector<ModuleBase::Matrix3> gmatc(symm.nrotk);

source/source_cell/module_symmetry/symm_analysis.cpp

@@ -3,7 +3,12 @@
 
 using namespace ModuleSymmetry;
 
-void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, std::ofstream& ofs_running,
+void Symmetry::analy_sys(const Lattice& lat, Atom* atoms,
+                         const int nat, const int ntype,
+                         const std::vector<int>& iat2it,
+                         const std::vector<int>& iat2ia,
+                         const ModuleBase::IntArray& itia2iat,
+                         std::ofstream& ofs_running,
                          const double symmetry_prec, const int nspin, const std::string& calculation,
                          const int* cal_symm_repr)
 {
@@ -30,9 +35,9 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms,
     // 1. copy data and allocate memory
     // --------------------------------
     // number of total atoms
-    this->nat = st.nat;
+    this->nat = nat;
     // number of atom species
-    this->ntype = st.ntype;
+    this->ntype = ntype;
 
     assert(ntype>0);
 
@@ -116,7 +121,7 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms,
 
             if (!pricell_loop && nspin == 2)
             {
-                this->analyze_magnetic_group(atoms, st, nrot_out, nrotk_out);
+                this->analyze_magnetic_group(atoms, iat2it, iat2ia, itia2iat, nrot_out, nrotk_out);
             }
             else
             {
@@ -288,13 +293,13 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms,
     // the magnetization (pseudovector), so they are not applied in k-reduction / density symmetrization.
     if (nspin == 4)
     {
-        this->analyze_magnetic_group_nspin4(atoms, st, latvec1);
+        this->analyze_magnetic_group_nspin4(atoms, iat2it, iat2ia, latvec1);
     }
 
     // Do this here for debug
     if (calculation == "relax")
     {
-        this->all_mbl = this->is_all_movable(atoms, st);
+        this->all_mbl = this->is_all_movable(atoms, iat2it, iat2ia);
         if (!this->all_mbl)
         {
             std::cout << "WARNING: Symmetry cannot be kept when not all atoms are movable.\n ";