From 1158c992f4c0f2630252111f89b685e558ef9807 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 6 Oct 2026 08:32:12 +0800 Subject: [PATCH 01/15] Refactor: remove the UnitCell Statistics layer and flattening reference members UnitCell stored its atom-count and index-map data in a nested `Statistics` struct and re-exported every field as a reference member (`int& nat = st.nat`, `int*& iat2it = st.iat2it`, ...). This created dual-view aliasing, an invisible raw-pointer ownership rule (fixtures had to `new[]` but never `delete[]`, and several unit tests double-freed as a result), and a latent shallow-copy double-free since the implicit copy constructor remained available while copy assignment was silently deleted by the reference members. Changes: - Repoint the only consumers of the `Statistics` type: module_symmetry and module_exx_symmetry interfaces now take `const UnitCell&` instead of `const Statistics&`; all call sites pass `ucell` (or `*ucell`). - Move the nine fields into UnitCell as direct members under the existing top-level names, so `ucell.nat`-style call sites are unchanged. Delete the `st` member, the reference members, and the `Statistics` struct. - Convert iat2it/iat2ia/iwt2iat/iwt2iw from raw int* to std::vector; migrate the ~100 test/support allocation sites from `new int[]` to resize()/assign(). Production sites (unitcell.cpp, cal_wfc.cpp) and the two force_pw.cpp device-upload call sites use .data(). - Delete UnitCell copy constructor/assignment and default the move operations, so accidental copies fail at compile time. - Update AGENTS.md: replace the obsolete fixture note with the new ownership/vector rule. Verification: - Full build succeeds (no compile errors). - ctest -R "unitcell|CELL_SYMMETRY|symmetry|hcontainer|dftu|operator_lcao| module_dm|pdos|hsr|restart_exx|gint|relax|esolver": 33/33 passed. - agent_governance_check: no issues. No INPUT parameter or user-facing behavior change; docs/parameters.yaml unaffected. The Lattice lat reference-member block is intentionally out of scope and can follow the same recipe in a later PR. --- AGENTS.md | 10 ++--- source/source_cell/cal_wfc.cpp | 6 +-- .../module_symmetry/irred_sec_bvk.cpp | 13 ++++--- .../module_symmetry/irreducible_sector.cpp | 35 ++++++++--------- .../module_symmetry/irreducible_sector.h | 10 ++--- .../module_symmetry/symm_analysis.cpp | 13 ++++--- .../module_symmetry/symm_magnetic.cpp | 25 ++++++------ .../module_symmetry/symm_rotation_k.cpp | 13 ++++--- .../module_symmetry/symm_rotation_k.h | 6 +-- .../source_cell/module_symmetry/symmetry.cpp | 9 +++-- source/source_cell/module_symmetry/symmetry.h | 10 +++-- .../test/symm_test_analysis.cpp | 4 +- .../module_symmetry/test/symm_test_symtrz.cpp | 4 +- source/source_cell/test/klist_test.cpp | 6 +-- source/source_cell/test/klist_test_para.cpp | 4 +- source/source_cell/test/little_group_test.cpp | 2 +- source/source_cell/test/qlist_test.cpp | 14 +++---- source/source_cell/test/unitcell_test.cpp | 8 ++-- source/source_cell/unitcell.cpp | 6 +-- source/source_cell/unitcell.h | 28 +++++++++----- source/source_cell/unitcell_data.h | 38 ------------------- source/source_esolver/esolver_fp.cpp | 4 +- source/source_esolver/esolver_gets.cpp | 2 +- .../source_esolver/esolver_lr_lcao_tddft.cpp | 2 +- .../source_esolver/test/esolver_dp_test.cpp | 4 +- .../module_dm/unittests/prepare_unitcell.h | 12 ++---- .../module_dm/unittests/test_cal_dm_r.cpp | 4 +- .../unittests/test_dm_constructor.cpp | 4 +- .../module_dm/unittests/test_dm_from_psi.cpp | 4 +- .../module_dm/unittests/test_dm_r_init.cpp | 4 +- source/source_estate/test/test_occ_mixer.cpp | 6 +-- .../module_gint/test/test_gint_common.cpp | 4 +- .../module_hcontainer/test/prepare_unitcell.h | 12 ++---- .../test/test_func_folding.cpp | 2 +- .../test/test_hcontainer.cpp | 2 +- .../test/test_hcontainer_complex.cpp | 2 +- .../test/test_hcontainer_output.cpp | 8 ++-- .../test/test_hcontainer_time.cpp | 2 +- .../module_hcontainer/test/test_transfer.cpp | 4 +- .../unittests/test_write_pdos_text.cpp | 6 +-- .../module_hs/unittests/test_hsr_writer.cpp | 2 +- .../source_io/test/restart_exx_csr_test.cpp | 4 +- .../source_lcao/module_dftu/dftu_nao_occ.cpp | 2 +- .../source_lcao/module_dftu/dftu_nao_op.cpp | 2 +- .../unittests/test_dftu_nao_fs_accum.cpp | 3 +- .../unittests/test_dftu_nao_ijr.cpp | 12 ++---- .../unittests/test_dftu_nao_op.cpp | 4 +- .../ri_benchmark/test/ri_benchmark_test.cpp | 2 +- .../test/test_ekinetic.cpp | 4 +- .../test/test_ekinetic_serial.cpp | 8 ++-- .../test/test_nonlocal.cpp | 4 +- .../test/test_overlap.cpp | 4 +- .../test/test_overlap_cd.cpp | 4 +- .../test/test_overlap_serial.cpp | 8 ++-- .../test/test_t_nl_cd.cpp | 4 +- source/source_lcao/module_rdmft/rdmft.cpp | 2 +- source/source_lcao/module_ri/exx_lri.hpp | 4 +- .../module_ri/exx_lri_interface.hpp | 2 +- .../module_exx_symmetry/symm_rotation.h | 15 ++++---- .../module_exx_symmetry/symm_rotation_r.hpp | 20 +++++----- .../symm_rotation_r_hcontainer.hpp | 22 +++++------ source/source_lcao/module_ri/rpa_lri.hpp | 2 +- .../test/test_init_dm_from_file.cpp | 4 +- .../test_output_hcontainer_consistency.cpp | 4 +- .../module_dfpt/test/dfpt_pw_data_test.cpp | 2 +- .../module_dfpt/test/dfpt_pw_run_test.cpp | 2 +- .../test_serial/dfpt_q0_serial_test.cpp | 6 +-- .../test_serial/dfpt_serial_fixture.cpp | 14 ++----- source/source_pw/module_pwdft/force_pw.cpp | 4 +- .../module_pwdft/test/dftu_base_test.cpp | 4 +- source/source_relax/unittests/for_test.h | 5 --- .../unittests/test_ions_move_bfgs2.cpp | 16 ++++---- .../unittests/test_ions_move_lbfgs.cpp | 4 -- .../unittests/test_relax_nsync.cpp | 15 ++------ .../unittests/test_relax_sync.cpp | 29 ++++++-------- 75 files changed, 261 insertions(+), 334 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 336ce0a4358..2e1fab066c3 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -115,11 +115,11 @@ 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`. +- `UnitCell` is non-copyable (copy constructor/assignment deleted) and movable; + its index maps (`iat2it`, `iat2ia`, `iwt2iat`, `iwt2iw`) are + `std::vector` members owned directly by `UnitCell`. Size them with + `resize()`/`assign()`; never `new[]`/`delete[]` them. Do not reintroduce a + nested data struct re-exported through reference members. ## Review And Exception Flow diff --git a/source/source_cell/cal_wfc.cpp b/source/source_cell/cal_wfc.cpp index b6307e5fc9e..5044ea33a47 100644 --- a/source/source_cell/cal_wfc.cpp +++ b/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); diff --git a/source/source_cell/module_symmetry/irred_sec_bvk.cpp b/source/source_cell/module_symmetry/irred_sec_bvk.cpp index 5d506cab434..115a8bbb3e6 100644 --- a/source/source_cell/module_symmetry/irred_sec_bvk.cpp +++ b/source/source_cell/module_symmetry/irred_sec_bvk.cpp @@ -1,4 +1,5 @@ #include "irreducible_sector.h" +#include "source_cell/unitcell.h" namespace ModuleSymmetry { ModuleBase::Matrix3 Irreducible_Sector::direct_to_cartesian(const ModuleBase::Matrix3& d, const ModuleBase::Matrix3& latvec)const @@ -55,7 +56,7 @@ 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 Irreducible_Sector::gen_symmetry_BvK(const ModuleSymmetry::Symmetry& symm, const Atom* atoms, const Lattice& lat, const UnitCell& ucell, const TC bvk_period) { ModuleBase::TITLE("Irreducible_Sector", "gen_symmetry_BvK"); auto set_matrix3 = [](const ModuleBase::Vector3& a1, const ModuleBase::Vector3& a2, const ModuleBase::Vector3& a3) @@ -78,11 +79,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 bvk_na(st.ntype); - std::vector bvk_istart(st.ntype, 0); + const int bvk_nat = ucell.nat * bvk_min_period[0] * bvk_min_period[1] * bvk_min_period[2]; + std::vector bvk_na(ucell.ntype); + std::vector bvk_istart(ucell.ntype, 0); int bvk_itmin_start = 0, bvk_itmin_type = 0; - for (int it = 0;it < st.ntype;++it) + for (int it = 0;it < ucell.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 +104,7 @@ namespace ModuleSymmetry s3 = a3 = lat.a3 * static_cast(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 < ucell.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) { diff --git a/source/source_cell/module_symmetry/irreducible_sector.cpp b/source/source_cell/module_symmetry/irreducible_sector.cpp index c6521c7ab52..49b6897eb19 100644 --- a/source/source_cell/module_symmetry/irreducible_sector.cpp +++ b/source/source_cell/module_symmetry/irreducible_sector.cpp @@ -1,4 +1,5 @@ #include "source_cell/module_symmetry/irreducible_sector.h" +#include "source_cell/unitcell.h" namespace ModuleSymmetry { // Raw-index dispatch shared by the real-space sector helpers, matching the convention used @@ -36,17 +37,17 @@ namespace ModuleSymmetry 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, + TCdouble Irreducible_Sector::get_aRb_direct(const Atom* atoms, const UnitCell& ucell, 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)); + return TCdouble(atoms[ucell.iat2it[iat1]].taud[ucell.iat2ia[iat1]] - atoms[ucell.iat2it[iat2]].taud[ucell.iat2ia[iat2]]) + (gauge == 'L' ? R : TCdouble(-R)); } - TCdouble Irreducible_Sector::get_aRb_direct(const Atom* atoms, const Statistics& st, + TCdouble Irreducible_Sector::get_aRb_direct(const Atom* atoms, const UnitCell& ucell, const int iat1, const int iat2, const TC& R, const char gauge) const { const TCdouble R_double(static_cast(R[0]), static_cast(R[1]), static_cast(R[2])); - return get_aRb_direct(atoms, st, iat1, iat2, R_double); + return get_aRb_direct(atoms, ucell, iat1, iat2, R_double); } inline void output_return_lattice(const std::vector>& return_lattice) @@ -103,27 +104,27 @@ 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 Irreducible_Sector::cal_return_lattice_all(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell) { ModuleBase::TITLE("Symmetry_rotation", "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(symm.nrotk + symm.nrotk_anti)); - for (int iat1 = 0;iat1 < st.nat;++iat1) + this->return_lattice_.resize(ucell.nat, std::vector(symm.nrotk + symm.nrotk_anti)); + for (int iat1 = 0;iat1 < ucell.nat;++iat1) { - int it = st.iat2it[iat1]; - int ia1 = st.iat2ia[iat1]; + int it = ucell.iat2it[iat1]; + int ia1 = ucell.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 = ucell.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 = ucell.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]); } } @@ -183,20 +184,20 @@ namespace ModuleSymmetry } } - void Irreducible_Sector::find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const Statistics& st, const std::vector& Rs, const TC& period, const Lattice& lat, const std::string& output_dir) + void Irreducible_Sector::find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const std::vector& 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, ucell); // 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> apR_all; - for (int iat1 = 0;iat1 < st.nat; iat1++) - for (int iat2 = 0; iat2 < st.nat; iat2++) + for (int iat1 = 0;iat1 < ucell.nat; iat1++) + for (int iat2 = 0; iat2 < ucell.nat; iat2++) for (auto& R : Rs) apR_all[{iat1, iat2}].insert(R); @@ -219,7 +220,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_symmetry_BvK(symm, atoms, lat, ucell, period); assert(!this->isymbvk_to_isym_.empty()); // std::vector in_2d_plain; // const bool judge_2d = (symm.real_brav == 4); @@ -274,7 +275,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(ucell.nat); // this->output_sector_star(); this->write_irreducible_sector(output_dir); } diff --git a/source/source_cell/module_symmetry/irreducible_sector.h b/source/source_cell/module_symmetry/irreducible_sector.h index f35be6a9d2f..e01e1e7fa88 100644 --- a/source/source_cell/module_symmetry/irreducible_sector.h +++ b/source/source_cell/module_symmetry/irreducible_sector.h @@ -61,7 +61,7 @@ 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_irreducible_sector(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const std::vector& Rs, const TC& period, const Lattice& lat, const std::string& output_dir = ""); const std::map>& get_irreducible_sector()const { return this->irreducible_sector_; } // const std::map>> convirt_irreducible_sector() {}; @@ -86,7 +86,7 @@ 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 UnitCell& ucell); //-------------------------------------------------------------------------------- /// The sub functions to find irreducible sector: {abR} @@ -96,8 +96,8 @@ namespace ModuleSymmetry 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; + TCdouble get_aRb_direct(const Atom* atoms, const UnitCell& ucell, const int iat1, const int iat2, const TC& R, const char gauge = 'R')const; + TCdouble get_aRb_direct(const Atom* atoms, const UnitCell& ucell, 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,7 @@ 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_symmetry_BvK(const Symmetry& symm, const Atom* atoms, const Lattice& lat, const UnitCell& ucell, const TC bvk_period); /// whether in 2D plain or not for each symmetry operation // std::vector in_plain(const ModuleSymmetry::Symmetry& symm, const ModuleBase::Matrix3& latvec)const; //-------------------------------------------------------------------------------- diff --git a/source/source_cell/module_symmetry/symm_analysis.cpp b/source/source_cell/module_symmetry/symm_analysis.cpp index ae5f4dd36a9..8c8d16c0fdb 100644 --- a/source/source_cell/module_symmetry/symm_analysis.cpp +++ b/source/source_cell/module_symmetry/symm_analysis.cpp @@ -1,9 +1,10 @@ #include "symmetry.h" +#include "source_cell/unitcell.h" #include "source_base/output.h" 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, const UnitCell& ucell, Atom* atoms, std::ofstream& ofs_running, const double symmetry_prec, const int nspin, const std::string& calculation, const int* cal_symm_repr) { @@ -30,9 +31,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 = ucell.nat; // number of atom species - this->ntype = st.ntype; + this->ntype = ucell.ntype; assert(ntype>0); @@ -122,7 +123,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, ucell, nrot_out, nrotk_out); } else { @@ -294,13 +295,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, ucell, 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, ucell); if (!this->all_mbl) { std::cout << "WARNING: Symmetry cannot be kept when not all atoms are movable.\n "; diff --git a/source/source_cell/module_symmetry/symm_magnetic.cpp b/source/source_cell/module_symmetry/symm_magnetic.cpp index c7439684de4..60f314616ee 100644 --- a/source/source_cell/module_symmetry/symm_magnetic.cpp +++ b/source/source_cell/module_symmetry/symm_magnetic.cpp @@ -3,11 +3,12 @@ using namespace ModuleSymmetry; #include "symm_rot_spin.h" #include "source_base/global_variable.h" +#include "source_cell/unitcell.h" #include #include -void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, int& nrot_out, int& nrotk_out) +void Symmetry::analyze_magnetic_group(const Atom* atoms, const UnitCell& ucell, int& nrot_out, int& nrotk_out) { // 1. classify atoms with different magmom // (use symmetry_prec to judge if two magmoms are the same) @@ -20,18 +21,18 @@ void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, i for (auto& mt : mag_type_atoms) { const int mag_iat = *mt.begin(); - const int mag_it = st.iat2it[mag_iat]; - const int mag_ia = st.iat2ia[mag_iat]; + const int mag_it = ucell.iat2it[mag_iat]; + const int mag_ia = ucell.iat2ia[mag_iat]; if (it == mag_it && this->equal(atoms[it].mag[ia], atoms[mag_it].mag[mag_ia])) { - mt.insert(st.itia2iat(it, ia)); + mt.insert(ucell.itia2iat(it, ia)); find = true; break; } } if (!find) { - mag_type_atoms.push_back(std::set({ st.itia2iat(it,ia) })); + mag_type_atoms.push_back(std::set({ ucell.itia2iat(it,ia) })); } } } @@ -58,7 +59,7 @@ void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, i { // this->newpos have been ordered by original structure(ntype, na), it cannot be directly used here. // we need to reset the calculate again the coordinate of the new structure. - const ModuleBase::Vector3 direct_tmp = atoms[st.iat2it[mag_iat]].tau[st.iat2ia[mag_iat]] * this->optlat.Inverse(); + const ModuleBase::Vector3 direct_tmp = atoms[ucell.iat2it[mag_iat]].tau[ucell.iat2ia[mag_iat]] * this->optlat.Inverse(); std::array direct = { direct_tmp.x, direct_tmp.y, direct_tmp.z }; for (int i = 0; i < 3; ++i) { @@ -78,7 +79,7 @@ void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, i } -void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const Statistics& st, const ModuleBase::Matrix3& latvec) +void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const UnitCell& ucell, const ModuleBase::Matrix3& latvec) { // Restrict the space group to the unitary magnetic subgroup (nspin=4 / SOC): // operation g survives if it preserves the magnetic configuration as a pseudovector, @@ -97,7 +98,7 @@ void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const Statistics bool has_moment = false; for (int iat = 0; iat < this->nat && !has_moment; ++iat) { - const ModuleBase::Vector3& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]]; + const ModuleBase::Vector3& m = atoms[ucell.iat2it[iat]].m_loc_[ucell.iat2ia[iat]]; if (!this->equal(m.x, 0.0) || !this->equal(m.y, 0.0) || !this->equal(m.z, 0.0)) { has_moment = true; } } std::vector anti; // operations that REVERSE the moment: Theta*g is a symmetry @@ -109,11 +110,11 @@ void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const Statistics bool ok = true; for (int iat = 0; iat < this->nat && ok; ++iat) { - const ModuleBase::Vector3& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]]; + const ModuleBase::Vector3& m = atoms[ucell.iat2it[iat]].m_loc_[ucell.iat2ia[iat]]; // pseudovector-rotated moment W*m (column-vector convention: m'^i = W_ij m^j) const ModuleBase::Vector3& mrot = W * m; const int jat = this->get_rotated_atom(isym, iat); - const ModuleBase::Vector3& mj = atoms[st.iat2it[jat]].m_loc_[st.iat2ia[jat]]; + const ModuleBase::Vector3& mj = atoms[ucell.iat2it[jat]].m_loc_[ucell.iat2ia[jat]]; if (!this->equal(mrot.x, mj.x) || !this->equal(mrot.y, mj.y) || !this->equal(mrot.z, mj.z)) { ok = false; } } if (ok) @@ -129,10 +130,10 @@ void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const Statistics bool anti_ok = true; for (int iat = 0; iat < this->nat && anti_ok; ++iat) { - const ModuleBase::Vector3& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]]; + const ModuleBase::Vector3& m = atoms[ucell.iat2it[iat]].m_loc_[ucell.iat2ia[iat]]; const ModuleBase::Vector3& mrot = W * m; const int jat = this->get_rotated_atom(isym, iat); - const ModuleBase::Vector3& mj = atoms[st.iat2it[jat]].m_loc_[st.iat2ia[jat]]; + const ModuleBase::Vector3& mj = atoms[ucell.iat2it[jat]].m_loc_[ucell.iat2ia[jat]]; if (!this->equal(mrot.x, -mj.x) || !this->equal(mrot.y, -mj.y) || !this->equal(mrot.z, -mj.z)) { anti_ok = false; } } if (anti_ok) { anti.push_back(isym); } diff --git a/source/source_cell/module_symmetry/symm_rotation_k.cpp b/source/source_cell/module_symmetry/symm_rotation_k.cpp index f645419b05f..b56bbf4c1b8 100644 --- a/source/source_cell/module_symmetry/symm_rotation_k.cpp +++ b/source/source_cell/module_symmetry/symm_rotation_k.cpp @@ -1,4 +1,5 @@ #include "symm_rotation_k.h" +#include "source_cell/unitcell.h" #include "source_base/constants.h" #include #include "source_base/parallel_reduce.h" @@ -112,7 +113,7 @@ namespace ModuleSymmetry for (const int op : needed) { this->Ms_[ik_ibz][op] = this->contruct_2d_rot_mat_ao( - ucell.symm, ucell.atoms, ucell.st, kvec_d_ibz_global[ik_ibz], op, pv, spin_U[op]); + ucell.symm, ucell.atoms, ucell, kvec_d_ibz_global[ik_ibz], op, pv, spin_U[op]); } } @@ -388,7 +389,7 @@ namespace ModuleSymmetry // 2d-block parallized rotation matrix in AO-representation, denoted as M. // finally we will use D(k)=M(R, k)^\dagger*D(Rk)*M(R, k) to D(k) from D(Rk) in cal_Ms. - std::vector> Symmetry_rotation_k::contruct_2d_rot_mat_ao(const Symmetry& symm, const Atom* atoms, const Statistics& cell_st, + std::vector> Symmetry_rotation_k::contruct_2d_rot_mat_ao(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const TCdouble& kvec_d_ibz, int isym, const Parallel_2D& pv, const SpinRotation::Su2& spin_U) const { const bool soc = (this->nspin_ == 4); @@ -402,12 +403,12 @@ namespace ModuleSymmetry return (is < nrotk_u) ? symm.get_rotated_atom(is, iat) : symm.get_rotated_atom_anti(is - nrotk_u, iat); }; - for (int iat1 = 0;iat1 < cell_st.nat;++iat1) + for (int iat1 = 0;iat1 < ucell.nat;++iat1) { - int it = cell_st.iat2it[iat1]; // it1=it2 - int ia1 = cell_st.iat2ia[iat1]; + int it = ucell.iat2it[iat1]; // it1=it2 + int ia1 = ucell.iat2ia[iat1]; int iat2 = rotated_atom(isym, iat1); //iat2=rot(iat1) - int ia2 = cell_st.iat2ia[iat2]; + int ia2 = ucell.iat2ia[iat2]; // cal phase factor from return lattice: exp(-ik_ibz*O) double arg = -2 * ModuleBase::PI * kvec_d_ibz * this->irs_.return_lattice_[iat1][isym]; std::complexphase_factor = std::complex(std::cos(arg), std::sin(arg)); diff --git a/source/source_cell/module_symmetry/symm_rotation_k.h b/source/source_cell/module_symmetry/symm_rotation_k.h index 48ee839f4ff..d261dbed044 100644 --- a/source/source_cell/module_symmetry/symm_rotation_k.h +++ b/source/source_cell/module_symmetry/symm_rotation_k.h @@ -42,10 +42,10 @@ namespace ModuleSymmetry const int& abfs_Lmax = this->abfs_Lmax_; //-------------------------------------------------------------------------------- // setters - void find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const Statistics& st, + void find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const std::vector& Rs, const TC& period, const Lattice& lat, const std::string& output_dir = "") { - this->irs_.find_irreducible_sector(symm, atoms, st, Rs, period, lat, output_dir); + this->irs_.find_irreducible_sector(symm, atoms, ucell, Rs, period, lat, output_dir); } void set_abfs_Lmax(const int l) { this->abfs_Lmax_ = l; } //-------------------------------------------------------------------------------- @@ -105,7 +105,7 @@ namespace ModuleSymmetry /// 2d-block parallized rotation matrix in AO-representation, denoted as M. /// finally we will use D(k)=M(R, k)^\dagger*D(Rk)*M(R, k) to recover D(k) from D(Rk). - std::vector> contruct_2d_rot_mat_ao(const Symmetry& symm, const Atom* atoms, const Statistics& cell_st, + std::vector> contruct_2d_rot_mat_ao(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const TCdouble& kvec_d_ibz, int isym, const Parallel_2D& pv, const SpinRotation::Su2& spin_U /*= SpinRotation::Su2{ 1.0, 0.0, 0.0, 1.0 }*/) const; diff --git a/source/source_cell/module_symmetry/symmetry.cpp b/source/source_cell/module_symmetry/symmetry.cpp index c22774aeb02..b8e99fd962d 100644 --- a/source/source_cell/module_symmetry/symmetry.cpp +++ b/source/source_cell/module_symmetry/symmetry.cpp @@ -1,6 +1,7 @@ #include #include #include "symmetry.h" +#include "source_cell/unitcell.h" namespace ModuleSymmetry { @@ -365,13 +366,13 @@ void Symmetry::get_optlat(ModuleBase::Vector3 &v1, ModuleBase::Vector3 &a1, @@ -367,7 +369,7 @@ class Symmetry : public Symmetry_Basic /// Analyze magnetic group without time-reversal symmetry /// (because currently the charge density symmetrization does not support it) /// Method: treat atoms with different magmom as atoms of different type - void analyze_magnetic_group(const Atom* atoms, const Statistics& st, int& nrot_out, int& nrotk_out); + void analyze_magnetic_group(const Atom* atoms, const UnitCell& ucell, int& nrot_out, int& nrotk_out); /// (nspin=4 / SOC) Restrict the already-built space group to the unitary magnetic /// subgroup: keep operation g only if it preserves the magnetization as a pseudovector, @@ -375,7 +377,7 @@ class Symmetry : public Symmetry_Basic /// that reverse the moment (which are only symmetries when combined with time reversal) /// from being applied in k-reduction and density symmetrization. /// Non-magnetic (m_i=0) keeps all operations. - void analyze_magnetic_group_nspin4(const Atom* atoms, const Statistics& st, const ModuleBase::Matrix3& latvec); + void analyze_magnetic_group_nspin4(const Atom* atoms, const UnitCell& ucell, const ModuleBase::Matrix3& latvec); }; /** diff --git a/source/source_cell/module_symmetry/test/symm_test_analysis.cpp b/source/source_cell/module_symmetry/test/symm_test_analysis.cpp index 51d6c452601..1a2a47011ad 100644 --- a/source/source_cell/module_symmetry/test/symm_test_analysis.cpp +++ b/source/source_cell/module_symmetry/test/symm_test_analysis.cpp @@ -46,7 +46,7 @@ TEST_F(SymmetryTest, AnalySys) ModuleSymmetry::Symmetry symm; construct_ucell(stru_lib[stru]); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); //1. ibrav std::string ref_point_group = stru_lib[stru].point_group; @@ -250,7 +250,7 @@ TEST_F(SymmetryTest, SG_Pricell) symm.epsilon = 1e-5; construct_ucell(supercell_lib[stru]); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); std::string ref_point_group = supercell_lib[stru].point_group; std::string cal_point_group = symm.pgname; diff --git a/source/source_cell/module_symmetry/test/symm_test_symtrz.cpp b/source/source_cell/module_symmetry/test/symm_test_symtrz.cpp index 79f32a9a9cd..9911fa94755 100644 --- a/source/source_cell/module_symmetry/test/symm_test_symtrz.cpp +++ b/source/source_cell/module_symmetry/test/symm_test_symtrz.cpp @@ -71,7 +71,7 @@ TEST_F(SymmetryTest, ForceSymmetry) ModuleSymmetry::Symmetry symm; construct_ucell(supercell_lib[stru]); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); ModuleBase::matrix force(ucell.nat, 3, true); //generate random number for force and restrict to [-100,100) @@ -103,7 +103,7 @@ TEST_F(SymmetryTest, StressSymmetry) ModuleSymmetry::Symmetry symm; construct_ucell(supercell_lib[stru]); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); ModuleBase::matrix stress(3, 3, true); //generate random number for stress and restrict to [-1e5,1e5) diff --git a/source/source_cell/test/klist_test.cpp b/source/source_cell/test/klist_test.cpp index a9a546b2196..ecedde8484d 100644 --- a/source/source_cell/test/klist_test.cpp +++ b/source/source_cell/test/klist_test.cpp @@ -844,7 +844,7 @@ TEST_F(KlistTest, IbzKpoint) construct_ucell(stru_lib[0]); ofs_running.open("tmp_klist_3"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); std::string k_file = "./support/KPT1"; kv->set_spin_mult(1); kv->read_kpoints_for_testing(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); @@ -871,7 +871,7 @@ TEST_F(KlistTest, IbzKpointIsMP) construct_ucell(stru_lib[0]); ofs_running.open("tmp_klist_4"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); std::string k_file = "./support/KPT1"; kv->set_spin_mult(1); kv->read_kpoints_for_testing(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); @@ -899,7 +899,7 @@ TEST_F(KlistTest, IbzKpointCustomWeights) construct_ucell(stru_lib[0]); ofs_running.open("tmp_klist_custom_weights"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); // Test 1: Non-MP k-points with uniform weights (KPT4) { diff --git a/source/source_cell/test/klist_test_para.cpp b/source/source_cell/test/klist_test_para.cpp index 0bd7552a853..a4fca253cd3 100644 --- a/source/source_cell/test/klist_test_para.cpp +++ b/source/source_cell/test/klist_test_para.cpp @@ -178,7 +178,7 @@ TEST_F(KlistParaTest, Set) GlobalV::ofs_running.open("tmp_klist_5"); } const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); // read KPT std::string k_file = "./support/KPT1"; // note: do NOT pre-set kv->spin_mult here; set() takes the physical @@ -302,7 +302,7 @@ TEST_F(KlistParaTest, SetAfterVC) GlobalV::ofs_running.open("tmp_klist_6"); } const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); // read KPT std::string k_file = "./support/KPT1"; // note: do NOT pre-set kv->spin_mult here; set() takes the physical diff --git a/source/source_cell/test/little_group_test.cpp b/source/source_cell/test/little_group_test.cpp index 153dc9f501d..a5c0454a7ed 100644 --- a/source/source_cell/test/little_group_test.cpp +++ b/source/source_cell/test/little_group_test.cpp @@ -120,7 +120,7 @@ class LittleGroupTest : public testing::Test } std::ofstream ofs_running("tmp_little_group"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); } void TearDown() override diff --git a/source/source_cell/test/qlist_test.cpp b/source/source_cell/test/qlist_test.cpp index 637898edf44..959fa254725 100644 --- a/source/source_cell/test/qlist_test.cpp +++ b/source/source_cell/test/qlist_test.cpp @@ -173,7 +173,7 @@ TEST_F(QListTest, GenerateMeshFullSymmetry) ofs_running.open("tmp_qlist_1"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {8, 8, 8}, true); @@ -211,7 +211,7 @@ TEST_F(QListTest, GenerateMeshSmallGrid) ofs_running.open("tmp_qlist_2"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {2, 2, 2}, true); @@ -235,7 +235,7 @@ TEST_F(QListTest, GammaOnlyGrid) ofs_running.open("tmp_qlist_3"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {1, 1, 1}, true); @@ -255,7 +255,7 @@ TEST_F(QListTest, IrrepPlaceholder) ofs_running.open("tmp_qlist_4"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {2, 2, 2}, true); @@ -282,7 +282,7 @@ TEST_F(QListTest, CartesianCoordinatesComputed) ofs_running.open("tmp_qlist_cart"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {2, 2, 2}, true); @@ -312,7 +312,7 @@ TEST_F(QListTest, UseIrrepsSwitch) ofs_running.open("tmp_qlist_irreps"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); // use_irreps = false: the q mesh is still reduced, but no irrep data qlist.generate_mesh(ucell, symm, {2, 2, 2}, false); @@ -331,7 +331,7 @@ TEST_F(QListTest, PrintQlists) ofs_running.open("tmp_qlist_print"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {1, 1, 1}, false); diff --git a/source/source_cell/test/unitcell_test.cpp b/source/source_cell/test/unitcell_test.cpp index 22ed1fd421c..9e8b1883cde 100644 --- a/source/source_cell/test/unitcell_test.cpp +++ b/source/source_cell/test/unitcell_test.cpp @@ -145,10 +145,10 @@ TEST_F(UcellTest, Constructor) EXPECT_EQ(ucell->nat, 0); EXPECT_EQ(ucell->namax, 0); EXPECT_EQ(ucell->nwmax, 0); - EXPECT_EQ(ucell->iat2it, nullptr); - EXPECT_EQ(ucell->iat2ia, nullptr); - EXPECT_EQ(ucell->iwt2iat, nullptr); - EXPECT_EQ(ucell->iwt2iw, nullptr); + EXPECT_TRUE(ucell->iat2it.empty()); + EXPECT_TRUE(ucell->iat2ia.empty()); + EXPECT_TRUE(ucell->iwt2iat.empty()); + EXPECT_TRUE(ucell->iwt2iw.empty()); EXPECT_DOUBLE_EQ(ucell->tpiba, 0.0); EXPECT_DOUBLE_EQ(ucell->tpiba2, 0.0); EXPECT_DOUBLE_EQ(ucell->omega, 0.0); diff --git a/source/source_cell/unitcell.cpp b/source/source_cell/unitcell.cpp index cc4db426946..2dfbd22fe76 100644 --- a/source/source_cell/unitcell.cpp +++ b/source/source_cell/unitcell.cpp @@ -42,10 +42,8 @@ UnitCell::~UnitCell() void UnitCell::set_iat2itia() { assert(nat > 0); - delete[] iat2it; - delete[] iat2ia; - this->iat2it = new int[nat]; - this->iat2ia = new int[nat]; + this->iat2it.resize(nat); + this->iat2ia.resize(nat); int iat = 0; for (int it = 0; it < ntype; it++) { for (int ia = 0; ia < atoms[it].na; ia++) { diff --git a/source/source_cell/unitcell.h b/source/source_cell/unitcell.h index 359456737a3..32588c86e3c 100644 --- a/source/source_cell/unitcell.h +++ b/source/source_cell/unitcell.h @@ -2,6 +2,7 @@ #define UNITCELL_H #include +#include #include "source_base/global_function.h" #include "source_cell/sep_cell.h" #include "source_cell/magnetism.h" @@ -17,6 +18,14 @@ class UnitCell : public BaseCell { UnitCell(); ~UnitCell(); + /// UnitCell owns dynamically-sized state (index maps, atom arrays, MPI-broadcast + /// buffers). Copying is disabled to make accidental shallow-copy bugs fail at + /// compile time; moving is allowed. + UnitCell(const UnitCell&) = delete; + UnitCell& operator=(const UnitCell&) = delete; + UnitCell(UnitCell&&) = default; + UnitCell& operator=(UnitCell&&) = default; + /// @brief Initialize basic cell parameters (latname, ntype, lmaxmax, init_vel) /// from INPUT and parse fixed_axes into lat_axis_free flags. void setup_from_input(const std::string& latname_in, @@ -138,16 +147,15 @@ class UnitCell : public BaseCell { /// @name Statistics /// @{ - 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; + int ntype = 0; ///< number of atom species in UnitCell + int nat = 0; ///< total number of atoms of all species in unitcell + std::vector iat2it; ///< iat==>it, distinguish a atom belong to which type + std::vector iat2ia; ///< iat==>ia + std::vector iwt2iat; ///< iwt ==> iat + std::vector 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 = 0; ///< the max na among all atom species + int nwmax = 0; ///< the max nw among all atom species /// @} Atom* atoms = nullptr; diff --git a/source/source_cell/unitcell_data.h b/source/source_cell/unitcell_data.h index 5e185773239..5fc98a74540 100644 --- a/source/source_cell/unitcell_data.h +++ b/source/source_cell/unitcell_data.h @@ -34,42 +34,4 @@ struct Lattice ModuleBase::Matrix3 invGGT = ModuleBase::Matrix3(); ///< inverse G }; -/** - * @brief Statistics data and index maps. - * - * Relationships between indices: - * - ntype, it: atom type index - * - nat, iat: total atom index - * - atoms[it].na, ia: atom index within type - * - atoms[it].nw, iw: orbital index within atom - * - * - if know it ==> atoms[it].na; atoms[it].nw - * - if know iat ==> it; ia - * - if know ia, must have known it ==> iat - * - if know iwt, must have known it, ia ==> iwt - */ -struct Statistics -{ - int ntype = 0; ///< number of atom species in UnitCell - int nat = 0; ///< total number of atoms of all species in unitcell - int* iat2it = nullptr; ///< iat==>it, distinguish a atom belong to which type - int* iat2ia = nullptr; ///< iat==>ia - int* iwt2iat = nullptr; ///< iwt ==> iat - int* iwt2iw = nullptr; ///< 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 = 0; ///< the max na among all atom species - int nwmax = 0; ///< the max nw among all atom species - - /** - * @brief Destructor. - */ - ~Statistics() - { - delete[] iat2it; - delete[] iat2ia; - delete[] iwt2iat; - delete[] iwt2iw; - } -}; - #endif diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 5a661ce4a35..07bf211b482 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -119,7 +119,7 @@ void ESolver_FP::before_all_runners(BaseCell& basecell, const Input_para& inp) if (ModuleSymmetry::Symmetry::symm_flag == 1) { const int cal_symm_repr[2] = {this->inp_->cal_symm_repr[0], this->inp_->cal_symm_repr[1]}; - ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, + ucell.symm.analy_sys(ucell.lat, ucell, ucell.atoms, GlobalV::ofs_running, this->inp_->symmetry_prec, inp.nspin, this->inp_->calculation, cal_symm_repr); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } @@ -212,7 +212,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) if (ModuleSymmetry::Symmetry::symm_flag == 1) { const int cal_symm_repr[2] = {this->inp_->cal_symm_repr[0], this->inp_->cal_symm_repr[1]}; - ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, + ucell.symm.analy_sys(ucell.lat, ucell, ucell.atoms, GlobalV::ofs_running, this->inp_->symmetry_prec, this->inp_->nspin, this->inp_->calculation, cal_symm_repr); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } diff --git a/source/source_esolver/esolver_gets.cpp b/source/source_esolver/esolver_gets.cpp index 5c27a14a18f..91c72cde48a 100644 --- a/source/source_esolver/esolver_gets.cpp +++ b/source/source_esolver/esolver_gets.cpp @@ -69,7 +69,7 @@ void ESolver_GetS::before_all_runners(BaseCell& basecell, const Input_para& inp) { const int cal_symm_repr[2] = {this->inp_->cal_symm_repr[0], this->inp_->cal_symm_repr[1]}; ucell.symm.analy_sys(ucell.lat, - ucell.st, + ucell, ucell.atoms, GlobalV::ofs_running, this->inp_->symmetry_prec, diff --git a/source/source_esolver/esolver_lr_lcao_tddft.cpp b/source/source_esolver/esolver_lr_lcao_tddft.cpp index e754ab02041..71ae7f1f4b7 100644 --- a/source/source_esolver/esolver_lr_lcao_tddft.cpp +++ b/source/source_esolver/esolver_lr_lcao_tddft.cpp @@ -374,7 +374,7 @@ void ModuleESolver::ESolver_LR::initialize_from_unitcell_(UnitCell& ucell if (ModuleSymmetry::Symmetry::symm_flag == 1) { const int cal_symm_repr[2] = {this->inp_->cal_symm_repr[0], this->inp_->cal_symm_repr[1]}; - ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, + ucell.symm.analy_sys(ucell.lat, ucell, ucell.atoms, GlobalV::ofs_running, this->inp_->symmetry_prec, this->inp_->nspin, this->inp_->calculation, cal_symm_repr); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } diff --git a/source/source_esolver/test/esolver_dp_test.cpp b/source/source_esolver/test/esolver_dp_test.cpp index a60b40a0783..4bd3c445015 100644 --- a/source/source_esolver/test/esolver_dp_test.cpp +++ b/source/source_esolver/test/esolver_dp_test.cpp @@ -26,10 +26,10 @@ class ESolverDPTest : public ::testing::Test { // Initialize variables before each test esolver = new ModuleESolver::ESolver_DP("./support/case_1.pb"); - ucell.iat2it = new int[2]; + ucell.iat2it.resize(2); ucell.iat2it[0] = 0; ucell.iat2it[1] = 1; - ucell.iat2ia = new int[2]; + ucell.iat2ia.resize(2); ucell.iat2ia[0] = 0; ucell.iat2ia[1] = 1; ucell.nat = 2; diff --git a/source/source_estate/module_dm/unittests/prepare_unitcell.h b/source/source_estate/module_dm/unittests/prepare_unitcell.h index 43e00e8f765..8faec588075 100644 --- a/source/source_estate/module_dm/unittests/prepare_unitcell.h +++ b/source/source_estate/module_dm/unittests/prepare_unitcell.h @@ -223,14 +223,10 @@ class UcellTestPrepare ucell.namax = namax; // cal_index assert(nlocal > 0); - delete[] ucell.iwt2iat; - delete[] ucell.iwt2iw; - delete[] ucell.iat2it; - delete[] ucell.iat2ia; - ucell.iwt2iat = new int[nlocal]; - ucell.iwt2iw = new int[nlocal]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; // set_iat2itia + ucell.iwt2iat.resize(nlocal); + ucell.iwt2iw.resize(nlocal); + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); // set_iat2itia ucell.itia2iat.create(ucell.ntype, ucell.namax); ucell.set_iat2iwt(1); int iat = 0; diff --git a/source/source_estate/module_dm/unittests/test_cal_dm_r.cpp b/source/source_estate/module_dm/unittests/test_cal_dm_r.cpp index dbad57eb1de..9176463cb3f 100644 --- a/source/source_estate/module_dm/unittests/test_cal_dm_r.cpp +++ b/source/source_estate/module_dm/unittests/test_cal_dm_r.cpp @@ -38,8 +38,8 @@ class DMTest : public testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_estate/module_dm/unittests/test_dm_constructor.cpp b/source/source_estate/module_dm/unittests/test_dm_constructor.cpp index f270439fe6f..a6b0062dd02 100644 --- a/source/source_estate/module_dm/unittests/test_dm_constructor.cpp +++ b/source/source_estate/module_dm/unittests/test_dm_constructor.cpp @@ -37,8 +37,8 @@ class DMTest : public testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_estate/module_dm/unittests/test_dm_from_psi.cpp b/source/source_estate/module_dm/unittests/test_dm_from_psi.cpp index 3988193ff72..49ac90cf4a2 100644 --- a/source/source_estate/module_dm/unittests/test_dm_from_psi.cpp +++ b/source/source_estate/module_dm/unittests/test_dm_from_psi.cpp @@ -38,8 +38,8 @@ class DMTest : public testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau = new ModuleBase::Vector3[ucell.nat]; ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_estate/module_dm/unittests/test_dm_r_init.cpp b/source/source_estate/module_dm/unittests/test_dm_r_init.cpp index 8c1a6b988e5..bce2a29a776 100644 --- a/source/source_estate/module_dm/unittests/test_dm_r_init.cpp +++ b/source/source_estate/module_dm/unittests/test_dm_r_init.cpp @@ -38,8 +38,8 @@ class DMTest : public testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_estate/test/test_occ_mixer.cpp b/source/source_estate/test/test_occ_mixer.cpp index 33f9a36e4e7..2d2f04e0a87 100644 --- a/source/source_estate/test/test_occ_mixer.cpp +++ b/source/source_estate/test/test_occ_mixer.cpp @@ -45,11 +45,9 @@ class OccMatMixerTest : public ::testing::Test cell.ntype = 1; cell.nat = 2; atoms_storage.resize(cell.ntype); - // UnitCell exposes raw pointers; Statistics::~Statistics() owns and - // delete[]s iat2it/iat2ia, so they must come from new[] here. cell.atoms = atoms_storage.data(); - cell.iat2it = new int[cell.nat]; - cell.iat2ia = new int[cell.nat]; + cell.iat2it.resize(cell.nat); + cell.iat2ia.resize(cell.nat); cell.itia2iat.create(cell.ntype, cell.nat); cell.atoms[0].na = 2; cell.atoms[0].nwl = 2; // max angular momentum present diff --git a/source/source_hamilt/module_gint/test/test_gint_common.cpp b/source/source_hamilt/module_gint/test/test_gint_common.cpp index 4d6a21607ca..478130d0e30 100644 --- a/source/source_hamilt/module_gint/test/test_gint_common.cpp +++ b/source/source_hamilt/module_gint/test/test_gint_common.cpp @@ -13,8 +13,8 @@ class GintCommonTest : public ::testing::Test ucell_.ntype = 1; ucell_.nat = 2; ucell_.atoms = new Atom[ucell_.ntype]; - ucell_.iat2it = new int[ucell_.nat]; - ucell_.iat2ia = new int[ucell_.nat]; + ucell_.iat2it.resize(ucell_.nat); + ucell_.iat2ia.resize(ucell_.nat); ucell_.atoms[0].na = 2; ucell_.atoms[0].nw = 2; diff --git a/source/source_hamilt/module_hcontainer/test/prepare_unitcell.h b/source/source_hamilt/module_hcontainer/test/prepare_unitcell.h index 834a177bef8..e8123a074dd 100644 --- a/source/source_hamilt/module_hcontainer/test/prepare_unitcell.h +++ b/source/source_hamilt/module_hcontainer/test/prepare_unitcell.h @@ -223,14 +223,10 @@ class UcellTestPrepare ucell.namax = namax; // cal_index assert(nlocal > 0); - delete[] ucell.iwt2iat; - delete[] ucell.iwt2iw; - delete[] ucell.iat2it; - delete[] ucell.iat2ia; - ucell.iwt2iat = new int[nlocal]; - ucell.iwt2iw = new int[nlocal]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; // set_iat2itia + ucell.iwt2iat.resize(nlocal); + ucell.iwt2iw.resize(nlocal); + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); // set_iat2itia ucell.itia2iat.create(ucell.ntype, ucell.namax); ucell.set_iat2iwt(1); int iat = 0; diff --git a/source/source_hamilt/module_hcontainer/test/test_func_folding.cpp b/source/source_hamilt/module_hcontainer/test/test_func_folding.cpp index 12f5c9b987e..3ee04d5e360 100644 --- a/source/source_hamilt/module_hcontainer/test/test_func_folding.cpp +++ b/source/source_hamilt/module_hcontainer/test/test_func_folding.cpp @@ -25,7 +25,7 @@ class FoldingTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) { ucell.iat2it[iat] = 0; diff --git a/source/source_hamilt/module_hcontainer/test/test_hcontainer.cpp b/source/source_hamilt/module_hcontainer/test/test_hcontainer.cpp index 991b074a4a0..dc2392aff66 100644 --- a/source/source_hamilt/module_hcontainer/test/test_hcontainer.cpp +++ b/source/source_hamilt/module_hcontainer/test/test_hcontainer.cpp @@ -25,7 +25,7 @@ class HContainerTest : public ::testing::Test ucell.ntype = 1; ucell.nat = 3; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) { ucell.iat2it[iat] = 0; diff --git a/source/source_hamilt/module_hcontainer/test/test_hcontainer_complex.cpp b/source/source_hamilt/module_hcontainer/test/test_hcontainer_complex.cpp index c7920e202dd..d0cf0ed52f0 100644 --- a/source/source_hamilt/module_hcontainer/test/test_hcontainer_complex.cpp +++ b/source/source_hamilt/module_hcontainer/test/test_hcontainer_complex.cpp @@ -25,7 +25,7 @@ class HContainerTest : public ::testing::Test ucell.ntype = 1; ucell.nat = 3; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) { ucell.iat2it[iat] = 0; diff --git a/source/source_hamilt/module_hcontainer/test/test_hcontainer_output.cpp b/source/source_hamilt/module_hcontainer/test/test_hcontainer_output.cpp index 443841ffc54..d37ac7e32c1 100644 --- a/source/source_hamilt/module_hcontainer/test/test_hcontainer_output.cpp +++ b/source/source_hamilt/module_hcontainer/test/test_hcontainer_output.cpp @@ -32,8 +32,8 @@ class OutputHContainerTest : public testing::Test ucell.ntype = 1; ucell.nat = 2; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) { ucell.iat2ia[iat] = iat; @@ -41,8 +41,8 @@ class OutputHContainerTest : public testing::Test } ucell.atoms[0].na = 2; ucell.atoms[0].nw = 2; - ucell.iwt2iat = new int[4]; - ucell.iwt2iw = new int[4]; + ucell.iwt2iat.resize(4); + ucell.iwt2iw.resize(4); ucell.itia2iat.create(ucell.ntype, ucell.nat); ucell.set_iat2iwt(1); ucell.itia2iat(0, 0) = 0; diff --git a/source/source_hamilt/module_hcontainer/test/test_hcontainer_time.cpp b/source/source_hamilt/module_hcontainer/test/test_hcontainer_time.cpp index 0ac84c718ba..6afe0976cb7 100644 --- a/source/source_hamilt/module_hcontainer/test/test_hcontainer_time.cpp +++ b/source/source_hamilt/module_hcontainer/test/test_hcontainer_time.cpp @@ -24,7 +24,7 @@ class HContainerTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) { ucell.iat2it[iat] = 0; diff --git a/source/source_hamilt/module_hcontainer/test/test_transfer.cpp b/source/source_hamilt/module_hcontainer/test/test_transfer.cpp index 385df86b6c1..26500c6d5fd 100644 --- a/source/source_hamilt/module_hcontainer/test/test_transfer.cpp +++ b/source/source_hamilt/module_hcontainer/test/test_transfer.cpp @@ -27,8 +27,8 @@ class TransferTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) { diff --git a/source/source_io/module_dos/unittests/test_write_pdos_text.cpp b/source/source_io/module_dos/unittests/test_write_pdos_text.cpp index 2ac56f31727..2ab6f9009b7 100644 --- a/source/source_io/module_dos/unittests/test_write_pdos_text.cpp +++ b/source/source_io/module_dos/unittests/test_write_pdos_text.cpp @@ -32,8 +32,8 @@ static UnitCell* make_si_ucell() atom.l_nchi = {1, 1}; // iat2it / iat2ia - ucell->iat2it = new int[1]; - ucell->iat2ia = new int[1]; + ucell->iat2it.resize(1); + ucell->iat2ia.resize(1); ucell->iat2it[0] = 0; ucell->iat2ia[0] = 0; @@ -57,8 +57,6 @@ class WritePdosTextTest : public ::testing::Test } void TearDown() override { - // iat2it/iat2ia are owned by UnitCell's internal Statistics member, - // whose destructor releases them; do not delete them here. delete ucell; } diff --git a/source/source_io/module_hs/unittests/test_hsr_writer.cpp b/source/source_io/module_hs/unittests/test_hsr_writer.cpp index f4a9e6c90d8..659456da411 100644 --- a/source/source_io/module_hs/unittests/test_hsr_writer.cpp +++ b/source/source_io/module_hs/unittests/test_hsr_writer.cpp @@ -379,7 +379,7 @@ TEST(HsrWriterIo, HContainerBinaryMpiGatherWritesCompleteFiles) ucell.set_atom_flag = true; ucell.atoms[0].na = 1; ucell.atoms[0].nw = 2; - ucell.iat2it = new int[1]; + ucell.iat2it.resize(1); ucell.iat2it[0] = 0; ucell.set_iat2iwt(1); const int* iat2iwt = ucell.get_iat2iwt(); diff --git a/source/source_io/test/restart_exx_csr_test.cpp b/source/source_io/test/restart_exx_csr_test.cpp index eee5fd96e50..4da9b0cef4e 100644 --- a/source/source_io/test/restart_exx_csr_test.cpp +++ b/source/source_io/test/restart_exx_csr_test.cpp @@ -15,8 +15,8 @@ void init_unitcell_for_ri(UnitCell& ucell) ucell.atoms[0].na = 1; ucell.atoms[0].nw = 2; ucell.atoms[0].stapos_wf = 0; - ucell.iat2it = new int[1]{0}; - ucell.iat2ia = new int[1]{0}; + ucell.iat2it = {0}; + ucell.iat2ia = {0}; } } // namespace diff --git a/source/source_lcao/module_dftu/dftu_nao_occ.cpp b/source/source_lcao/module_dftu/dftu_nao_occ.cpp index e535ed59ec2..13a57778f67 100644 --- a/source/source_lcao/module_dftu/dftu_nao_occ.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_occ.cpp @@ -362,7 +362,7 @@ void cal_occ_mat_k(const Parallel_Orbitals* pv, if (dftu_spacegroup_symmetry && !dftu_occ_symrot_built) { const std::array period{ kv.nmp[0], kv.nmp[1], kv.nmp[2] }; - dftu_occ_symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st, + dftu_occ_symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell, ModuleSymmetry::Symmetry_rotation_k::get_bvk_cells(period), period, ucell.lat); dftu_occ_symrot.cal_Ms(kv, ucell, *pv, nspin); dftu_occ_symrot_built = true; diff --git a/source/source_lcao/module_dftu/dftu_nao_op.cpp b/source/source_lcao/module_dftu/dftu_nao_op.cpp index 022672f386d..6a61dc1cb2d 100644 --- a/source/source_lcao/module_dftu/dftu_nao_op.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_op.cpp @@ -141,7 +141,7 @@ void hamilt::DFTU_onsite>::contributeHR() { const std::array period{ this->kv_->nmp[0], this->kv_->nmp[1], this->kv_->nmp[2] }; // for return_lattice to calculate Ms - this->symrot_.find_irreducible_sector(this->ucell->symm, this->ucell->atoms, this->ucell->st, + this->symrot_.find_irreducible_sector(this->ucell->symm, this->ucell->atoms, *this->ucell, ModuleSymmetry::Symmetry_rotation_k::get_bvk_cells(period), period, this->ucell->lat); this->symrot_.cal_Ms(*this->kv_, *this->ucell, *pv, this->nspin); this->symrot_built_ = true; diff --git a/source/source_lcao/module_dftu/unittests/test_dftu_nao_fs_accum.cpp b/source/source_lcao/module_dftu/unittests/test_dftu_nao_fs_accum.cpp index f5f01e828a0..ac1a7ba014f 100644 --- a/source/source_lcao/module_dftu/unittests/test_dftu_nao_fs_accum.cpp +++ b/source/source_lcao/module_dftu/unittests/test_dftu_nao_fs_accum.cpp @@ -33,8 +33,7 @@ class DftuFsAccumTest : public ::testing::Test pv.set_serial(3, 3); // iwt2iat: orbital 0,1 -> atom 0; orbital 2 -> atom 1 - // Must use new[] because Statistics destructor will delete[] it. - ucell.iwt2iat = new int[3]{0, 0, 1}; + ucell.iwt2iat = {0, 0, 1}; } Parallel_Orbitals pv; diff --git a/source/source_lcao/module_dftu/unittests/test_dftu_nao_ijr.cpp b/source/source_lcao/module_dftu/unittests/test_dftu_nao_ijr.cpp index d70848daf18..0083c0d7dfb 100644 --- a/source/source_lcao/module_dftu/unittests/test_dftu_nao_ijr.cpp +++ b/source/source_lcao/module_dftu/unittests/test_dftu_nao_ijr.cpp @@ -69,16 +69,14 @@ class AccumulateHrIat0Test : public ::testing::Test ucell.ntype = 1; ucell.nat = nat; ucell.atoms = atoms_buf; - iat2it_buf.assign(nat, 0); - iat2ia_buf.resize(nat); - ucell.iat2it = iat2it_buf.data(); - ucell.iat2ia = iat2ia_buf.data(); + ucell.iat2it.assign(nat, 0); + ucell.iat2ia.resize(nat); ucell.atoms[0].tau.resize(nat); ucell.lat0 = 1.0; ucell.itia2iat.create(ucell.ntype, nat); for (int iat = 0; iat < nat; iat++) { - iat2ia_buf[iat] = iat; + ucell.iat2ia[iat] = iat; ucell.atoms[0].tau[iat] = ModuleBase::Vector3(0.0, 0.0, 0.0); ucell.itia2iat(0, iat) = iat; } @@ -102,8 +100,6 @@ class AccumulateHrIat0Test : public ::testing::Test HR.reset(); paraV.reset(); ucell.atoms = nullptr; - ucell.iat2it = nullptr; - ucell.iat2ia = nullptr; } AdjacentAtomInfo make_adjs() const @@ -142,8 +138,6 @@ class AccumulateHrIat0Test : public ::testing::Test int my_rank = 0; UnitCell ucell; Atom atoms_buf[1]; - std::vector iat2it_buf; - std::vector iat2ia_buf; std::unique_ptr paraV; std::unique_ptr> HR; }; diff --git a/source/source_lcao/module_dftu/unittests/test_dftu_nao_op.cpp b/source/source_lcao/module_dftu/unittests/test_dftu_nao_op.cpp index 0a87a8a1578..132131147c1 100644 --- a/source/source_lcao/module_dftu/unittests/test_dftu_nao_op.cpp +++ b/source/source_lcao/module_dftu/unittests/test_dftu_nao_op.cpp @@ -31,8 +31,8 @@ class DFTUTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.lat0 = 1.0; ucell.itia2iat.create(ucell.ntype, ucell.nat); diff --git a/source/source_lcao/module_lr/ri_benchmark/test/ri_benchmark_test.cpp b/source/source_lcao/module_lr/ri_benchmark/test/ri_benchmark_test.cpp index 2a23cc7b420..407697dd869 100644 --- a/source/source_lcao/module_lr/ri_benchmark/test/ri_benchmark_test.cpp +++ b/source/source_lcao/module_lr/ri_benchmark/test/ri_benchmark_test.cpp @@ -20,7 +20,7 @@ Sep_Cell::Sep_Cell() noexcept {} Sep_Cell::~Sep_Cell() noexcept {} UnitCell::UnitCell() { atoms = new Atom[1]; - iat2it = new int[1]; iat2it[0] = 0; + iat2it = {0}; iat2iwt.resize(1, 0); } UnitCell::~UnitCell() { diff --git a/source/source_lcao/module_operator_lcao/test/test_ekinetic.cpp b/source/source_lcao/module_operator_lcao/test/test_ekinetic.cpp index 9edc265ccbf..d241ca82e03 100644 --- a/source/source_lcao/module_operator_lcao/test/test_ekinetic.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_ekinetic.cpp @@ -33,8 +33,8 @@ class EKineticTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_lcao/module_operator_lcao/test/test_ekinetic_serial.cpp b/source/source_lcao/module_operator_lcao/test/test_ekinetic_serial.cpp index 87ce65d4342..9d035bd4bd0 100644 --- a/source/source_lcao/module_operator_lcao/test/test_ekinetic_serial.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_ekinetic_serial.cpp @@ -33,8 +33,8 @@ class EKineticTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) @@ -191,8 +191,8 @@ TEST_F(EKineticTest, singleAtom) ucell_single.ntype = 1; ucell_single.nat = 1; ucell_single.atoms = new Atom[1]; - ucell_single.iat2it = new int[1]; - ucell_single.iat2ia = new int[1]; + ucell_single.iat2it.resize(1); + ucell_single.iat2ia.resize(1); ucell_single.atoms[0].tau.resize(1); ucell_single.itia2iat.create(1, 1); ucell_single.iat2it[0] = 0; diff --git a/source/source_lcao/module_operator_lcao/test/test_nonlocal.cpp b/source/source_lcao/module_operator_lcao/test/test_nonlocal.cpp index 5eea2ba121a..d604dab5f77 100644 --- a/source/source_lcao/module_operator_lcao/test/test_nonlocal.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_nonlocal.cpp @@ -38,8 +38,8 @@ class NonlocalTest : public ::testing::Test ucell.infoNL.reset(lcao_nl); ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_lcao/module_operator_lcao/test/test_overlap.cpp b/source/source_lcao/module_operator_lcao/test/test_overlap.cpp index e3c017857b2..aa6a14ae5a0 100644 --- a/source/source_lcao/module_operator_lcao/test/test_overlap.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_overlap.cpp @@ -33,8 +33,8 @@ class OverlapTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_lcao/module_operator_lcao/test/test_overlap_cd.cpp b/source/source_lcao/module_operator_lcao/test/test_overlap_cd.cpp index ad2e8daf7e1..eca29b00dda 100644 --- a/source/source_lcao/module_operator_lcao/test/test_overlap_cd.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_overlap_cd.cpp @@ -32,8 +32,8 @@ class OverlapTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_lcao/module_operator_lcao/test/test_overlap_serial.cpp b/source/source_lcao/module_operator_lcao/test/test_overlap_serial.cpp index 98b8e4d7419..9a69b82bb7f 100644 --- a/source/source_lcao/module_operator_lcao/test/test_overlap_serial.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_overlap_serial.cpp @@ -33,8 +33,8 @@ class OverlapTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) @@ -198,8 +198,8 @@ TEST_F(OverlapTest, singleAtom) ucell_single.ntype = 1; ucell_single.nat = 1; ucell_single.atoms = new Atom[1]; - ucell_single.iat2it = new int[1]; - ucell_single.iat2ia = new int[1]; + ucell_single.iat2it.resize(1); + ucell_single.iat2ia.resize(1); ucell_single.atoms[0].tau.resize(1); ucell_single.itia2iat.create(1, 1); ucell_single.iat2it[0] = 0; diff --git a/source/source_lcao/module_operator_lcao/test/test_t_nl_cd.cpp b/source/source_lcao/module_operator_lcao/test/test_t_nl_cd.cpp index 5b9e6c1c572..a83235b2816 100644 --- a/source/source_lcao/module_operator_lcao/test/test_t_nl_cd.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_t_nl_cd.cpp @@ -35,8 +35,8 @@ class TNLTest : public ::testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); for (int iat = 0; iat < ucell.nat; iat++) diff --git a/source/source_lcao/module_rdmft/rdmft.cpp b/source/source_lcao/module_rdmft/rdmft.cpp index 22883e46d6d..16c94e978dd 100644 --- a/source/source_lcao/module_rdmft/rdmft.cpp +++ b/source/source_lcao/module_rdmft/rdmft.cpp @@ -157,7 +157,7 @@ void RDMFT::init(Parallel_Orbitals& ParaV_in, if (exx_spacegroup_symmetry) { const std::array& period = RI_Util::get_Born_vonKarmen_period(*kv); - this->symrot_exx.find_irreducible_sector(ucell->symm, ucell->atoms, ucell->st, + this->symrot_exx.find_irreducible_sector(ucell->symm, ucell->atoms, *ucell, RI_Util::get_Born_von_Karmen_cells(period), period, ucell->lat, PARAM.globalv.global_out_dir); this->symrot_exx.cal_Ms(*kv, *ucell, *ParaV, nspin); } diff --git a/source/source_lcao/module_ri/exx_lri.hpp b/source/source_lcao/module_ri/exx_lri.hpp index 04a134a7748..e202851ca9e 100644 --- a/source/source_lcao/module_ri/exx_lri.hpp +++ b/source/source_lcao/module_ri/exx_lri.hpp @@ -822,7 +822,7 @@ void Exx_LRI::cal_exx_elec(const std::vectorexx_lri.post_2D.set_tensors_map2(this->exx_lri.Hs); // rotate locally without repeat - Hs_a2D = p_symrot->restore_HR(ucell.symm, ucell.atoms, ucell.st, 'H', Hs_a2D); + Hs_a2D = p_symrot->restore_HR(ucell.symm, ucell.atoms, ucell, 'H', Hs_a2D); // cal energy using full Hs without repeat this->exx_lri.energy = this->exx_lri.post_2D.cal_energy( this->exx_lri.post_2D.saves["Ds_" + suffix], @@ -864,7 +864,7 @@ void Exx_LRI::cal_exx_elec_soc( // pass 2: spinor-coupled rotation of the 4 channels from the irreducible sector to the full BZ std::array>>, 4> Hs_full = - p_symrot->restore_HR_nspin4(ucell.symm, ucell.atoms, ucell.st, 'H', Hs_irr); + p_symrot->restore_HR_nspin4(ucell.symm, ucell.atoms, ucell, 'H', Hs_irr); // pass 3: per-channel energy (full Hs, no repeat), then gather the repeated full Hs for abacus for (int is = 0; is < 4; ++is) diff --git a/source/source_lcao/module_ri/exx_lri_interface.hpp b/source/source_lcao/module_ri/exx_lri_interface.hpp index c09592f5fca..20b3c39aee1 100644 --- a/source/source_lcao/module_ri/exx_lri_interface.hpp +++ b/source/source_lcao/module_ri/exx_lri_interface.hpp @@ -91,7 +91,7 @@ void Exx_LRI_Interface::exx_before_all_runners( { const std::array& period = RI_Util::get_Born_vonKarmen_period(kv); this->symrot_.find_irreducible_sector( - ucell.symm, ucell.atoms, ucell.st, + ucell.symm, ucell.atoms, ucell, RI_Util::get_Born_von_Karmen_cells(period), period, ucell.lat, PARAM.globalv.global_out_dir); this->symrot_.set_abfs_Lmax(Exx_Abfs::Construct_Orbs::get_Lmax(this->exx_ptr->abfs)); this->symrot_.cal_Ms(kv, ucell, pv, PARAM.inp.nspin); diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h index 16d7ba91928..4dfded3ddbf 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h @@ -6,6 +6,7 @@ #include #include "source_hamilt/module_hcontainer/hcontainer.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" +#include "source_cell/unitcell.h" namespace ModuleSymmetry { @@ -31,11 +32,11 @@ namespace ModuleSymmetry /// Given H(R) in the irreduceble sector, calculate H(R) for all the atompairs and cells. template // RI::Tensor type std::map, RI::Tensor>> restore_HR( - const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, const std::map, RI::Tensor>>& HR_irreduceble)const; template // HContainer type void restore_HR( - const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, const hamilt::HContainer& HR_irreduceble, hamilt::HContainer& HR_rotated)const; /// (nspin=4) spinor overload: rotate all 4 spin channels of H(R) together. On top of the /// orbital rotation T1^dagger(.)T2 (mode 'H') / T1^T(.)T2^* (mode 'D') applied to every @@ -47,20 +48,20 @@ namespace ModuleSymmetry /// remap of the 4 channels applied after the SU(2) mixing; see symm_rotation_r.hpp. template // RI::Tensor type std::array, RI::Tensor>>, 4> restore_HR_nspin4( - const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, const std::array, RI::Tensor>>, 4>& HR_irreducible_soc)const; //-------------------------------------------------------------------------------- /// test functions /// test H(R) rotation: giver a full H(R), pick out H(R) in the irreducible sector, rotate it, and compare with the original full H(R) template // RI::Tensor type, using col-major implementation - void test_HR_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + void test_HR_rotation(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, const std::map, RI::Tensor>>& HR_full); template // test the rotation of RI coefficients - void test_Cs_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, + void test_Cs_rotation(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const std::map, RI::Tensor>>& Cs_full)const; template // HContainer type, using row-major implementation - void test_HR_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + void test_HR_rotation(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, const hamilt::HContainer& HR_full); template // HContainer type void print_HR(const std::map, RI::Tensor>>& HR, const std::string name, const double& threshold = 0.0); @@ -84,7 +85,7 @@ namespace ModuleSymmetry void rotate_atompair_serial(Tdata* TAT, const Tdata* A, const int& nw1, const int& nw2, const int isym, const Atom& a1, const Atom& a2, const char mode)const; template // HContainer type, pblas - void rotate_atompair_parallel(const TR* Alocal_in, const int isym, const Atom* atoms, const Statistics& st, + void rotate_atompair_parallel(const TR* Alocal_in, const int isym, const Atom* atoms, const UnitCell& ucell, const Tap& ap_in, const Tap& ap_out, const char mode, const Parallel_Orbitals& pv, TR* Alocal_out, const bool output = false)const; /// rotate a 3-dim C tensor in RI technique diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp index 381aaf5f435..c896f085bb8 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp @@ -77,7 +77,7 @@ namespace ModuleSymmetry template std::map, RI::Tensor>> Symmetry_rotation::restore_HR( - const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, const std::map, RI::Tensor>>& HR_irreducible) const { ModuleBase::TITLE("Symmetry_rotation", "restore_HR"); @@ -102,7 +102,7 @@ namespace ModuleSymmetry const int& ap1 = apR.first.first; const int& ap2 = apR.first.second; const TC& R = apR.second; - HR_full[ap1][{ap2, R}] = rotate_atompair_serial(tmp2.second, isym, atoms[st.iat2it[irap1]], atoms[st.iat2it[irap2]], mode); + HR_full[ap1][{ap2, R}] = rotate_atompair_serial(tmp2.second, isym, atoms[ucell.iat2it[irap1]], atoms[ucell.iat2it[irap2]], mode); } } else { std::cout << "Warning: not found: irreducible atom pair =(" << irap1 << "," << irap2 << "), irR=(" << irR[0] << "," << irR[1] << "," << irR[2] << ")\n";} @@ -121,7 +121,7 @@ namespace ModuleSymmetry // // H_12(R)=T^\dagger(V)H_1'2'(VR+O_1-O_2)T(V) // if (HR_irreducible.find(irap.first) != HR_irreducible.end() && HR_irreducible.at(irap.first).find({ irap.second, irR }) != HR_irreducible.at(irap.first).end()) // HR_full[ap.first][{ap.second, R}] = rotate_atompair_serial(HR_irreducible.at(irap.first).at({ irap.second, irR }), - // isym, atoms[st.iat2it[irap.first]], atoms[st.iat2it[irap.second]], mode); + // isym, atoms[ucell.iat2it[irap.first]], atoms[ucell.iat2it[irap.second]], mode); // else // std::cout << "not found: current atom pair =(" << ap.first << "," << ap.second << "), R=(" << R[0] << "," << R[1] << "," << R[2] << "), irreducible atom pair =(" << irap.first << "," << irap.second << "), irR=(" << irR[0] << "," << irR[1] << "," << irR[2] << ")\n"; // } @@ -171,7 +171,7 @@ namespace ModuleSymmetry // when it is applied in either direction. template std::array, RI::Tensor>>, 4> Symmetry_rotation::restore_HR_nspin4( - const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, const std::array, RI::Tensor>>, 4>& HR_irreducible_soc)const { ModuleBase::TITLE("Symmetry_rotation", "restore_HR_nspin4"); @@ -203,8 +203,8 @@ namespace ModuleSymmetry std::cout << "Warning: not found: irreducible atom pair =(" << irap1 << "," << irap2 << "), irR=(" << irR[0] << "," << irR[1] << "," << irR[2] << ")\n"; continue; } - const Atom& a1 = atoms[st.iat2it[irap1]]; - const Atom& a2 = atoms[st.iat2it[irap2]]; + const Atom& a1 = atoms[ucell.iat2it[irap1]]; + const Atom& a2 = atoms[ucell.iat2it[irap2]]; // gather the 4 irreducible spin-channel blocks (zero-filled when absent) std::array, 4> Hir; for (int is = 0;is < 4;++is) @@ -406,7 +406,7 @@ namespace ModuleSymmetry } template - void Symmetry_rotation::test_HR_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + void Symmetry_rotation::test_HR_rotation(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, const std::map, RI::Tensor>>& HR_full) { ModuleBase::TITLE("Symmetry_rotation", "test_HR_rotation"); @@ -425,7 +425,7 @@ namespace ModuleSymmetry } } // 2. rotate - std::map, RI::Tensor>> HR_rotated = restore_HR(symm, atoms, st, mode, HR_irreducible); + std::map, RI::Tensor>> HR_rotated = restore_HR(symm, atoms, ucell, mode, HR_irreducible); // 3. compare for (auto& HR_ia1 : HR_rotated) { @@ -452,7 +452,7 @@ namespace ModuleSymmetry } template - void Symmetry_rotation::test_Cs_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, + void Symmetry_rotation::test_Cs_rotation(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const std::map, RI::Tensor>>& Cs_full)const { for (auto& sector_pair : this->irs_.full_map_to_irreducible_sector_) @@ -468,7 +468,7 @@ namespace ModuleSymmetry const RI::Tensor& Cs_ir = Cs_full.at(irapR.first.first).at({ irapR.first.second,irapR.second }); const RI::Tensor& Cs_ref = Cs_full.at(apR.first.first).at({ apR.first.second,apR.second }); const RI::Tensor& Cs_rot = this->rotate_singleC_serial(Cs_ir, isym, - atoms[st.iat2it[irapR.first.first]], atoms[st.iat2it[irapR.first.second]], irapR.first.first); + atoms[ucell.iat2it[irapR.first.first]], atoms[ucell.iat2it[irapR.first.second]], irapR.first.first); print_tensor3(Cs_rot, "Cs_rot"); print_tensor3(Cs_ref, "Cs_ref"); print_tensor3(Cs_ir, "Cs_irreducible"); diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r_hcontainer.hpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r_hcontainer.hpp index 7e96480dabd..61d5e344490 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r_hcontainer.hpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r_hcontainer.hpp @@ -36,7 +36,7 @@ namespace ModuleSymmetry template // HContainer type void Symmetry_rotation::restore_HR( - const Symmetry& symm, const Atom* atoms, const Statistics& st, const char mode, + const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, const hamilt::HContainer& HR_irreduceble, hamilt::HContainer& HR_rotated)const { @@ -59,7 +59,7 @@ namespace ModuleSymmetry const int R_hc = ap_hc.find_R(R[0], R[1], R[2]); if (R_hc < 0) continue; TR* ijR_ptr = ap_hc.get_pointer(R_hc); - rotate_atompair_parallel(irijR_ptr, isym, atoms, st, irap, ap, mode, *irap_hc.get_paraV(), ijR_ptr); + rotate_atompair_parallel(irijR_ptr, isym, atoms, ucell, irap, ap, mode, *irap_hc.get_paraV(), ijR_ptr); } ModuleBase::timer::end("Symmetry_rotation", "restore_HR"); } @@ -77,13 +77,13 @@ namespace ModuleSymmetry obj[starti + i + (startj + j) * nr] = block(j, i); }; template - void Symmetry_rotation::rotate_atompair_parallel(const TR* Alocal_in, const int isym, const Atom* atoms, const Statistics& st, + void Symmetry_rotation::rotate_atompair_parallel(const TR* Alocal_in, const int isym, const Atom* atoms, const UnitCell& ucell, const Tap& ap_in, const Tap& ap_out, const char mode, const Parallel_Orbitals& pv, TR* Alocal_out, const bool output)const { // all the matrices are row-major (col-contiguous) int iat1 = ap_in.first, iat2 = ap_in.second; - int it1 = st.iat2it[iat1], it2 = st.iat2it[iat2]; - int ia1 = st.iat2ia[iat1], ia2 = st.iat2ia[iat2]; + int it1 = ucell.iat2it[iat1], it2 = ucell.iat2it[iat2]; + int ia1 = ucell.iat2ia[iat1], ia2 = ucell.iat2ia[iat2]; // contruct T matrix std::vector T1, T2; auto set_rotation_matrix = [&](const int& it, const int& ia, std::vector& T)->int @@ -149,8 +149,8 @@ namespace ModuleSymmetry // copy back to Alocal_out iat1 = ap_out.first, iat2 = ap_out.second; - it1 = st.iat2it[iat1], it2 = st.iat2it[iat2]; - ia1 = st.iat2ia[iat1], ia2 = st.iat2ia[iat2]; + it1 = ucell.iat2it[iat1], it2 = ucell.iat2it[iat2]; + ia1 = ucell.iat2ia[iat1], ia2 = ucell.iat2ia[iat2]; abr1 = pv.atom_begin_row[iat1], abc2 = pv.atom_begin_col[iat2]; iw1start = atoms[it1].stapos_wf + ia1 * atoms[it1].nw; iw2start = atoms[it2].stapos_wf + ia2 * atoms[it2].nw; @@ -164,7 +164,7 @@ namespace ModuleSymmetry } template - void Symmetry_rotation::test_HR_rotation(const Symmetry& symm, const Atom* atoms, const Statistics& st, + void Symmetry_rotation::test_HR_rotation(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, const hamilt::HContainer& HR_full) { ModuleBase::TITLE("Symmetry_rotation", "test_HR_rotation"); @@ -210,11 +210,11 @@ namespace ModuleSymmetry //2. rotate hamilt::HContainer HR_rotated(HR_full); HR_rotated.set_zero(); - this->restore_HR(symm, atoms, st, mode, HR_irreducible, HR_rotated); + this->restore_HR(symm, atoms, ucell, mode, HR_irreducible, HR_rotated); //3. compare - for (int iat1 = 0;iat1 < st.nat;++iat1) + for (int iat1 = 0;iat1 < ucell.nat;++iat1) { - for (int iat2 = 0;iat2 < st.nat;++iat2) + for (int iat2 = 0;iat2 < ucell.nat;++iat2) { const hamilt::AtomPair& ap_full = HR_full.get_atom_pair(iat1, iat2); const hamilt::AtomPair& ap_rotated = HR_rotated.get_atom_pair(iat1, iat2); diff --git a/source/source_lcao/module_ri/rpa_lri.hpp b/source/source_lcao/module_ri/rpa_lri.hpp index 5f19a66120b..8c4e3485465 100644 --- a/source/source_lcao/module_ri/rpa_lri.hpp +++ b/source/source_lcao/module_ri/rpa_lri.hpp @@ -167,7 +167,7 @@ void RPA_LRI::cal_postSCF_exx(const module_dm::DensityMatrix { const std::array period = RI_Util::get_Born_vonKarmen_period(kv); const auto& Rs = RI_Util::get_Born_von_Karmen_cells(period); - symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st, Rs, period, ucell.lat, PARAM.globalv.global_out_dir); + symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell, Rs, period, ucell.lat, PARAM.globalv.global_out_dir); // set Lmax of the rotation matrices to max(l_ao, l_abf), to support rotation under ABF // NOTE: Using Exx_Abfs::Construct_Orbs::get_Lmax() to compute Lmax from the actual ABFs // instead of relying on exx_cut_coulomb->abfs_Lmax() (not yet initialized) or diff --git a/source/source_lcao/test/test_init_dm_from_file.cpp b/source/source_lcao/test/test_init_dm_from_file.cpp index fb2ab91d67d..99647533ab5 100644 --- a/source/source_lcao/test/test_init_dm_from_file.cpp +++ b/source/source_lcao/test/test_init_dm_from_file.cpp @@ -39,8 +39,8 @@ class InitDMFileTest : public testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.atoms[0].taud.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); diff --git a/source/source_lcao/test/test_output_hcontainer_consistency.cpp b/source/source_lcao/test/test_output_hcontainer_consistency.cpp index 075495d4f1b..bed577dc630 100644 --- a/source/source_lcao/test/test_output_hcontainer_consistency.cpp +++ b/source/source_lcao/test/test_output_hcontainer_consistency.cpp @@ -41,8 +41,8 @@ class OutputHContainerTest : public testing::Test ucell.ntype = 1; ucell.nat = test_size; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.atoms[0].taud.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); diff --git a/source/source_pw/module_dfpt/test/dfpt_pw_data_test.cpp b/source/source_pw/module_dfpt/test/dfpt_pw_data_test.cpp index dcfc50dfbd7..8e879a638b0 100644 --- a/source/source_pw/module_dfpt/test/dfpt_pw_data_test.cpp +++ b/source/source_pw/module_dfpt/test/dfpt_pw_data_test.cpp @@ -50,7 +50,7 @@ class DFPT_PW_DataTest : public DFPTStruTestFixture ofs_running.open("tmp_dfpt_qlist"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {2, 2, 2}, true); data.init(&qlist, 1, 2, 3, 0, 1, 1, nullptr); } diff --git a/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp b/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp index 881eca4d94f..30cc2e65b63 100644 --- a/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp +++ b/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp @@ -46,7 +46,7 @@ class DFPT_PWRunTest : public DFPTStruTestFixture ofs_running.open("tmp_dfpt_run"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); ucell.symm = symm; } diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_q0_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_q0_serial_test.cpp index d7f12199275..a02dd2a66c1 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_q0_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_q0_serial_test.cpp @@ -556,10 +556,8 @@ TEST_F(DFPTQ0SerialTest, StarRotationCyclicGroup) ucell_.atoms[0].taud[i] = ModuleBase::Vector3(t0[i], t0[(i + 1) % 3], t0[(i + 2) % 3]); ucell_.atoms[0].tau[i] = ucell_.atoms[0].taud[i] * a_; } - delete[] ucell_.iat2it; - delete[] ucell_.iat2ia; - ucell_.iat2it = new int[3]; - ucell_.iat2ia = new int[3]; + ucell_.iat2it.resize(3); + ucell_.iat2ia.resize(3); for (int i = 0; i < 3; ++i) { ucell_.iat2it[i] = 0; diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp index 13d843f72ed..dfbc4f70923 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp @@ -24,8 +24,8 @@ void DFPTSerialBase::SetUpCell() ucell_.tpiba = ModuleBase::TWO_PI / lat0_; ucell_.tpiba2 = ucell_.tpiba * ucell_.tpiba; ucell_.omega = a_ * a_ * a_ * lat0_ * lat0_ * lat0_; - ucell_.iat2it = new int[1]; - ucell_.iat2ia = new int[1]; + ucell_.iat2it.resize(1); + ucell_.iat2ia.resize(1); ucell_.iat2it[0] = 0; ucell_.iat2ia[0] = 0; MakeCoulombAtom(); @@ -61,10 +61,6 @@ void DFPTSerialBase::TearDown() { delete[] ucell_.atoms; ucell_.atoms = nullptr; - delete[] ucell_.iat2it; - ucell_.iat2it = nullptr; - delete[] ucell_.iat2ia; - ucell_.iat2ia = nullptr; } void DFPTSerialBase::MakeCoulombAtom() @@ -143,10 +139,8 @@ void DFPTSerialBase::MakeTwoAtomCell() at.ncpp.nh = 0; at.mass = (it == 0) ? 12.0 : 4.0; } - delete[] ucell_.iat2it; - delete[] ucell_.iat2ia; - ucell_.iat2it = new int[2]; - ucell_.iat2ia = new int[2]; + ucell_.iat2it.resize(2); + ucell_.iat2ia.resize(2); ucell_.iat2it[0] = 0; ucell_.iat2ia[0] = 0; ucell_.iat2it[1] = 1; diff --git a/source/source_pw/module_pwdft/force_pw.cpp b/source/source_pw/module_pwdft/force_pw.cpp index d37e2379954..71dcec97601 100644 --- a/source/source_pw/module_pwdft/force_pw.cpp +++ b/source/source_pw/module_pwdft/force_pw.cpp @@ -367,7 +367,7 @@ void Forces::cal_force_loc(const UnitCell& ucell, resmem_var_op()(forcelc_d, this->nat * 3); resmem_var_op()(vloc_d, vloc.nr * vloc.nc); - syncmem_int_h2d_op()(iat2it_d, ucell.iat2it, this->nat); + syncmem_int_h2d_op()(iat2it_d, ucell.iat2it.data(), this->nat); syncmem_int_h2d_op()(ig2gg_d, rho_basis->ig2igg, rho_basis->npw); syncmem_var_h2d_op()(gcar_d, gcar_h.data(), rho_basis->npw * 3); syncmem_var_h2d_op()(tau_d, tau_h.data(), this->nat * 3); @@ -554,7 +554,7 @@ void Forces::cal_force_ew(const UnitCell& ucell, resmem_complex_op()(aux_d, rho_basis->npw); resmem_var_op()(forceion_d, this->nat * 3); - syncmem_int_h2d_op()(iat2it_d, ucell.iat2it, this->nat); + syncmem_int_h2d_op()(iat2it_d, ucell.iat2it.data(), this->nat); syncmem_var_h2d_op()(gcar_d, gcar_h.data(), rho_basis->npw * 3); syncmem_var_h2d_op()(tau_d, tau_h.data(), this->nat * 3); syncmem_var_h2d_op()(it_fact_d, it_fact_h.data(), ucell.ntype); diff --git a/source/source_pw/module_pwdft/test/dftu_base_test.cpp b/source/source_pw/module_pwdft/test/dftu_base_test.cpp index 2281b7f2186..cb2c3dc3a0d 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -32,8 +32,8 @@ class DFTUBaseTest : public testing::Test ucell.ntype = 1; ucell.nat = 1; ucell.atoms = new Atom[ucell.ntype]; - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.atoms[0].tau.resize(ucell.nat); ucell.atoms[0].taud.resize(ucell.nat); ucell.itia2iat.create(ucell.ntype, ucell.nat); diff --git a/source/source_relax/unittests/for_test.h b/source/source_relax/unittests/for_test.h index 043c88ed470..698f57f3d4f 100644 --- a/source/source_relax/unittests/for_test.h +++ b/source/source_relax/unittests/for_test.h @@ -24,11 +24,6 @@ UnitCell::UnitCell() namax = 0; nwmax = 0; - iat2it = nullptr; - iat2ia = nullptr; - iwt2iat = nullptr; - iwt2iw = nullptr; - itia2iat.create(1, 1); latvec = ModuleBase::Matrix3(); diff --git a/source/source_relax/unittests/test_ions_move_bfgs2.cpp b/source/source_relax/unittests/test_ions_move_bfgs2.cpp index 76dad4a9f8b..8d176b780f7 100644 --- a/source/source_relax/unittests/test_ions_move_bfgs2.cpp +++ b/source/source_relax/unittests/test_ions_move_bfgs2.cpp @@ -226,8 +226,8 @@ TEST_F(BFGSTest, UpdateUsesAbsoluteDisplacementThreshold) ucell.atoms = new Atom[ucell.ntype]; ucell.atoms[0].na = 1; ucell.atoms[0].mbl = std::vector>(1, {1, 1, 1}); - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); ucell.iat2it[0] = 0; ucell.iat2ia[0] = 0; @@ -266,8 +266,8 @@ TEST_F(BFGSTest, GetPosAndPostaud) ucell.atoms[0].tau[1].x = 1.0; ucell.atoms[0].tau[1].y = 0.0; ucell.atoms[0].tau[1].z = 0.0; // allocate mapping arrays - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); int k = 0; for (int it = 0; it < ucell.ntype; ++it) { for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { @@ -301,8 +301,8 @@ TEST_F(BFGSTest, CalculateLargestGrad) ucell.atoms[0].mbl = std::vector>(2, {1, 1, 1}); // mapping arrays - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); int k = 0; for (int it = 0; it < ucell.ntype; ++it) { for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { @@ -343,8 +343,8 @@ TEST_F(BFGSTest, RelaxStepBasic) ucell.atoms[0].tau = std::vector>(2); ucell.atoms[0].taud = std::vector>(2); ucell.atoms[0].mbl = std::vector>(2, {1, 1, 1}); - ucell.iat2it = new int[ucell.nat]; - ucell.iat2ia = new int[ucell.nat]; + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); int k = 0; for (int it = 0; it < ucell.ntype; ++it) { for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { diff --git a/source/source_relax/unittests/test_ions_move_lbfgs.cpp b/source/source_relax/unittests/test_ions_move_lbfgs.cpp index e7f882719b7..93275d92b9e 100644 --- a/source/source_relax/unittests/test_ions_move_lbfgs.cpp +++ b/source/source_relax/unittests/test_ions_move_lbfgs.cpp @@ -38,10 +38,6 @@ UnitCell::UnitCell() ntype = 1; nat = 1; - iat2it = nullptr; - iat2ia = nullptr; - iwt2iat = nullptr; - iwt2iw = nullptr; itia2iat.create(1, 1); atoms = new Atom[ntype]; diff --git a/source/source_relax/unittests/test_relax_nsync.cpp b/source/source_relax/unittests/test_relax_nsync.cpp index 82a900e204c..823fa901649 100644 --- a/source/source_relax/unittests/test_relax_nsync.cpp +++ b/source/source_relax/unittests/test_relax_nsync.cpp @@ -49,10 +49,6 @@ UnitCell::UnitCell() ntype = 1; nat = 1; - iat2it = nullptr; - iat2ia = nullptr; - iwt2iat = nullptr; - iwt2iw = nullptr; itia2iat.create(1, 1); atoms = new Atom[ntype]; @@ -181,13 +177,10 @@ class IonCellOptimizerTest : public ::testing::Test ucell.lat_axis_free[0] = 1; ucell.lat_axis_free[1] = 1; ucell.lat_axis_free[2] = 1; - // The mock UnitCell leaves iat2it/iat2ia null; the L-BFGS update loop - // dereferences them from the second iteration on, so give the single - // atom a valid type/index. Ownership passes to UnitCell's internal - // Statistics member, whose destructor releases them -- do NOT delete - // them here. - ucell.iat2it = new int[natom]; - ucell.iat2ia = new int[natom]; + // The L-BFGS update loop dereferences iat2it/iat2ia from the second + // iteration on, so give the single atom a valid type/index. + ucell.iat2it.resize(natom); + ucell.iat2ia.resize(natom); ucell.iat2it[0] = 0; ucell.iat2ia[0] = 0; } diff --git a/source/source_relax/unittests/test_relax_sync.cpp b/source/source_relax/unittests/test_relax_sync.cpp index b65d8bd9285..e645ed92f6e 100644 --- a/source/source_relax/unittests/test_relax_sync.cpp +++ b/source/source_relax/unittests/test_relax_sync.cpp @@ -43,8 +43,8 @@ class Test_SETGRAD : public testing::Test ucell.omega = 1.0; ucell.lat0 = 1.0; - ucell.iat2it = new int[nat]; - ucell.iat2ia = new int[nat]; + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); ucell.atoms[0].mbl.resize(nat); ucell.atoms[0].taud.resize(nat); ucell.atoms[0].tau.resize(nat); @@ -291,8 +291,8 @@ class Test_RELAX : public testing::Test ucell.omega = 452.590903143121; ucell.lat0 = 1.8897259886; - ucell.iat2it = new int[nat]; - ucell.iat2ia = new int[nat]; + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); ucell.iat2it[0] = 0; ucell.iat2it[1] = 1; ucell.iat2it[2] = 2; @@ -409,8 +409,8 @@ TEST(RelaxSyncSummary, ConvergedPrintsSummary) ucell.atoms[0].label = "Si"; ucell.omega = 1.0; ucell.lat0 = 1.0; - ucell.iat2it = new int[nat]; - ucell.iat2ia = new int[nat]; + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); ucell.iat2it[0] = 0; ucell.iat2ia[0] = 0; ucell.atoms[0].mbl.resize(nat); @@ -469,8 +469,8 @@ TEST(RelaxSyncSummary, TwoStepConvergence) ucell.atoms[0].label = "Si"; ucell.omega = 1.0; ucell.lat0 = 1.0; - ucell.iat2it = new int[nat]; - ucell.iat2ia = new int[nat]; + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); ucell.iat2it[0] = 0; ucell.iat2ia[0] = 0; ucell.atoms[0].mbl.resize(nat); @@ -539,8 +539,8 @@ TEST(RelaxSyncSummary, CellRelaxConvergedPrintsStressHistory) ucell.atoms[0].label = "Si"; ucell.omega = 1.0; ucell.lat0 = 1.0; - ucell.iat2it = new int[nat]; - ucell.iat2ia = new int[nat]; + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); ucell.iat2it[0] = 0; ucell.iat2ia[0] = 0; ucell.atoms[0].mbl.resize(nat); @@ -614,8 +614,8 @@ void make_two_type_cell(UnitCell& ucell, const int mbl_flat[9]) ucell.atoms[t].vel.resize(na); } // iat: 0 -> (type0, ia0); 1 -> (type1, ia0); 2 -> (type1, ia1) - ucell.iat2it = new int[3]; - ucell.iat2ia = new int[3]; + ucell.iat2it.resize(3); + ucell.iat2ia.resize(3); ucell.iat2it[0] = 0; ucell.iat2it[1] = 1; ucell.iat2it[2] = 1; ucell.iat2ia[0] = 0; ucell.iat2ia[1] = 0; ucell.iat2ia[2] = 1; @@ -640,11 +640,6 @@ void make_two_type_cell(UnitCell& ucell, const int mbl_flat[9]) void free_two_type_cell(UnitCell& ucell) { delete[] ucell.atoms; - // iat2it / iat2ia are owned by the mock; release explicitly. - delete[] ucell.iat2it; - delete[] ucell.iat2ia; - ucell.iat2it = nullptr; - ucell.iat2ia = nullptr; } } // namespace From 6c7e0998b7dcd1ef3870f9299da2630e72eacd6f Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 6 Oct 2026 09:17:52 +0800 Subject: [PATCH 02/15] Fix GPU build: use .data() for iat2it cudaMemcpy in gint_gpu_vars The Statistics refactor converted UnitCell::iat2it from int* to std::vector. The CUDA-only translation unit gint_gpu_vars.cpp was not compiled by the local CPU build, so this call site was missed: cudaMemcpy expects a const void* source but received the vector object. Pass ucell.iat2it.data() instead, matching the force_pw.cpp device-upload sites. --- source/source_hamilt/module_gint/kernel/gint_gpu_vars.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/source_hamilt/module_gint/kernel/gint_gpu_vars.cpp b/source/source_hamilt/module_gint/kernel/gint_gpu_vars.cpp index 1b893b36c40..beba4946e99 100644 --- a/source/source_hamilt/module_gint/kernel/gint_gpu_vars.cpp +++ b/source/source_hamilt/module_gint/kernel/gint_gpu_vars.cpp @@ -96,7 +96,7 @@ GintGpuVars::GintGpuVars(std::shared_ptr biggrid_info, CHECK_CUDA(cudaMemcpy(mgrids_pos_d, mgrids_pos_h.data(), sizeof(double3) * mgrid_num, cudaMemcpyHostToDevice)); CHECK_CUDA(cudaMalloc((void**)&iat2it_d, sizeof(int) * ucell.nat)); - CHECK_CUDA(cudaMemcpy(iat2it_d, ucell.iat2it, sizeof(int) * ucell.nat, cudaMemcpyHostToDevice)); + CHECK_CUDA(cudaMemcpy(iat2it_d, ucell.iat2it.data(), sizeof(int) * ucell.nat, cudaMemcpyHostToDevice)); } GintGpuVars::~GintGpuVars() From a7c403dfca01ca38dbee02563cd29167b713f633 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 9 Oct 2026 13:09:14 +0800 Subject: [PATCH 03/15] Delete UnitCell move operations to avoid double free The defaulted move constructor shallow-copied the owning raw pointer 'atoms' together with 'set_atom_flag' without resetting the source, so both the moved-to and moved-from objects ran delete[] on the same Atom array. The defaulted move assignment was implicitly deleted because the reference aliases into 'lat' (Coordinate, a1, latvec, ...) cannot be re-bound, triggering -Wdefaulted-function-deleted. Explicitly delete both move operations so any accidental move fails at compile time; correct move semantics require an owning container for 'atoms' and removal of the reference aliases in later refactoring steps. --- source/source_cell/unitcell.h | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/source/source_cell/unitcell.h b/source/source_cell/unitcell.h index 32588c86e3c..de3c04b8e02 100644 --- a/source/source_cell/unitcell.h +++ b/source/source_cell/unitcell.h @@ -18,13 +18,17 @@ class UnitCell : public BaseCell { UnitCell(); ~UnitCell(); - /// UnitCell owns dynamically-sized state (index maps, atom arrays, MPI-broadcast - /// buffers). Copying is disabled to make accidental shallow-copy bugs fail at - /// compile time; moving is allowed. + /// Copy and move are both disabled. The defaulted move constructor + /// shallow-copies the owning pointer 'atoms' and 'set_atom_flag' without + /// resetting the source (double free), and the reference aliases into + /// 'lat' (Coordinate, a1, latvec, ...) make the defaulted move assignment + /// implicitly deleted (-Wdefaulted-function-deleted). Correct move + /// semantics require an owning container for 'atoms' and removing the + /// aliases; until then any accidental move must fail at compile time. UnitCell(const UnitCell&) = delete; UnitCell& operator=(const UnitCell&) = delete; - UnitCell(UnitCell&&) = default; - UnitCell& operator=(UnitCell&&) = default; + UnitCell(UnitCell&&) = delete; + UnitCell& operator=(UnitCell&&) = delete; /// @brief Initialize basic cell parameters (latname, ntype, lmaxmax, init_vel) /// from INPUT and parse fixed_axes into lat_axis_free flags. From ba17904ae22baf10f2156339026af46cddc713e9 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 9 Oct 2026 13:21:54 +0800 Subject: [PATCH 04/15] Add compile-time guard for UnitCell non-copyable/non-movable contract Add static_asserts on std::is_copy/move_constructible/assignable so that re-enabling either operation during future UnitCell refactoring fails the translation unit at compile time instead of reintroducing the double free and dangling 'lat' reference aliases. --- source/source_cell/test/unitcell_test.cpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/source/source_cell/test/unitcell_test.cpp b/source/source_cell/test/unitcell_test.cpp index 9e8b1883cde..edeef913aa0 100644 --- a/source/source_cell/test/unitcell_test.cpp +++ b/source/source_cell/test/unitcell_test.cpp @@ -15,6 +15,7 @@ #include "source_cell/update_cell.h" #include "source_cell/read_stru.h" #include +#include #include #include @@ -133,6 +134,22 @@ class UcellTest : public ::testing::Test std::string output; }; +/// Compile-time guard: UnitCell owns the raw pointer 'atoms' tracked by +/// 'set_atom_flag' and exposes reference aliases into its 'lat' member, so +/// neither copy nor move can be implemented correctly until that ownership +/// model is refactored. Re-enabling either operation must fail this TU. +TEST(UnitCellTypeTraits, NotCopyableOrMovable) +{ + static_assert(!std::is_copy_constructible::value, + "UnitCell must not be copy constructible"); + static_assert(!std::is_copy_assignable::value, + "UnitCell must not be copy assignable"); + static_assert(!std::is_move_constructible::value, + "UnitCell must not be move constructible"); + static_assert(!std::is_move_assignable::value, + "UnitCell must not be move assignable"); +} + using UcellDeathTest = UcellTest; TEST_F(UcellTest, Constructor) From c1e17f1ce5e970821fa2ad102fb5a528948483d3 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 9 Oct 2026 13:37:57 +0800 Subject: [PATCH 05/15] update AGENTS.md --- AGENTS.md | 5 ----- 1 file changed, 5 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 2e1fab066c3..3a65ba552c6 100644 --- a/AGENTS.md +++ b/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. -- `UnitCell` is non-copyable (copy constructor/assignment deleted) and movable; - its index maps (`iat2it`, `iat2ia`, `iwt2iat`, `iwt2iw`) are - `std::vector` members owned directly by `UnitCell`. Size them with - `resize()`/`assign()`; never `new[]`/`delete[]` them. Do not reintroduce a - nested data struct re-exported through reference members. ## Review And Exception Flow From da443daebdcb795e58d36e6917223e47b0d1b17b Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 9 Oct 2026 13:53:29 +0800 Subject: [PATCH 06/15] Decouple Symmetry::analy_sys from UnitCell Replace the const UnitCell& parameter of analy_sys (and its three private helpers is_all_movable, analyze_magnetic_group, analyze_magnetic_group_nspin4) with only the index data they actually read: nat, ntype, iat2it, iat2ia, itia2iat. This breaks the UnitCell <-> Symmetry type-level dependency cycle: symmetry.h no longer names UnitCell at all, and the three implementation files drop their unitcell.h includes. The two other classes in the module that still need UnitCell (Irreducible_Sector, Symmetry_rotation_k) carry their own forward declaration now. All production and test call sites pass the five data members explicitly; behavior is unchanged. --- .../module_symmetry/irreducible_sector.h | 2 ++ .../module_symmetry/symm_analysis.cpp | 18 ++++++---- .../module_symmetry/symm_magnetic.cpp | 32 ++++++++++------- .../module_symmetry/symm_rotation_k.h | 2 ++ .../source_cell/module_symmetry/symmetry.cpp | 10 +++--- source/source_cell/module_symmetry/symmetry.h | 35 ++++++++++++++----- .../test/symm_test_analysis.cpp | 8 +++-- .../module_symmetry/test/symm_test_symtrz.cpp | 8 +++-- source/source_cell/test/klist_test.cpp | 12 +++++-- source/source_cell/test/klist_test_para.cpp | 8 +++-- source/source_cell/test/little_group_test.cpp | 4 ++- source/source_cell/test/qlist_test.cpp | 28 +++++++++++---- source/source_esolver/esolver_fp.cpp | 8 +++-- source/source_esolver/esolver_gets.cpp | 6 +++- .../source_esolver/esolver_lr_lcao_tddft.cpp | 4 ++- .../module_dfpt/test/dfpt_pw_data_test.cpp | 4 ++- .../module_dfpt/test/dfpt_pw_run_test.cpp | 4 ++- 17 files changed, 136 insertions(+), 57 deletions(-) diff --git a/source/source_cell/module_symmetry/irreducible_sector.h b/source/source_cell/module_symmetry/irreducible_sector.h index e01e1e7fa88..7e3d42ea672 100644 --- a/source/source_cell/module_symmetry/irreducible_sector.h +++ b/source/source_cell/module_symmetry/irreducible_sector.h @@ -9,6 +9,8 @@ #include "source_cell/module_symmetry/symmetry.h" #include "source_cell/klist.h" +class UnitCell; + namespace ModuleSymmetry { using Tap = std::pair; diff --git a/source/source_cell/module_symmetry/symm_analysis.cpp b/source/source_cell/module_symmetry/symm_analysis.cpp index 8c8d16c0fdb..a70db80f524 100644 --- a/source/source_cell/module_symmetry/symm_analysis.cpp +++ b/source/source_cell/module_symmetry/symm_analysis.cpp @@ -1,10 +1,14 @@ #include "symmetry.h" -#include "source_cell/unitcell.h" #include "source_base/output.h" using namespace ModuleSymmetry; -void Symmetry::analy_sys(const Lattice& lat, const UnitCell& ucell, Atom* atoms, std::ofstream& ofs_running, +void Symmetry::analy_sys(const Lattice& lat, Atom* atoms, + const int nat, const int ntype, + const std::vector& iat2it, + const std::vector& 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) { @@ -31,9 +35,9 @@ void Symmetry::analy_sys(const Lattice& lat, const UnitCell& ucell, Atom* atoms, // 1. copy data and allocate memory // -------------------------------- // number of total atoms - this->nat = ucell.nat; + this->nat = nat; // number of atom species - this->ntype = ucell.ntype; + this->ntype = ntype; assert(ntype>0); @@ -123,7 +127,7 @@ void Symmetry::analy_sys(const Lattice& lat, const UnitCell& ucell, Atom* atoms, if (!pricell_loop && nspin == 2) { - this->analyze_magnetic_group(atoms, ucell, nrot_out, nrotk_out); + this->analyze_magnetic_group(atoms, iat2it, iat2ia, itia2iat, nrot_out, nrotk_out); } else { @@ -295,13 +299,13 @@ void Symmetry::analy_sys(const Lattice& lat, const UnitCell& ucell, Atom* atoms, // the magnetization (pseudovector), so they are not applied in k-reduction / density symmetrization. if (nspin == 4) { - this->analyze_magnetic_group_nspin4(atoms, ucell, 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, ucell); + 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 "; diff --git a/source/source_cell/module_symmetry/symm_magnetic.cpp b/source/source_cell/module_symmetry/symm_magnetic.cpp index 60f314616ee..1c495ac93a8 100644 --- a/source/source_cell/module_symmetry/symm_magnetic.cpp +++ b/source/source_cell/module_symmetry/symm_magnetic.cpp @@ -3,12 +3,15 @@ using namespace ModuleSymmetry; #include "symm_rot_spin.h" #include "source_base/global_variable.h" -#include "source_cell/unitcell.h" #include #include -void Symmetry::analyze_magnetic_group(const Atom* atoms, const UnitCell& ucell, int& nrot_out, int& nrotk_out) +void Symmetry::analyze_magnetic_group(const Atom* atoms, + const std::vector& iat2it, + const std::vector& iat2ia, + const ModuleBase::IntArray& itia2iat, + int& nrot_out, int& nrotk_out) { // 1. classify atoms with different magmom // (use symmetry_prec to judge if two magmoms are the same) @@ -21,18 +24,18 @@ void Symmetry::analyze_magnetic_group(const Atom* atoms, const UnitCell& ucell, for (auto& mt : mag_type_atoms) { const int mag_iat = *mt.begin(); - const int mag_it = ucell.iat2it[mag_iat]; - const int mag_ia = ucell.iat2ia[mag_iat]; + const int mag_it = iat2it[mag_iat]; + const int mag_ia = iat2ia[mag_iat]; if (it == mag_it && this->equal(atoms[it].mag[ia], atoms[mag_it].mag[mag_ia])) { - mt.insert(ucell.itia2iat(it, ia)); + mt.insert(itia2iat(it, ia)); find = true; break; } } if (!find) { - mag_type_atoms.push_back(std::set({ ucell.itia2iat(it,ia) })); + mag_type_atoms.push_back(std::set({ itia2iat(it, ia) })); } } } @@ -59,7 +62,7 @@ void Symmetry::analyze_magnetic_group(const Atom* atoms, const UnitCell& ucell, { // this->newpos have been ordered by original structure(ntype, na), it cannot be directly used here. // we need to reset the calculate again the coordinate of the new structure. - const ModuleBase::Vector3 direct_tmp = atoms[ucell.iat2it[mag_iat]].tau[ucell.iat2ia[mag_iat]] * this->optlat.Inverse(); + const ModuleBase::Vector3 direct_tmp = atoms[iat2it[mag_iat]].tau[iat2ia[mag_iat]] * this->optlat.Inverse(); std::array direct = { direct_tmp.x, direct_tmp.y, direct_tmp.z }; for (int i = 0; i < 3; ++i) { @@ -79,7 +82,10 @@ void Symmetry::analyze_magnetic_group(const Atom* atoms, const UnitCell& ucell, } -void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const UnitCell& ucell, const ModuleBase::Matrix3& latvec) +void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, + const std::vector& iat2it, + const std::vector& iat2ia, + const ModuleBase::Matrix3& latvec) { // Restrict the space group to the unitary magnetic subgroup (nspin=4 / SOC): // operation g survives if it preserves the magnetic configuration as a pseudovector, @@ -98,7 +104,7 @@ void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const UnitCell& bool has_moment = false; for (int iat = 0; iat < this->nat && !has_moment; ++iat) { - const ModuleBase::Vector3& m = atoms[ucell.iat2it[iat]].m_loc_[ucell.iat2ia[iat]]; + const ModuleBase::Vector3& m = atoms[iat2it[iat]].m_loc_[iat2ia[iat]]; if (!this->equal(m.x, 0.0) || !this->equal(m.y, 0.0) || !this->equal(m.z, 0.0)) { has_moment = true; } } std::vector anti; // operations that REVERSE the moment: Theta*g is a symmetry @@ -110,11 +116,11 @@ void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const UnitCell& bool ok = true; for (int iat = 0; iat < this->nat && ok; ++iat) { - const ModuleBase::Vector3& m = atoms[ucell.iat2it[iat]].m_loc_[ucell.iat2ia[iat]]; + const ModuleBase::Vector3& m = atoms[iat2it[iat]].m_loc_[iat2ia[iat]]; // pseudovector-rotated moment W*m (column-vector convention: m'^i = W_ij m^j) const ModuleBase::Vector3& mrot = W * m; const int jat = this->get_rotated_atom(isym, iat); - const ModuleBase::Vector3& mj = atoms[ucell.iat2it[jat]].m_loc_[ucell.iat2ia[jat]]; + const ModuleBase::Vector3& mj = atoms[iat2it[jat]].m_loc_[iat2ia[jat]]; if (!this->equal(mrot.x, mj.x) || !this->equal(mrot.y, mj.y) || !this->equal(mrot.z, mj.z)) { ok = false; } } if (ok) @@ -130,10 +136,10 @@ void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const UnitCell& bool anti_ok = true; for (int iat = 0; iat < this->nat && anti_ok; ++iat) { - const ModuleBase::Vector3& m = atoms[ucell.iat2it[iat]].m_loc_[ucell.iat2ia[iat]]; + const ModuleBase::Vector3& m = atoms[iat2it[iat]].m_loc_[iat2ia[iat]]; const ModuleBase::Vector3& mrot = W * m; const int jat = this->get_rotated_atom(isym, iat); - const ModuleBase::Vector3& mj = atoms[ucell.iat2it[jat]].m_loc_[ucell.iat2ia[jat]]; + const ModuleBase::Vector3& mj = atoms[iat2it[jat]].m_loc_[iat2ia[jat]]; if (!this->equal(mrot.x, -mj.x) || !this->equal(mrot.y, -mj.y) || !this->equal(mrot.z, -mj.z)) { anti_ok = false; } } if (anti_ok) { anti.push_back(isym); } diff --git a/source/source_cell/module_symmetry/symm_rotation_k.h b/source/source_cell/module_symmetry/symm_rotation_k.h index d261dbed044..e85ddb21da5 100644 --- a/source/source_cell/module_symmetry/symm_rotation_k.h +++ b/source/source_cell/module_symmetry/symm_rotation_k.h @@ -5,6 +5,8 @@ #include "source_base/complexmatrix.h" #include "source_cell/module_symmetry/symm_rot_spin.h" +class UnitCell; + namespace ModuleSymmetry { /// @brief k-space AO-representation symmetry restoration: reconstructs D(k) at every diff --git a/source/source_cell/module_symmetry/symmetry.cpp b/source/source_cell/module_symmetry/symmetry.cpp index b8e99fd962d..09cf9836a49 100644 --- a/source/source_cell/module_symmetry/symmetry.cpp +++ b/source/source_cell/module_symmetry/symmetry.cpp @@ -1,7 +1,6 @@ #include #include #include "symmetry.h" -#include "source_cell/unitcell.h" namespace ModuleSymmetry { @@ -366,13 +365,14 @@ void Symmetry::get_optlat(ModuleBase::Vector3 &v1, ModuleBase::Vector3& iat2it, + const std::vector& iat2ia) const { bool all_mbl = true; - for (int iat = 0;iat < ucell.nat;++iat) + for (int iat = 0;iat < this->nat;++iat) { - int it = ucell.iat2it[iat]; - int ia = ucell.iat2ia[iat]; + int it = iat2it[iat]; + int ia = iat2ia[iat]; if (!atoms[it].mbl[ia].x || !atoms[it].mbl[ia].y || !atoms[it].mbl[ia].z) { all_mbl = false; diff --git a/source/source_cell/module_symmetry/symmetry.h b/source/source_cell/module_symmetry/symmetry.h index 0cce7004599..9a44d720533 100644 --- a/source/source_cell/module_symmetry/symmetry.h +++ b/source/source_cell/module_symmetry/symmetry.h @@ -12,8 +12,6 @@ #include "source_base/constants.h" #include "symm_basic.h" -class UnitCell; - namespace ModuleSymmetry { @@ -41,15 +39,24 @@ class Symmetry : public Symmetry_Basic /// @brief analyze the symmetry of the system /// @param lat structure of lattice - /// @param st /// @param atoms all atoms - /// @param ofs_running + /// @param nat total number of atoms + /// @param ntype number of atom species + /// @param iat2it map from global atom index iat to its species index it + /// @param iat2ia map from global atom index iat to its index ia within the species + /// @param itia2iat map from (it, ia) to the global atom index iat + /// @param ofs_running /// @param symmetry_prec precision for symmetry analysis /// @param nspin number of spin components /// @param calculation calculation type (scf, relax, cell-relax, etc.) /// @param cal_symm_repr control for symmetry representation output [0]=flag, [1]=precision /// get the symmetry information of the system, gmatries (rotation 3*3 matrixs), gtrans (transfer a collections vector3), etc. - void analy_sys(const Lattice& lat, const UnitCell& ucell, Atom* atoms, std::ofstream& ofs_running, + void analy_sys(const Lattice& lat, Atom* atoms, + const int nat, const int ntype, + const std::vector& iat2it, + const std::vector& 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); @@ -351,7 +358,10 @@ class Symmetry : public Symmetry_Basic 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 UnitCell& ucell)const; + /// @param iat2it map from global atom index iat to its species index it + /// @param iat2ia map from global atom index iat to its index ia within the species + bool is_all_movable(const Atom* atoms, const std::vector& iat2it, + const std::vector& iat2ia) const; // to be called in lattice_type void get_shortest_latvec(ModuleBase::Vector3 &a1, @@ -366,10 +376,14 @@ class Symmetry : public Symmetry_Basic /// If not all the same, primitive cells should not be looped in rhog_symmetry. bool magmom_same_check(const Atom* atoms)const; - /// Analyze magnetic group without time-reversal symmetry + /// Analyze magnetic group without time-reversal symmetry /// (because currently the charge density symmetrization does not support it) /// Method: treat atoms with different magmom as atoms of different type - void analyze_magnetic_group(const Atom* atoms, const UnitCell& ucell, int& nrot_out, int& nrotk_out); + void analyze_magnetic_group(const Atom* atoms, + const std::vector& iat2it, + const std::vector& iat2ia, + const ModuleBase::IntArray& itia2iat, + int& nrot_out, int& nrotk_out); /// (nspin=4 / SOC) Restrict the already-built space group to the unitary magnetic /// subgroup: keep operation g only if it preserves the magnetization as a pseudovector, @@ -377,7 +391,10 @@ class Symmetry : public Symmetry_Basic /// that reverse the moment (which are only symmetries when combined with time reversal) /// from being applied in k-reduction and density symmetrization. /// Non-magnetic (m_i=0) keeps all operations. - void analyze_magnetic_group_nspin4(const Atom* atoms, const UnitCell& ucell, const ModuleBase::Matrix3& latvec); + void analyze_magnetic_group_nspin4(const Atom* atoms, + const std::vector& iat2it, + const std::vector& iat2ia, + const ModuleBase::Matrix3& latvec); }; /** diff --git a/source/source_cell/module_symmetry/test/symm_test_analysis.cpp b/source/source_cell/module_symmetry/test/symm_test_analysis.cpp index 1a2a47011ad..98f873e294a 100644 --- a/source/source_cell/module_symmetry/test/symm_test_analysis.cpp +++ b/source/source_cell/module_symmetry/test/symm_test_analysis.cpp @@ -46,7 +46,9 @@ TEST_F(SymmetryTest, AnalySys) ModuleSymmetry::Symmetry symm; construct_ucell(stru_lib[stru]); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); //1. ibrav std::string ref_point_group = stru_lib[stru].point_group; @@ -250,7 +252,9 @@ TEST_F(SymmetryTest, SG_Pricell) symm.epsilon = 1e-5; construct_ucell(supercell_lib[stru]); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); std::string ref_point_group = supercell_lib[stru].point_group; std::string cal_point_group = symm.pgname; diff --git a/source/source_cell/module_symmetry/test/symm_test_symtrz.cpp b/source/source_cell/module_symmetry/test/symm_test_symtrz.cpp index 9911fa94755..845a446a880 100644 --- a/source/source_cell/module_symmetry/test/symm_test_symtrz.cpp +++ b/source/source_cell/module_symmetry/test/symm_test_symtrz.cpp @@ -71,7 +71,9 @@ TEST_F(SymmetryTest, ForceSymmetry) ModuleSymmetry::Symmetry symm; construct_ucell(supercell_lib[stru]); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); ModuleBase::matrix force(ucell.nat, 3, true); //generate random number for force and restrict to [-100,100) @@ -103,7 +105,9 @@ TEST_F(SymmetryTest, StressSymmetry) ModuleSymmetry::Symmetry symm; construct_ucell(supercell_lib[stru]); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); ModuleBase::matrix stress(3, 3, true); //generate random number for stress and restrict to [-1e5,1e5) diff --git a/source/source_cell/test/klist_test.cpp b/source/source_cell/test/klist_test.cpp index ecedde8484d..d512debf5c9 100644 --- a/source/source_cell/test/klist_test.cpp +++ b/source/source_cell/test/klist_test.cpp @@ -844,7 +844,9 @@ TEST_F(KlistTest, IbzKpoint) construct_ucell(stru_lib[0]); ofs_running.open("tmp_klist_3"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); std::string k_file = "./support/KPT1"; kv->set_spin_mult(1); kv->read_kpoints_for_testing(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); @@ -871,7 +873,9 @@ TEST_F(KlistTest, IbzKpointIsMP) construct_ucell(stru_lib[0]); ofs_running.open("tmp_klist_4"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); std::string k_file = "./support/KPT1"; kv->set_spin_mult(1); kv->read_kpoints_for_testing(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); @@ -899,7 +903,9 @@ TEST_F(KlistTest, IbzKpointCustomWeights) construct_ucell(stru_lib[0]); ofs_running.open("tmp_klist_custom_weights"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); // Test 1: Non-MP k-points with uniform weights (KPT4) { diff --git a/source/source_cell/test/klist_test_para.cpp b/source/source_cell/test/klist_test_para.cpp index a4fca253cd3..dfe44ad42e4 100644 --- a/source/source_cell/test/klist_test_para.cpp +++ b/source/source_cell/test/klist_test_para.cpp @@ -178,7 +178,9 @@ TEST_F(KlistParaTest, Set) GlobalV::ofs_running.open("tmp_klist_5"); } const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); // read KPT std::string k_file = "./support/KPT1"; // note: do NOT pre-set kv->spin_mult here; set() takes the physical @@ -302,7 +304,9 @@ TEST_F(KlistParaTest, SetAfterVC) GlobalV::ofs_running.open("tmp_klist_6"); } const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); // read KPT std::string k_file = "./support/KPT1"; // note: do NOT pre-set kv->spin_mult here; set() takes the physical diff --git a/source/source_cell/test/little_group_test.cpp b/source/source_cell/test/little_group_test.cpp index a5c0454a7ed..1f0c6c6b2b3 100644 --- a/source/source_cell/test/little_group_test.cpp +++ b/source/source_cell/test/little_group_test.cpp @@ -120,7 +120,9 @@ class LittleGroupTest : public testing::Test } std::ofstream ofs_running("tmp_little_group"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); } void TearDown() override diff --git a/source/source_cell/test/qlist_test.cpp b/source/source_cell/test/qlist_test.cpp index 959fa254725..1344a692ebc 100644 --- a/source/source_cell/test/qlist_test.cpp +++ b/source/source_cell/test/qlist_test.cpp @@ -173,7 +173,9 @@ TEST_F(QListTest, GenerateMeshFullSymmetry) ofs_running.open("tmp_qlist_1"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {8, 8, 8}, true); @@ -211,7 +213,9 @@ TEST_F(QListTest, GenerateMeshSmallGrid) ofs_running.open("tmp_qlist_2"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {2, 2, 2}, true); @@ -235,7 +239,9 @@ TEST_F(QListTest, GammaOnlyGrid) ofs_running.open("tmp_qlist_3"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {1, 1, 1}, true); @@ -255,7 +261,9 @@ TEST_F(QListTest, IrrepPlaceholder) ofs_running.open("tmp_qlist_4"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {2, 2, 2}, true); @@ -282,7 +290,9 @@ TEST_F(QListTest, CartesianCoordinatesComputed) ofs_running.open("tmp_qlist_cart"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {2, 2, 2}, true); @@ -312,7 +322,9 @@ TEST_F(QListTest, UseIrrepsSwitch) ofs_running.open("tmp_qlist_irreps"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); // use_irreps = false: the q mesh is still reduced, but no irrep data qlist.generate_mesh(ucell, symm, {2, 2, 2}, false); @@ -331,7 +343,9 @@ TEST_F(QListTest, PrintQlists) ofs_running.open("tmp_qlist_print"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {1, 1, 1}, false); diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 07bf211b482..962c4869118 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -119,7 +119,9 @@ void ESolver_FP::before_all_runners(BaseCell& basecell, const Input_para& inp) if (ModuleSymmetry::Symmetry::symm_flag == 1) { const int cal_symm_repr[2] = {this->inp_->cal_symm_repr[0], this->inp_->cal_symm_repr[1]}; - ucell.symm.analy_sys(ucell.lat, ucell, ucell.atoms, GlobalV::ofs_running, + ucell.symm.analy_sys(ucell.lat, ucell.atoms, + ucell.nat, ucell.ntype, ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + GlobalV::ofs_running, this->inp_->symmetry_prec, inp.nspin, this->inp_->calculation, cal_symm_repr); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } @@ -212,7 +214,9 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) if (ModuleSymmetry::Symmetry::symm_flag == 1) { const int cal_symm_repr[2] = {this->inp_->cal_symm_repr[0], this->inp_->cal_symm_repr[1]}; - ucell.symm.analy_sys(ucell.lat, ucell, ucell.atoms, GlobalV::ofs_running, + ucell.symm.analy_sys(ucell.lat, ucell.atoms, + ucell.nat, ucell.ntype, ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + GlobalV::ofs_running, this->inp_->symmetry_prec, this->inp_->nspin, this->inp_->calculation, cal_symm_repr); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } diff --git a/source/source_esolver/esolver_gets.cpp b/source/source_esolver/esolver_gets.cpp index 91c72cde48a..70dfa9f130c 100644 --- a/source/source_esolver/esolver_gets.cpp +++ b/source/source_esolver/esolver_gets.cpp @@ -69,8 +69,12 @@ void ESolver_GetS::before_all_runners(BaseCell& basecell, const Input_para& inp) { const int cal_symm_repr[2] = {this->inp_->cal_symm_repr[0], this->inp_->cal_symm_repr[1]}; ucell.symm.analy_sys(ucell.lat, - ucell, ucell.atoms, + ucell.nat, + ucell.ntype, + ucell.iat2it, + ucell.iat2ia, + ucell.itia2iat, GlobalV::ofs_running, this->inp_->symmetry_prec, inp.nspin, diff --git a/source/source_esolver/esolver_lr_lcao_tddft.cpp b/source/source_esolver/esolver_lr_lcao_tddft.cpp index 71ae7f1f4b7..bfb9bf2bb00 100644 --- a/source/source_esolver/esolver_lr_lcao_tddft.cpp +++ b/source/source_esolver/esolver_lr_lcao_tddft.cpp @@ -374,7 +374,9 @@ void ModuleESolver::ESolver_LR::initialize_from_unitcell_(UnitCell& ucell if (ModuleSymmetry::Symmetry::symm_flag == 1) { const int cal_symm_repr[2] = {this->inp_->cal_symm_repr[0], this->inp_->cal_symm_repr[1]}; - ucell.symm.analy_sys(ucell.lat, ucell, ucell.atoms, GlobalV::ofs_running, + ucell.symm.analy_sys(ucell.lat, ucell.atoms, + ucell.nat, ucell.ntype, ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + GlobalV::ofs_running, this->inp_->symmetry_prec, this->inp_->nspin, this->inp_->calculation, cal_symm_repr); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } diff --git a/source/source_pw/module_dfpt/test/dfpt_pw_data_test.cpp b/source/source_pw/module_dfpt/test/dfpt_pw_data_test.cpp index 8e879a638b0..3a89b83ba29 100644 --- a/source/source_pw/module_dfpt/test/dfpt_pw_data_test.cpp +++ b/source/source_pw/module_dfpt/test/dfpt_pw_data_test.cpp @@ -50,7 +50,9 @@ class DFPT_PW_DataTest : public DFPTStruTestFixture ofs_running.open("tmp_dfpt_qlist"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); qlist.generate_mesh(ucell, symm, {2, 2, 2}, true); data.init(&qlist, 1, 2, 3, 0, 1, 1, nullptr); } diff --git a/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp b/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp index 30cc2e65b63..c2b5646666f 100644 --- a/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp +++ b/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp @@ -46,7 +46,9 @@ class DFPT_PWRunTest : public DFPTStruTestFixture ofs_running.open("tmp_dfpt_run"); ModuleSymmetry::Symmetry symm; const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, "scf", cal_symm_repr); ucell.symm = symm; } From f2b7f5d015ef6b7c371097a3eaa3fd852e8e405a Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 9 Oct 2026 14:25:33 +0800 Subject: [PATCH 07/15] Decouple Irred_Sector and Symmetry_rotation_k from UnitCell Replace the const UnitCell& parameters of find_irred_sector, cal_return_lattice_all, gen_symm_bvk (Irred_Sector) and cal_Ms / contruct_2d_rot_mat_ao (Symmetry_rotation_k) with the data they actually read: nat, ntype, the iat2it/iat2ia maps, and explicit symm/atoms/lat/latvec/lmax arguments. This removes the possibility of passing mismatched atoms/lat alongside one cell's index maps and makes both headers free of any UnitCell name; the three implementation files drop their unitcell.h includes. Drop the unused get_aRb_direct helpers (no callers anywhere in the tree). Rename for brevity: Irreducible_Sector -> Irred_Sector, irreducible_sector.{h,cpp} -> irred_sec.{h,cpp} (matching the existing irred_sec_bvk.cpp and the irs_ member), find/get_irreducible_sector -> find/get_irred_sector, and gen_symmetry_BvK -> gen_symm_bvk. Comments, log strings and the output file name keep the standard "BvK" spelling. All three production call sites (DFTU x2, RDMFT) and the EXX info printer are updated; behavior is unchanged. --- .../module_symmetry/CMakeLists.txt | 2 +- .../{irreducible_sector.cpp => irred_sec.cpp} | 67 +++++++++---------- .../{irreducible_sector.h => irred_sec.h} | 23 +++---- .../module_symmetry/irred_sec_bvk.cpp | 25 +++---- .../module_symmetry/symm_rotation_k.cpp | 40 +++++------ .../module_symmetry/symm_rotation_k.h | 31 ++++++--- .../source_lcao/module_dftu/dftu_nao_occ.cpp | 7 +- .../source_lcao/module_dftu/dftu_nao_op.cpp | 8 ++- source/source_lcao/module_rdmft/rdmft.cpp | 8 ++- .../module_exx_symmetry/symm_rot_out.cpp | 2 +- 10 files changed, 117 insertions(+), 96 deletions(-) rename source/source_cell/module_symmetry/{irreducible_sector.cpp => irred_sec.cpp} (85%) rename source/source_cell/module_symmetry/{irreducible_sector.h => irred_sec.h} (87%) diff --git a/source/source_cell/module_symmetry/CMakeLists.txt b/source/source_cell/module_symmetry/CMakeLists.txt index 5b2a03cd967..cd6eec7056c 100644 --- a/source/source_cell/module_symmetry/CMakeLists.txt +++ b/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 ) diff --git a/source/source_cell/module_symmetry/irreducible_sector.cpp b/source/source_cell/module_symmetry/irred_sec.cpp similarity index 85% rename from source/source_cell/module_symmetry/irreducible_sector.cpp rename to source/source_cell/module_symmetry/irred_sec.cpp index 49b6897eb19..0eda5200c3d 100644 --- a/source/source_cell/module_symmetry/irreducible_sector.cpp +++ b/source/source_cell/module_symmetry/irred_sec.cpp @@ -1,5 +1,11 @@ -#include "source_cell/module_symmetry/irreducible_sector.h" -#include "source_cell/unitcell.h" +#include "source_cell/module_symmetry/irred_sec.h" +#include "source_base/global_variable.h" + +#include +#include +#include +#include + namespace ModuleSymmetry { // Raw-index dispatch shared by the real-space sector helpers, matching the convention used @@ -17,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(x + symm.epsilon) : static_cast(x - symm.epsilon); }; @@ -29,7 +35,7 @@ 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), @@ -37,19 +43,6 @@ namespace ModuleSymmetry 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 UnitCell& ucell, - const int iat1, const int iat2, const TCdouble& R, const char gauge)const - { - return TCdouble(atoms[ucell.iat2it[iat1]].taud[ucell.iat2ia[iat1]] - atoms[ucell.iat2it[iat2]].taud[ucell.iat2ia[iat2]]) + (gauge == 'L' ? R : TCdouble(-R)); - } - - TCdouble Irreducible_Sector::get_aRb_direct(const Atom* atoms, const UnitCell& ucell, - const int iat1, const int iat2, const TC& R, const char gauge) const - { - const TCdouble R_double(static_cast(R[0]), static_cast(R[1]), static_cast(R[2])); - return get_aRb_direct(atoms, ucell, iat1, iat2, R_double); - } - inline void output_return_lattice(const std::vector>& return_lattice) { std::cout << "return lattice:" << std::endl; @@ -65,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 { @@ -104,27 +97,28 @@ 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 UnitCell& ucell) + void Irred_Sector::cal_return_lattice_all(const Symmetry& symm, const Atom* atoms, + const int nat, const std::vector& iat2it, const std::vector& 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(ucell.nat, std::vector(symm.nrotk + symm.nrotk_anti)); - for (int iat1 = 0;iat1 < ucell.nat;++iat1) + this->return_lattice_.resize(nat, std::vector(symm.nrotk + symm.nrotk_anti)); + for (int iat1 = 0;iat1 < nat;++iat1) { - int it = ucell.iat2it[iat1]; - int ia1 = ucell.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 = ucell.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 = ucell.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]); } } @@ -132,7 +126,7 @@ namespace ModuleSymmetry // 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_) @@ -147,7 +141,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_) @@ -170,7 +164,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()) { @@ -184,20 +178,23 @@ namespace ModuleSymmetry } } - void Irreducible_Sector::find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const std::vector& 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& iat2it, const std::vector& iat2ia, + const std::vector& 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, ucell); + 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> apR_all; - for (int iat1 = 0;iat1 < ucell.nat; iat1++) - for (int iat2 = 0; iat2 < ucell.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); @@ -220,7 +217,7 @@ namespace ModuleSymmetry // get symmetry of BvK supercell if (this->isymbvk_to_isym_.empty()) - this->gen_symmetry_BvK(symm, atoms, lat, ucell, period); + this->gen_symm_bvk(symm, atoms, lat, nat, ntype, period); assert(!this->isymbvk_to_isym_.empty()); // std::vector in_2d_plain; // const bool judge_2d = (symm.real_brav == 4); @@ -275,7 +272,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(ucell.nat); + // this->output_full_map_to_irreducible_sector(nat); // this->output_sector_star(); this->write_irreducible_sector(output_dir); } diff --git a/source/source_cell/module_symmetry/irreducible_sector.h b/source/source_cell/module_symmetry/irred_sec.h similarity index 87% rename from source/source_cell/module_symmetry/irreducible_sector.h rename to source/source_cell/module_symmetry/irred_sec.h index 7e3d42ea672..2e8808e4ac0 100644 --- a/source/source_cell/module_symmetry/irreducible_sector.h +++ b/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 #include @@ -9,8 +9,6 @@ #include "source_cell/module_symmetry/symmetry.h" #include "source_cell/klist.h" -class UnitCell; - namespace ModuleSymmetry { using Tap = std::pair; @@ -18,7 +16,7 @@ namespace ModuleSymmetry using TapR = std::pair; using TCdouble = ModuleBase::Vector3; - class Irreducible_Sector + class Irred_Sector { public: struct ap_less_func @@ -63,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 UnitCell& ucell, + void find_irred_sector(const Symmetry& symm, const Atom* atoms, + const int nat, const int ntype, + const std::vector& iat2it, const std::vector& iat2ia, const std::vector& Rs, const TC& period, const Lattice& lat, const std::string& output_dir = ""); - const std::map>& get_irreducible_sector()const { return this->irreducible_sector_; } + const std::map>& get_irred_sector()const { return this->irreducible_sector_; } // const std::map>> convirt_irreducible_sector() {}; //-------------------------------------------------------------------------------- @@ -88,7 +88,8 @@ namespace ModuleSymmetry return this->return_lattice_[iat][isym]; } protected: - void cal_return_lattice_all(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell); + void cal_return_lattice_all(const Symmetry& symm, const Atom* atoms, + const int nat, const std::vector& iat2it, const std::vector& iat2ia); //-------------------------------------------------------------------------------- /// The sub functions to find irreducible sector: {abR} @@ -97,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 UnitCell& ucell, const int iat1, const int iat2, const TC& R, const char gauge = 'R')const; - TCdouble get_aRb_direct(const Atom* atoms, const UnitCell& ucell, 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; @@ -117,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 UnitCell& ucell, 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 in_plain(const ModuleSymmetry::Symmetry& symm, const ModuleBase::Matrix3& latvec)const; //-------------------------------------------------------------------------------- diff --git a/source/source_cell/module_symmetry/irred_sec_bvk.cpp b/source/source_cell/module_symmetry/irred_sec_bvk.cpp index 115a8bbb3e6..714e9cf5063 100644 --- a/source/source_cell/module_symmetry/irred_sec_bvk.cpp +++ b/source/source_cell/module_symmetry/irred_sec_bvk.cpp @@ -1,8 +1,8 @@ -#include "irreducible_sector.h" -#include "source_cell/unitcell.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; } @@ -44,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; @@ -56,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 UnitCell& ucell, 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& a1, const ModuleBase::Vector3& a2, const ModuleBase::Vector3& 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 @@ -79,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 = ucell.nat * bvk_min_period[0] * bvk_min_period[1] * bvk_min_period[2]; - std::vector bvk_na(ucell.ntype); - std::vector bvk_istart(ucell.ntype, 0); + const int bvk_nat = nat * bvk_min_period[0] * bvk_min_period[1] * bvk_min_period[2]; + std::vector bvk_na(ntype); + std::vector bvk_istart(ntype, 0); int bvk_itmin_start = 0, bvk_itmin_type = 0; - for (int it = 0;it < ucell.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]; @@ -104,7 +105,7 @@ namespace ModuleSymmetry s3 = a3 = lat.a3 * static_cast(bvk_min_period[2]); ModuleBase::Matrix3 bvk_min_lat = set_matrix3(s1, s2, s3); int at = 0; - for (int it = 0; it < ucell.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) { @@ -185,7 +186,7 @@ namespace ModuleSymmetry return; } - // std::vector Irreducible_Sector::in_plain(const ModuleSymmetry::Symmetry& symm, const ModuleBase::Matrix3& latvec)const + // std::vector 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 gmatc(symm.nrotk); diff --git a/source/source_cell/module_symmetry/symm_rotation_k.cpp b/source/source_cell/module_symmetry/symm_rotation_k.cpp index b56bbf4c1b8..aeee3c3763b 100644 --- a/source/source_cell/module_symmetry/symm_rotation_k.cpp +++ b/source/source_cell/module_symmetry/symm_rotation_k.cpp @@ -1,5 +1,4 @@ #include "symm_rotation_k.h" -#include "source_cell/unitcell.h" #include "source_base/constants.h" #include #include "source_base/parallel_reduce.h" @@ -23,29 +22,31 @@ namespace ModuleSymmetry return cells; } - void Symmetry_rotation_k::cal_Ms(const K_Vectors& kv, - const UnitCell& ucell, const Parallel_2D& pv, const int nspin) + void Symmetry_rotation_k::cal_Ms(const K_Vectors& kv, const Symmetry& symm, const Atom* atoms, + const ModuleBase::Matrix3& latvec, const int lmax, + const int nat, const std::vector& iat2it, const std::vector& iat2ia, + const Parallel_2D& pv, const int nspin) { ModuleBase::TITLE("Symmetry_rotation_k", "cal_Ms"); ModuleBase::timer::start("Symmetry_rotation_k", "cal_Ms"); this->nspin_ = nspin; - this->nsym_ = ucell.symm.nrotk; - this->nanti_ = ucell.symm.nrotk_anti; - this->magnetic_nspin4_ = ucell.symm.magnetic_nspin4; - this->eps_ = ucell.symm.epsilon; + this->nsym_ = symm.nrotk; + this->nanti_ = symm.nrotk_anti; + this->magnetic_nspin4_ = symm.magnetic_nspin4; + this->eps_ = symm.epsilon; if (this->irs_.invmap_.empty()) { - this->irs_.invmap_.resize(ucell.symm.nrotk); - ucell.symm.gmatrix_invmap(ucell.symm.gmatrix, ucell.symm.nrotk, this->irs_.invmap_.data()); + this->irs_.invmap_.resize(symm.nrotk); + symm.gmatrix_invmap(symm.gmatrix, symm.nrotk, this->irs_.invmap_.data()); } // 1. calculate the rotation matrix in real spherical harmonics representation for each symmetry operation: [T_l (isym)]_mm' const int nop_tot = this->nsym_ + this->nanti_; std::vector gmatc(nop_tot); - for (int i = 0;i < nsym_;++i) { gmatc[i] = this->irs_.direct_to_cartesian(ucell.symm.gmatrix[i], ucell.latvec); } + for (int i = 0;i < nsym_;++i) { gmatc[i] = this->irs_.direct_to_cartesian(symm.gmatrix[i], latvec); } for (int j = 0;j < this->nanti_;++j) - { gmatc[nsym_ + j] = this->irs_.direct_to_cartesian(ucell.symm.gmatrix_anti[j], ucell.latvec); } - this->cal_rotmat_Slm(gmatc.data(), std::max(this->abfs_Lmax_, ucell.lmax), nop_tot); + { gmatc[nsym_ + j] = this->irs_.direct_to_cartesian(symm.gmatrix_anti[j], latvec); } + this->cal_rotmat_Slm(gmatc.data(), std::max(this->abfs_Lmax_, lmax), nop_tot); // 1.5 (nspin=4) the SU(2) spin-1/2 rotation U(isym) for each symmetry operation. The AO // rotation matrix M becomes the spinor operator T(isym) (x) U(isym) so that the same @@ -101,7 +102,7 @@ namespace ModuleSymmetry } for (int op = 0; op < nsym_; ++op) { - const ModuleBase::Vector3 delta = kvec_d_ibz_global[ik_ibz] * ucell.symm.kgmatrix[op] - kvec_d_ibz_global[ik_ibz]; + const ModuleBase::Vector3 delta = kvec_d_ibz_global[ik_ibz] * symm.kgmatrix[op] - kvec_d_ibz_global[ik_ibz]; if (std::abs(delta.x - std::round(delta.x)) < this->eps_ && std::abs(delta.y - std::round(delta.y)) < this->eps_ && std::abs(delta.z - std::round(delta.z)) < this->eps_) @@ -113,7 +114,7 @@ namespace ModuleSymmetry for (const int op : needed) { this->Ms_[ik_ibz][op] = this->contruct_2d_rot_mat_ao( - ucell.symm, ucell.atoms, ucell, kvec_d_ibz_global[ik_ibz], op, pv, spin_U[op]); + symm, atoms, nat, iat2it, iat2ia, kvec_d_ibz_global[ik_ibz], op, pv, spin_U[op]); } } @@ -389,7 +390,8 @@ namespace ModuleSymmetry // 2d-block parallized rotation matrix in AO-representation, denoted as M. // finally we will use D(k)=M(R, k)^\dagger*D(Rk)*M(R, k) to D(k) from D(Rk) in cal_Ms. - std::vector> Symmetry_rotation_k::contruct_2d_rot_mat_ao(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, + std::vector> Symmetry_rotation_k::contruct_2d_rot_mat_ao(const Symmetry& symm, const Atom* atoms, + const int nat, const std::vector& iat2it, const std::vector& iat2ia, const TCdouble& kvec_d_ibz, int isym, const Parallel_2D& pv, const SpinRotation::Su2& spin_U) const { const bool soc = (this->nspin_ == 4); @@ -403,12 +405,12 @@ namespace ModuleSymmetry return (is < nrotk_u) ? symm.get_rotated_atom(is, iat) : symm.get_rotated_atom_anti(is - nrotk_u, iat); }; - for (int iat1 = 0;iat1 < ucell.nat;++iat1) + for (int iat1 = 0;iat1 < nat;++iat1) { - int it = ucell.iat2it[iat1]; // it1=it2 - int ia1 = ucell.iat2ia[iat1]; + int it = iat2it[iat1]; // it1=it2 + int ia1 = iat2ia[iat1]; int iat2 = rotated_atom(isym, iat1); //iat2=rot(iat1) - int ia2 = ucell.iat2ia[iat2]; + int ia2 = iat2ia[iat2]; // cal phase factor from return lattice: exp(-ik_ibz*O) double arg = -2 * ModuleBase::PI * kvec_d_ibz * this->irs_.return_lattice_[iat1][isym]; std::complexphase_factor = std::complex(std::cos(arg), std::sin(arg)); diff --git a/source/source_cell/module_symmetry/symm_rotation_k.h b/source/source_cell/module_symmetry/symm_rotation_k.h index e85ddb21da5..402b57d312c 100644 --- a/source/source_cell/module_symmetry/symm_rotation_k.h +++ b/source/source_cell/module_symmetry/symm_rotation_k.h @@ -1,12 +1,10 @@ #ifndef SYMM_ROTATION_K_H #define SYMM_ROTATION_K_H -#include "irreducible_sector.h" +#include "irred_sec.h" #include "source_base/parallel_2d.h" #include "source_base/complexmatrix.h" #include "source_cell/module_symmetry/symm_rot_spin.h" -class UnitCell; - namespace ModuleSymmetry { /// @brief k-space AO-representation symmetry restoration: reconstructs D(k) at every @@ -28,7 +26,7 @@ namespace ModuleSymmetry //-------------------------------------------------------------------------------- // getters - const std::map>& get_irreducible_sector()const { return this->irs_.get_irreducible_sector(); } + const std::map>& get_irred_sector()const { return this->irs_.get_irred_sector(); } TCdouble get_return_lattice(const Symmetry& symm, const ModuleBase::Matrix3& gmatd, const TCdouble gtransd, const TCdouble& posd_a1, const TCdouble& posd_a2)const @@ -44,10 +42,13 @@ namespace ModuleSymmetry const int& abfs_Lmax = this->abfs_Lmax_; //-------------------------------------------------------------------------------- // setters - void find_irreducible_sector(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, + void find_irred_sector(const Symmetry& symm, const Atom* atoms, + const int nat, const int ntype, + const std::vector& iat2it, const std::vector& iat2ia, const std::vector& Rs, const TC& period, const Lattice& lat, const std::string& output_dir = "") { - this->irs_.find_irreducible_sector(symm, atoms, ucell, Rs, period, lat, output_dir); + this->irs_.find_irred_sector(symm, atoms, nat, ntype, iat2it, iat2ia, + Rs, period, lat, output_dir); } void set_abfs_Lmax(const int l) { this->abfs_Lmax_ = l; } //-------------------------------------------------------------------------------- @@ -56,12 +57,21 @@ namespace ModuleSymmetry /// The top-level calculation interface of this class. calculate the rotation matrix in AO representation: M /// only need once call in each ion step (decided by the configuration) /// @param kstars equal k points to each ibz-kpont, corresponding to a certain symmetry operations. + /// @param symm analyzed symmetry of the cell (gmatrix, gmatrix_anti, kgmatrix, ...) + /// @param atoms all atoms + /// @param latvec Cartesian lattice vectors, as rows of a 3x3 matrix + /// @param lmax maximum angular momentum of the atomic orbitals + /// @param nat total number of atoms + /// @param iat2it map from global atom index iat to its species index it + /// @param iat2ia map from global atom index iat to its index ia within the species /// @param nspin stored as a member so restore_dm()/contruct_2d_rot_mat_ao() do not each /// need to read the global nspin config setting (keeps this LibRI-free class /// free of a module_parameter link dependency; every existing caller already /// has nspin in scope). - void cal_Ms(const K_Vectors& kv, - const UnitCell& ucell, const Parallel_2D& pv, const int nspin); + void cal_Ms(const K_Vectors& kv, const Symmetry& symm, const Atom* atoms, + const ModuleBase::Matrix3& latvec, const int lmax, + const int nat, const std::vector& iat2it, const std::vector& iat2ia, + const Parallel_2D& pv, const int nspin); /// Use calculated M matrix to recover D(k) from D(k_ibz): D(k) = M(R, k)^\dagger D(k_ibz) M(R, k) /// the link "ik_ibz-isym-ik" can be found in kstars: k_bz = gmat[isym](k) @@ -107,7 +117,8 @@ namespace ModuleSymmetry /// 2d-block parallized rotation matrix in AO-representation, denoted as M. /// finally we will use D(k)=M(R, k)^\dagger*D(Rk)*M(R, k) to recover D(k) from D(Rk). - std::vector> contruct_2d_rot_mat_ao(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, + std::vector> contruct_2d_rot_mat_ao(const Symmetry& symm, const Atom* atoms, + const int nat, const std::vector& iat2it, const std::vector& iat2ia, const TCdouble& kvec_d_ibz, int isym, const Parallel_2D& pv, const SpinRotation::Su2& spin_U /*= SpinRotation::Su2{ 1.0, 0.0, 0.0, 1.0 }*/) const; @@ -173,7 +184,7 @@ namespace ModuleSymmetry std::vector spin_U_; /// irreducible sector - Irreducible_Sector irs_; + Irred_Sector irs_; }; } #endif // SYMM_ROTATION_K_H diff --git a/source/source_lcao/module_dftu/dftu_nao_occ.cpp b/source/source_lcao/module_dftu/dftu_nao_occ.cpp index 13a57778f67..47ed453fb20 100644 --- a/source/source_lcao/module_dftu/dftu_nao_occ.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_occ.cpp @@ -362,9 +362,12 @@ void cal_occ_mat_k(const Parallel_Orbitals* pv, if (dftu_spacegroup_symmetry && !dftu_occ_symrot_built) { const std::array period{ kv.nmp[0], kv.nmp[1], kv.nmp[2] }; - dftu_occ_symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell, + dftu_occ_symrot.find_irred_sector(ucell.symm, ucell.atoms, + ucell.nat, ucell.ntype, ucell.iat2it, ucell.iat2ia, ModuleSymmetry::Symmetry_rotation_k::get_bvk_cells(period), period, ucell.lat); - dftu_occ_symrot.cal_Ms(kv, ucell, *pv, nspin); + dftu_occ_symrot.cal_Ms(kv, ucell.symm, ucell.atoms, + ucell.latvec, ucell.lmax, ucell.nat, ucell.iat2it, ucell.iat2ia, + *pv, nspin); dftu_occ_symrot_built = true; } diff --git a/source/source_lcao/module_dftu/dftu_nao_op.cpp b/source/source_lcao/module_dftu/dftu_nao_op.cpp index 6a61dc1cb2d..b670df11e87 100644 --- a/source/source_lcao/module_dftu/dftu_nao_op.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_op.cpp @@ -141,9 +141,13 @@ void hamilt::DFTU_onsite>::contributeHR() { const std::array period{ this->kv_->nmp[0], this->kv_->nmp[1], this->kv_->nmp[2] }; // for return_lattice to calculate Ms - this->symrot_.find_irreducible_sector(this->ucell->symm, this->ucell->atoms, *this->ucell, + this->symrot_.find_irred_sector(this->ucell->symm, this->ucell->atoms, + this->ucell->nat, this->ucell->ntype, this->ucell->iat2it, this->ucell->iat2ia, ModuleSymmetry::Symmetry_rotation_k::get_bvk_cells(period), period, this->ucell->lat); - this->symrot_.cal_Ms(*this->kv_, *this->ucell, *pv, this->nspin); + this->symrot_.cal_Ms(*this->kv_, this->ucell->symm, this->ucell->atoms, + this->ucell->latvec, this->ucell->lmax, + this->ucell->nat, this->ucell->iat2it, this->ucell->iat2ia, + *pv, this->nspin); this->symrot_built_ = true; } const int nspin0 = (this->nspin == 2) ? 2 : 1; diff --git a/source/source_lcao/module_rdmft/rdmft.cpp b/source/source_lcao/module_rdmft/rdmft.cpp index 16c94e978dd..5567140ebc2 100644 --- a/source/source_lcao/module_rdmft/rdmft.cpp +++ b/source/source_lcao/module_rdmft/rdmft.cpp @@ -157,9 +157,13 @@ void RDMFT::init(Parallel_Orbitals& ParaV_in, if (exx_spacegroup_symmetry) { const std::array& period = RI_Util::get_Born_vonKarmen_period(*kv); - this->symrot_exx.find_irreducible_sector(ucell->symm, ucell->atoms, *ucell, + this->symrot_exx.find_irred_sector(ucell->symm, ucell->atoms, + ucell->nat, ucell->ntype, ucell->iat2it, ucell->iat2ia, RI_Util::get_Born_von_Karmen_cells(period), period, ucell->lat, PARAM.globalv.global_out_dir); - this->symrot_exx.cal_Ms(*kv, *ucell, *ParaV, nspin); + this->symrot_exx.cal_Ms(*kv, ucell->symm, ucell->atoms, + ucell->latvec, ucell->lmax, + ucell->nat, ucell->iat2it, ucell->iat2ia, + *ParaV, nspin); } if (this->exx_info_->info_ri.real_number) diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rot_out.cpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rot_out.cpp index 5dbecc1a40f..70971e431b5 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rot_out.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rot_out.cpp @@ -17,7 +17,7 @@ namespace ModuleSymmetry ModuleBase::TITLE("ModuleSymmetry", "print_symrot_info_R"); std::ofstream ofs(PARAM.globalv.global_out_dir + "symrot_R.txt"); // Print the irreducible sector (to be optimized) - ofs << "Number of irreducible sector: " << symrot.get_irreducible_sector().size() << std::endl; + ofs << "Number of irreducible sector: " << symrot.get_irred_sector().size() << std::endl; ofs << "Lmax of AOs: " << lmax_ao << "\n"; ofs << "Lmax of ABFs: " << symrot.abfs_Lmax << "\n"; // print AO rotation matrix T From 8bffde967d5e6872d69e35596baedf55a5e4722e Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 9 Oct 2026 16:00:17 +0800 Subject: [PATCH 08/15] Decouple EXX Symmetry_rotation from UnitCell and drop dead test helpers The EXX/RPA/RDMFT rotation templates only consume the atom index maps, so pass iat2it/iat2ia explicitly instead of const UnitCell&, and remove the unitcell.h include from symm_rotation.h (forward-declare UnitCell for the non-template print_symrot_info_k). TUs that dereference UnitCell now include unitcell.h directly. Delete unreachable helpers: test_HR_rotation (both overloads), test_Cs_rotation, get_Rs_from_adjacent_list, get_Rs_from_BvK. Also fix the three missed .hpp call sites of the Irred_Sector rename (exx_lri_interface.hpp, exx_lri.hpp, rpa_lri.hpp). Verified: make -j 30 builds abacus_max_para with 0 warnings; MODULE_RI_EXX_SYMMETRY_rotation test target also builds. --- source/source_lcao/module_ri/exx_lri.hpp | 7 +- .../module_ri/exx_lri_interface.hpp | 9 +- .../module_exx_symmetry/symm_rot_out.cpp | 1 + .../module_exx_symmetry/symm_rotation.cpp | 51 +---------- .../module_exx_symmetry/symm_rotation.h | 37 +++----- .../module_exx_symmetry/symm_rotation_r.hpp | 83 ++---------------- .../symm_rotation_r_hcontainer.hpp | 85 ++----------------- .../test/symm_rotation_test.cpp | 1 + source/source_lcao/module_ri/rpa_lri.hpp | 8 +- 9 files changed, 48 insertions(+), 234 deletions(-) diff --git a/source/source_lcao/module_ri/exx_lri.hpp b/source/source_lcao/module_ri/exx_lri.hpp index e202851ca9e..335419d9b4e 100644 --- a/source/source_lcao/module_ri/exx_lri.hpp +++ b/source/source_lcao/module_ri/exx_lri.hpp @@ -7,6 +7,7 @@ #define EXX_LRI_HPP #include "exx_lri.h" +#include "source_cell/unitcell.h" #include "ri_2d_comm.h" #include "ri_util.h" #include "exx_lri_detail.h" @@ -789,7 +790,7 @@ void Exx_LRI::cal_exx_elec(const std::vector, std::set>> judge = RI_2D_Comm::get_2D_judge(ucell,pv); if(p_symrot) - { this->exx_lri.set_symmetry(true, p_symrot->get_irreducible_sector()); } + { this->exx_lri.set_symmetry(true, p_symrot->get_irred_sector()); } else { this->exx_lri.set_symmetry(false, {}); } @@ -822,7 +823,7 @@ void Exx_LRI::cal_exx_elec(const std::vectorexx_lri.post_2D.set_tensors_map2(this->exx_lri.Hs); // rotate locally without repeat - Hs_a2D = p_symrot->restore_HR(ucell.symm, ucell.atoms, ucell, 'H', Hs_a2D); + Hs_a2D = p_symrot->restore_HR(ucell.symm, ucell.atoms, ucell.iat2it, 'H', Hs_a2D); // cal energy using full Hs without repeat this->exx_lri.energy = this->exx_lri.post_2D.cal_energy( this->exx_lri.post_2D.saves["Ds_" + suffix], @@ -864,7 +865,7 @@ void Exx_LRI::cal_exx_elec_soc( // pass 2: spinor-coupled rotation of the 4 channels from the irreducible sector to the full BZ std::array>>, 4> Hs_full = - p_symrot->restore_HR_nspin4(ucell.symm, ucell.atoms, ucell, 'H', Hs_irr); + p_symrot->restore_HR_nspin4(ucell.symm, ucell.atoms, ucell.iat2it, 'H', Hs_irr); // pass 3: per-channel energy (full Hs, no repeat), then gather the repeated full Hs for abacus for (int is = 0; is < 4; ++is) diff --git a/source/source_lcao/module_ri/exx_lri_interface.hpp b/source/source_lcao/module_ri/exx_lri_interface.hpp index 20b3c39aee1..c12d80b2653 100644 --- a/source/source_lcao/module_ri/exx_lri_interface.hpp +++ b/source/source_lcao/module_ri/exx_lri_interface.hpp @@ -1,6 +1,7 @@ #ifndef EXX_LRI_INTERFACE_HPP #define EXX_LRI_INTERFACE_HPP #include "exx_lri_interface.h" +#include "source_cell/unitcell.h" #include "source_base/formatter.h" #include "source_base/parallel_common.h" #include "source_estate/elecstate_lcao.h" @@ -90,11 +91,13 @@ void Exx_LRI_Interface::exx_before_all_runners( if (this->exx_spacegroup_symmetry) { const std::array& period = RI_Util::get_Born_vonKarmen_period(kv); - this->symrot_.find_irreducible_sector( - ucell.symm, ucell.atoms, ucell, + this->symrot_.find_irred_sector( + ucell.symm, ucell.atoms, + ucell.nat, ucell.ntype, ucell.iat2it, ucell.iat2ia, RI_Util::get_Born_von_Karmen_cells(period), period, ucell.lat, PARAM.globalv.global_out_dir); this->symrot_.set_abfs_Lmax(Exx_Abfs::Construct_Orbs::get_Lmax(this->exx_ptr->abfs)); - this->symrot_.cal_Ms(kv, ucell, pv, PARAM.inp.nspin); + this->symrot_.cal_Ms(kv, ucell.symm, ucell.atoms, ucell.latvec, ucell.lmax, + ucell.nat, ucell.iat2it, ucell.iat2ia, pv, PARAM.inp.nspin); } } diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rot_out.cpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rot_out.cpp index 70971e431b5..bffd3118cbf 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rot_out.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rot_out.cpp @@ -1,4 +1,5 @@ #include "./symm_rotation.h" +#include "source_cell/unitcell.h" namespace ModuleSymmetry { std::string mat3_fmt(const ModuleBase::Matrix3& m) diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp index b8261ad6fe1..d7bc618ceb6 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp @@ -1,13 +1,6 @@ #include "symm_rotation.h" -#include "source_base/constants.h" -#include "source_io/module_parameter/parameter.h" -#include -#include "source_base/parallel_reduce.h" -#include "source_base/module_external/scalapack_connector.h" -#include "source_base/tool_title.h" -#include "source_base/timer.h" -#include "source_base/mathzone.h" -#include "source_lcao/module_ri/ri_util.h" + +#include namespace ModuleSymmetry { @@ -18,44 +11,4 @@ namespace ModuleSymmetry for (auto& abfs_T : abfs_l_nchi) { this->abfs_Lmax_ = std::max(this->abfs_Lmax_, static_cast(abfs_T.size()) - 1); } } - std::vector Symmetry_rotation::get_Rs_from_adjacent_list(const UnitCell& ucell, - const Grid_Driver& gd, - const Parallel_Orbitals& pv) const - { - // find the union set of Rs for all the atom pairs - std::set Rs_set; - for (int iat1 = 0;iat1 < ucell.nat;++iat1) - { - auto tau1 = ucell.get_tau(iat1); - int it1 = ucell.iat2it[iat1], ia1 = ucell.iat2ia[iat1]; - AdjacentAtomInfo adjs; - gd.Find_atom(ucell, tau1, it1, ia1, &adjs); - for (int ad = 0; ad < adjs.adj_num + 1; ++ad) - { - const int it2 = adjs.ntype[ad]; - const int ia2 = adjs.natom[ad]; - int iat2 = ucell.itia2iat(it2, ia2); - if (pv.get_nrow_atom(iat1) && pv.get_ncol_atom(iat2)) - { - const ModuleBase::Vector3& R_index = adjs.box[ad]; - if (ucell.cal_dtau(iat1, iat2, R_index).norm() * ucell.lat0 - < ucell.atoms[it1].Rcut + ucell.atoms[it2].Rcut) { - Rs_set.insert({ R_index.x, R_index.y, R_index.z }); -} - } - } - } - // set to vector - std::vector Rs(Rs_set.size()); - for (auto& R : Rs_set) { Rs.push_back(R); -} - return Rs; - } - - std::vector Symmetry_rotation::get_Rs_from_BvK(const K_Vectors& kv) const - { - const TC& period = RI_Util::get_Born_vonKarmen_period(kv); - return RI_Util::get_Born_von_Karmen_cells(period); - } - } diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h index 4dfded3ddbf..2adbb5561b2 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.h @@ -5,8 +5,8 @@ #include "source_basis/module_ao/parallel_orbitals.h" #include #include "source_hamilt/module_hcontainer/hcontainer.h" -#include "source_cell/module_neighbor/sltk_grid_driver.h" -#include "source_cell/unitcell.h" + +class UnitCell; namespace ModuleSymmetry { @@ -30,13 +30,17 @@ namespace ModuleSymmetry //-------------------------------------------------------------------------------- /// The main functions to rotate matrices /// Given H(R) in the irreduceble sector, calculate H(R) for all the atompairs and cells. + /// @param iat2it map from global atom index to its species index template // RI::Tensor type std::map, RI::Tensor>> restore_HR( - const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, + const Symmetry& symm, const Atom* atoms, const std::vector& iat2it, const char mode, const std::map, RI::Tensor>>& HR_irreduceble)const; + /// @param iat2it map from global atom index to its species index + /// @param iat2ia map from global atom index to its index within the species template // HContainer type void restore_HR( - const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, + const Symmetry& symm, const Atom* atoms, + const std::vector& iat2it, const std::vector& iat2ia, const char mode, const hamilt::HContainer& HR_irreduceble, hamilt::HContainer& HR_rotated)const; /// (nspin=4) spinor overload: rotate all 4 spin channels of H(R) together. On top of the /// orbital rotation T1^dagger(.)T2 (mode 'H') / T1^T(.)T2^* (mode 'D') applied to every @@ -48,31 +52,17 @@ namespace ModuleSymmetry /// remap of the 4 channels applied after the SU(2) mixing; see symm_rotation_r.hpp. template // RI::Tensor type std::array, RI::Tensor>>, 4> restore_HR_nspin4( - const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, + const Symmetry& symm, const Atom* atoms, const std::vector& iat2it, const char mode, const std::array, RI::Tensor>>, 4>& HR_irreducible_soc)const; //-------------------------------------------------------------------------------- - /// test functions - /// test H(R) rotation: giver a full H(R), pick out H(R) in the irreducible sector, rotate it, and compare with the original full H(R) - template // RI::Tensor type, using col-major implementation - void test_HR_rotation(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, - const std::map, RI::Tensor>>& HR_full); - template // test the rotation of RI coefficients - void test_Cs_rotation(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, - const std::map, RI::Tensor>>& Cs_full)const; - template // HContainer type, using row-major implementation - void test_HR_rotation(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, - const hamilt::HContainer& HR_full); - template // HContainer type + template // HContainer type void print_HR(const std::map, RI::Tensor>>& HR, const std::string name, const double& threshold = 0.0); //-------------------------------------------------------------------------------- private: //-------------------------------------------------------------------------------- - std::vector get_Rs_from_BvK(const K_Vectors& kv)const; - std::vector get_Rs_from_adjacent_list(const UnitCell& ucell, - const Grid_Driver& gd, - const Parallel_Orbitals& pv) const; + //-------------------------------------------------------------------------------- /// The sub functions to rotate matrices @@ -84,8 +74,9 @@ namespace ModuleSymmetry template // pointer type, blas void rotate_atompair_serial(Tdata* TAT, const Tdata* A, const int& nw1, const int& nw2, const int isym, const Atom& a1, const Atom& a2, const char mode)const; - template // HContainer type, pblas - void rotate_atompair_parallel(const TR* Alocal_in, const int isym, const Atom* atoms, const UnitCell& ucell, + template // HContainer type, pblas + void rotate_atompair_parallel(const TR* Alocal_in, const int isym, const Atom* atoms, + const std::vector& iat2it, const std::vector& iat2ia, const Tap& ap_in, const Tap& ap_out, const char mode, const Parallel_Orbitals& pv, TR* Alocal_out, const bool output = false)const; /// rotate a 3-dim C tensor in RI technique diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp index c896f085bb8..86a20e6364f 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r.hpp @@ -77,7 +77,7 @@ namespace ModuleSymmetry template std::map, RI::Tensor>> Symmetry_rotation::restore_HR( - const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, + const Symmetry& symm, const Atom* atoms, const std::vector& iat2it, const char mode, const std::map, RI::Tensor>>& HR_irreducible) const { ModuleBase::TITLE("Symmetry_rotation", "restore_HR"); @@ -102,7 +102,7 @@ namespace ModuleSymmetry const int& ap1 = apR.first.first; const int& ap2 = apR.first.second; const TC& R = apR.second; - HR_full[ap1][{ap2, R}] = rotate_atompair_serial(tmp2.second, isym, atoms[ucell.iat2it[irap1]], atoms[ucell.iat2it[irap2]], mode); + HR_full[ap1][{ap2, R}] = rotate_atompair_serial(tmp2.second, isym, atoms[iat2it[irap1]], atoms[iat2it[irap2]], mode); } } else { std::cout << "Warning: not found: irreducible atom pair =(" << irap1 << "," << irap2 << "), irR=(" << irR[0] << "," << irR[1] << "," << irR[2] << ")\n";} @@ -171,7 +171,7 @@ namespace ModuleSymmetry // when it is applied in either direction. template std::array, RI::Tensor>>, 4> Symmetry_rotation::restore_HR_nspin4( - const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, + const Symmetry& symm, const Atom* atoms, const std::vector& iat2it, const char mode, const std::array, RI::Tensor>>, 4>& HR_irreducible_soc)const { ModuleBase::TITLE("Symmetry_rotation", "restore_HR_nspin4"); @@ -203,8 +203,8 @@ namespace ModuleSymmetry std::cout << "Warning: not found: irreducible atom pair =(" << irap1 << "," << irap2 << "), irR=(" << irR[0] << "," << irR[1] << "," << irR[2] << ")\n"; continue; } - const Atom& a1 = atoms[ucell.iat2it[irap1]]; - const Atom& a2 = atoms[ucell.iat2it[irap2]]; + const Atom& a1 = atoms[iat2it[irap1]]; + const Atom& a2 = atoms[iat2it[irap2]]; // gather the 4 irreducible spin-channel blocks (zero-filled when absent) std::array, 4> Hir; for (int is = 0;is < 4;++is) @@ -405,77 +405,4 @@ namespace ModuleSymmetry } } - template - void Symmetry_rotation::test_HR_rotation(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, - const std::map, RI::Tensor>>& HR_full) - { - ModuleBase::TITLE("Symmetry_rotation", "test_HR_rotation"); - - // 1. pick out H(R) in the irreducible sector from full H(R) - std::map, RI::Tensor>> HR_irreducible; - for (auto& irap_Rs : this->irs_.irreducible_sector_) - { - const Tap& irap = irap_Rs.first; - for (auto& irR : irap_Rs.second) - { - const std::pair a2_irR = { irap.second, irR }; - HR_irreducible[irap.first][a2_irR] = (HR_full.at(irap.first).count(a2_irR) != 0) ? - HR_full.at(irap.first).at(a2_irR) - : RI::Tensor(HR_full.at(irap.first).begin()->second.shape); - } - } - // 2. rotate - std::map, RI::Tensor>> HR_rotated = restore_HR(symm, atoms, ucell, mode, HR_irreducible); - // 3. compare - for (auto& HR_ia1 : HR_rotated) - { - int iat1 = HR_ia1.first; - for (auto& HR_ia12R : HR_ia1.second) - { - int iat2 = HR_ia12R.first.first; - TC R = HR_ia12R.first.second; - const RI::Tensor& HR_rot = HR_ia12R.second; - if (HR_full.at(iat1).count({ iat2, R }) == 0)// rot back but not found - { - std::cout << "R_rot not found in atom pair (" << iat1 << ", " << iat2 << "): R=(" << R[0] << "," << R[1] << "," << R[2] << "):\n"; - continue; - } - const RI::Tensor& HR_ref = HR_full.at(iat1).at({ iat2, R }); - assert(HR_rot.shape[0] == HR_ref.shape[0]); - assert(HR_rot.shape[1] == HR_ref.shape[1]); - // output - std::cout << "atom pair (" << iat1 << ", " << iat2 << "), R=(" << R[0] << "," << R[1] << "," << R[2] << "):\n"; - print_tensor(HR_rot, std::string("R_rot").insert(0, 1, mode)); - print_tensor(HR_ref, std::string("R_ref").insert(0, 1, mode)); - } - } - } - - template - void Symmetry_rotation::test_Cs_rotation(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, - const std::map, RI::Tensor>>& Cs_full)const - { - for (auto& sector_pair : this->irs_.full_map_to_irreducible_sector_) - { - const TapR& apR = sector_pair.first; - const int& isym = sector_pair.second.first; - const TapR& irapR = sector_pair.second.second; - // if (apR.first != irapR.first || apR.second != irapR.second) - if (apR.first != irapR.first) - { - std::cout << "irapR=(" << irapR.first.first << "," << irapR.first.second << "), (" << irapR.second[0] << "," << irapR.second[1] << "," << irapR.second[2] << "):\n"; - std::cout << "apR=(" << apR.first.first << "," << apR.first.second << "), (" << apR.second[0] << "," << apR.second[1] << "," << apR.second[2] << "):\n"; - const RI::Tensor& Cs_ir = Cs_full.at(irapR.first.first).at({ irapR.first.second,irapR.second }); - const RI::Tensor& Cs_ref = Cs_full.at(apR.first.first).at({ apR.first.second,apR.second }); - const RI::Tensor& Cs_rot = this->rotate_singleC_serial(Cs_ir, isym, - atoms[ucell.iat2it[irapR.first.first]], atoms[ucell.iat2it[irapR.first.second]], irapR.first.first); - print_tensor3(Cs_rot, "Cs_rot"); - print_tensor3(Cs_ref, "Cs_ref"); - print_tensor3(Cs_ir, "Cs_irreducible"); - exit(0); - } - } - - } - } \ No newline at end of file diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r_hcontainer.hpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r_hcontainer.hpp index 61d5e344490..ad0f7274952 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r_hcontainer.hpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation_r_hcontainer.hpp @@ -36,7 +36,8 @@ namespace ModuleSymmetry template // HContainer type void Symmetry_rotation::restore_HR( - const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, const char mode, + const Symmetry& symm, const Atom* atoms, + const std::vector& iat2it, const std::vector& iat2ia, const char mode, const hamilt::HContainer& HR_irreduceble, hamilt::HContainer& HR_rotated)const { @@ -59,7 +60,7 @@ namespace ModuleSymmetry const int R_hc = ap_hc.find_R(R[0], R[1], R[2]); if (R_hc < 0) continue; TR* ijR_ptr = ap_hc.get_pointer(R_hc); - rotate_atompair_parallel(irijR_ptr, isym, atoms, ucell, irap, ap, mode, *irap_hc.get_paraV(), ijR_ptr); + rotate_atompair_parallel(irijR_ptr, isym, atoms, iat2it, iat2ia, irap, ap, mode, *irap_hc.get_paraV(), ijR_ptr); } ModuleBase::timer::end("Symmetry_rotation", "restore_HR"); } @@ -77,13 +78,14 @@ namespace ModuleSymmetry obj[starti + i + (startj + j) * nr] = block(j, i); }; template - void Symmetry_rotation::rotate_atompair_parallel(const TR* Alocal_in, const int isym, const Atom* atoms, const UnitCell& ucell, + void Symmetry_rotation::rotate_atompair_parallel(const TR* Alocal_in, const int isym, const Atom* atoms, + const std::vector& iat2it, const std::vector& iat2ia, const Tap& ap_in, const Tap& ap_out, const char mode, const Parallel_Orbitals& pv, TR* Alocal_out, const bool output)const { // all the matrices are row-major (col-contiguous) int iat1 = ap_in.first, iat2 = ap_in.second; - int it1 = ucell.iat2it[iat1], it2 = ucell.iat2it[iat2]; - int ia1 = ucell.iat2ia[iat1], ia2 = ucell.iat2ia[iat2]; + int it1 = iat2it[iat1], it2 = iat2it[iat2]; + int ia1 = iat2ia[iat1], ia2 = iat2ia[iat2]; // contruct T matrix std::vector T1, T2; auto set_rotation_matrix = [&](const int& it, const int& ia, std::vector& T)->int @@ -149,8 +151,8 @@ namespace ModuleSymmetry // copy back to Alocal_out iat1 = ap_out.first, iat2 = ap_out.second; - it1 = ucell.iat2it[iat1], it2 = ucell.iat2it[iat2]; - ia1 = ucell.iat2ia[iat1], ia2 = ucell.iat2ia[iat2]; + it1 = iat2it[iat1], it2 = iat2it[iat2]; + ia1 = iat2ia[iat1], ia2 = iat2ia[iat2]; abr1 = pv.atom_begin_row[iat1], abc2 = pv.atom_begin_col[iat2]; iw1start = atoms[it1].stapos_wf + ia1 * atoms[it1].nw; iw2start = atoms[it2].stapos_wf + ia2 * atoms[it2].nw; @@ -162,73 +164,4 @@ namespace ModuleSymmetry = TAT[(pv.local2global_row(i + abr1) - iw1start) * atoms[it2].nw + (pv.local2global_col(j + abc2) - iw2start)]; } } - - template - void Symmetry_rotation::test_HR_rotation(const Symmetry& symm, const Atom* atoms, const UnitCell& ucell, - const char mode, const hamilt::HContainer& HR_full) - { - ModuleBase::TITLE("Symmetry_rotation", "test_HR_rotation"); - auto get_irreducible_ijR_info = [&HR_full, this]() -> std::vector - { - std::vector irreducible_ijR_info; - irreducible_ijR_info.push_back(this->irs_.irreducible_sector_.size()); - for (auto& irap_irR : this->irs_.irreducible_sector_) - { - const int iat1 = irap_irR.first.first, iat2 = irap_irR.first.second; - irreducible_ijR_info.insert(irreducible_ijR_info.end(), { iat1, iat2, 0 }); - int nR = 0; - for (auto& R : irap_irR.second) - if (HR_full.get_atom_pair(iat1, iat2).find_R(R[0], R[1], R[2]) >= 0) - { - irreducible_ijR_info.insert(irreducible_ijR_info.end(), R.begin(), R.end()); - ++nR; - } - irreducible_ijR_info[irreducible_ijR_info.size() - 3 * nR - 1] = nR; - } - return irreducible_ijR_info; - }; - - const Parallel_Orbitals* pv = HR_full.get_atom_pair(0, 0).get_paraV(); - // 1. pick out H(R) in the irreducible sector from full H(R) - const std::vector& irreducible_ijR_info = get_irreducible_ijR_info(); - hamilt::HContainer HR_irreducible(pv, nullptr, &irreducible_ijR_info); - HR_irreducible.set_zero(); - for (auto& irap_irR : this->irs_.irreducible_sector_) - { - const int iat1 = irap_irR.first.first, iat2 = irap_irR.first.second; - const hamilt::AtomPair& irap = HR_irreducible.get_atom_pair(iat1, iat2); - const hamilt::AtomPair& ap_full = HR_full.get_atom_pair(iat1, iat2); - for (auto& R : irap_irR.second) - if (irap.find_R(R[0], R[1], R[2]) >= 0) - { // out-of-range R can be added to irreducible sector to be calculated, but not in this test for lack of reference - TR* irptr = irap.get_HR_values(R[0], R[1], R[2]).get_pointer(); - TR* ptr = ap_full.get_HR_values(R[0], R[1], R[2]).get_pointer(); - std::copy(ptr, ptr + pv->get_nrow_atom(iat1) * pv->get_ncol_atom(iat2), irptr); - } - } - - //2. rotate - hamilt::HContainer HR_rotated(HR_full); - HR_rotated.set_zero(); - this->restore_HR(symm, atoms, ucell, mode, HR_irreducible, HR_rotated); - //3. compare - for (int iat1 = 0;iat1 < ucell.nat;++iat1) - { - for (int iat2 = 0;iat2 < ucell.nat;++iat2) - { - const hamilt::AtomPair& ap_full = HR_full.get_atom_pair(iat1, iat2); - const hamilt::AtomPair& ap_rotated = HR_rotated.get_atom_pair(iat1, iat2); - assert(ap_full.get_R_size() == ap_rotated.get_R_size()); - for (int irR = 0;irR < ap_full.get_R_size();++irR) - { - const TR* full_ptr = ap_full.get_pointer(irR); - int* R = ap_full.get_R_index(irR); - const TR* rotated_ptr = ap_rotated.get_HR_values(R[0], R[1], R[2]).get_pointer(); - GlobalV::ofs_running << "atom pair: (" << iat1 << ", " << iat2 << "), R: (" << R[0] << " " << R[1] << " " << R[2] << ")\n"; - print_atompair_local(*pv, iat1, iat2, rotated_ptr, std::string("R_rot").insert(0, 1, mode)); - print_atompair_local(*pv, iat1, iat2, full_ptr, std::string("R_ref").insert(0, 1, mode)); - } - } - } - } } \ No newline at end of file diff --git a/source/source_lcao/module_ri/module_exx_symmetry/test/symm_rotation_test.cpp b/source/source_lcao/module_ri/module_exx_symmetry/test/symm_rotation_test.cpp index 483df5d894a..ed032ec60d1 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/test/symm_rotation_test.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/test/symm_rotation_test.cpp @@ -1,5 +1,6 @@ #include "mpi.h" #include "../symm_rotation.h" +#include "source_cell/unitcell.h" #include "gtest/gtest.h" #define DOUBLETHRESHOLD 1e-8 diff --git a/source/source_lcao/module_ri/rpa_lri.hpp b/source/source_lcao/module_ri/rpa_lri.hpp index 8c4e3485465..e30641d8f44 100644 --- a/source/source_lcao/module_ri/rpa_lri.hpp +++ b/source/source_lcao/module_ri/rpa_lri.hpp @@ -16,6 +16,7 @@ #include #include #include "source_lcao/module_ri/module_exx_symmetry/symm_rotation.h" +#include "source_cell/unitcell.h" #include "rpa_lri.h" #include "exx_lri.h" @@ -167,13 +168,16 @@ void RPA_LRI::cal_postSCF_exx(const module_dm::DensityMatrix { const std::array period = RI_Util::get_Born_vonKarmen_period(kv); const auto& Rs = RI_Util::get_Born_von_Karmen_cells(period); - symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell, Rs, period, ucell.lat, PARAM.globalv.global_out_dir); + symrot.find_irred_sector(ucell.symm, ucell.atoms, + ucell.nat, ucell.ntype, ucell.iat2it, ucell.iat2ia, + Rs, period, ucell.lat, PARAM.globalv.global_out_dir); // set Lmax of the rotation matrices to max(l_ao, l_abf), to support rotation under ABF // NOTE: Using Exx_Abfs::Construct_Orbs::get_Lmax() to compute Lmax from the actual ABFs // instead of relying on exx_cut_coulomb->abfs_Lmax() (not yet initialized) or // this->info.abfs_Lmax (defaults to 0). This ensures correct Lmax for symmetry rotation. symrot.set_abfs_Lmax(Exx_Abfs::Construct_Orbs::get_Lmax(abfs_for_lmax)); - symrot.cal_Ms(kv, ucell, parav, PARAM.inp.nspin); + symrot.cal_Ms(kv, ucell.symm, ucell.atoms, ucell.latvec, ucell.lmax, + ucell.nat, ucell.iat2it, ucell.iat2ia, parav, PARAM.inp.nspin); // output Ts (symrot_R.txt) and Ms (symrot_k.txt) ModuleSymmetry::print_symrot_info_R(symrot, ucell.symm, ucell.lmax, Rs); ModuleSymmetry::print_symrot_info_k(symrot, kv, ucell); From e37a07fb55e5f743e06ff2a7482aaadd3f1a5a21 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 9 Oct 2026 16:17:47 +0800 Subject: [PATCH 09/15] Point integration tests to local build path Set the ABACUS executable in Autotest.sh and general_info to the local build_max_para_test/abacus_max_para binary for local test runs. --- tests/integrate/Autotest.sh | 2 +- tests/integrate/general_info | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/integrate/Autotest.sh b/tests/integrate/Autotest.sh index 48b6d2308f7..626fe461403 100755 --- a/tests/integrate/Autotest.sh +++ b/tests/integrate/Autotest.sh @@ -1,7 +1,7 @@ #!/bin/bash # ABACUS executable path -abacus=abacus +abacus=/home/510Group/6_abacus_mc/abacus-mc-2/build_max_para_test/abacus_max_para # number of MPI processes np=4 nt=$OMP_NUM_THREADS # number of OpenMP threads, default is $OMP_NUM_THREADS diff --git a/tests/integrate/general_info b/tests/integrate/general_info index c0d397efdc5..e31797aca11 100644 --- a/tests/integrate/general_info +++ b/tests/integrate/general_info @@ -1,3 +1,3 @@ -EXEC abacus +EXEC /home/510Group/6_abacus_mc/abacus-mc-2/build_max_para_test/abacus_max_para CHECKACCURACY 2 NUMBEROFPROCESS 4 From c346731c24bd339f03c910777a0b0ad1644763ba Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 9 Oct 2026 16:20:38 +0800 Subject: [PATCH 10/15] Fix Intel Makefile build after irred_sec rename Makefile.Objects still listed irreducible_sector.o, which no longer exists after renaming irreducible_sector.{h,cpp} to irred_sec.{h,cpp}; update it to irred_sec.o so the mpiicpx Makefile build links again. --- source/Makefile.Objects | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 62df335c194..7c7c93b79e6 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -627,7 +627,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\ From 91ff5982daa279b746f8ff3786e9bae20239973c Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 9 Oct 2026 16:58:16 +0800 Subject: [PATCH 11/15] Fix iat2it/iat2ia allocation in relax_sync unit test The upstream merge brought in a test that assigned new int[nat] to ucell.iat2it/iat2ia, but these are std::vector members on this branch. Use resize() instead, matching the other fixtures in this file and removing the leaked allocation. Verified: MODULE_RELAX_relax_sync_test target builds with 0 errors. --- source/source_relax/unittests/test_relax_sync.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/source_relax/unittests/test_relax_sync.cpp b/source/source_relax/unittests/test_relax_sync.cpp index a1cb0dda25d..0405bcf59bc 100644 --- a/source/source_relax/unittests/test_relax_sync.cpp +++ b/source/source_relax/unittests/test_relax_sync.cpp @@ -725,8 +725,8 @@ TEST(RelaxSyncMaxStep, GeometryMovedAndNotConverged) ucell.atoms[0].label = "Si"; ucell.omega = 1.0; ucell.lat0 = 1.0; - ucell.iat2it = new int[nat]; - ucell.iat2ia = new int[nat]; + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); ucell.iat2it[0] = 0; ucell.iat2ia[0] = 0; ucell.atoms[0].mbl.resize(nat); From b6be72d9319e445b90e1595c2552ea3dadfa54c9 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 9 Oct 2026 20:40:00 +0800 Subject: [PATCH 12/15] Revert general_info EXEC to 'abacus' The local absolute build path broke CUDA CI jobs that do not have this directory; restore the upstream default of resolving abacus from PATH. --- tests/integrate/general_info | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/integrate/general_info b/tests/integrate/general_info index e31797aca11..c0d397efdc5 100644 --- a/tests/integrate/general_info +++ b/tests/integrate/general_info @@ -1,3 +1,3 @@ -EXEC /home/510Group/6_abacus_mc/abacus-mc-2/build_max_para_test/abacus_max_para +EXEC abacus CHECKACCURACY 2 NUMBEROFPROCESS 4 From cac9d064172b8593ec0b118d7d56a28deb835c2f Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Fri, 9 Oct 2026 21:37:58 +0800 Subject: [PATCH 13/15] Revert Autotest.sh executable to 'abacus' Complete the partial revert in b6be72d931: the hardcoded local path build_max_para_test/abacus_max_para was still set in Autotest.sh, so CUDA CI jobs (which lack that local CPU build) failed to launch every integration test case. Restore the upstream default and resolve the executable via PATH, ABACUS_EXE, or the -a option for local runs. --- tests/integrate/Autotest.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/integrate/Autotest.sh b/tests/integrate/Autotest.sh index 626fe461403..48b6d2308f7 100755 --- a/tests/integrate/Autotest.sh +++ b/tests/integrate/Autotest.sh @@ -1,7 +1,7 @@ #!/bin/bash # ABACUS executable path -abacus=/home/510Group/6_abacus_mc/abacus-mc-2/build_max_para_test/abacus_max_para +abacus=abacus # number of MPI processes np=4 nt=$OMP_NUM_THREADS # number of OpenMP threads, default is $OMP_NUM_THREADS From 79b64c04bda4ed1727552229f0894e2950ec6367 Mon Sep 17 00:00:00 2001 From: Critsium Date: Sat, 10 Oct 2026 16:25:31 +0800 Subject: [PATCH 14/15] Fix stale Statistics/Irreducible_Sector uses after merging #8099 The previous merge of upstream/develop brought in #8099, whose code still used the pre-refactor API and broke every CI build: - symm_rotation_k.cpp: reset_symmetry() assigned Irreducible_Sector(), renamed to Irred_Sector on this branch. - symm_test_analysis.cpp and test_symm_rotation.cpp: replace ucell.st, find_irreducible_sector() and cal_Ms(kv, cell, ...) with the new explicit nat/ntype/iat2it/iat2ia/itia2iat signatures, and assign the index maps as std::vector instead of new int[]. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../module_symmetry/symm_rotation_k.cpp | 2 +- .../test/symm_test_analysis.cpp | 10 +++--- .../test/test_symm_rotation.cpp | 36 ++++++++++++------- 3 files changed, 30 insertions(+), 18 deletions(-) diff --git a/source/source_cell/module_symmetry/symm_rotation_k.cpp b/source/source_cell/module_symmetry/symm_rotation_k.cpp index eae54d83869..ac145a67e16 100644 --- a/source/source_cell/module_symmetry/symm_rotation_k.cpp +++ b/source/source_cell/module_symmetry/symm_rotation_k.cpp @@ -12,7 +12,7 @@ namespace ModuleSymmetry { void Symmetry_rotation_k::reset_symmetry() { - this->irs_ = Irreducible_Sector(); + this->irs_ = Irred_Sector(); this->Ms_.clear(); this->little_groups_.clear(); this->spin_U_.clear(); diff --git a/source/source_cell/module_symmetry/test/symm_test_analysis.cpp b/source/source_cell/module_symmetry/test/symm_test_analysis.cpp index ad02f425073..389342700bd 100644 --- a/source/source_cell/module_symmetry/test/symm_test_analysis.cpp +++ b/source/source_cell/module_symmetry/test/symm_test_analysis.cpp @@ -55,8 +55,9 @@ TEST_F(SymmetryTest, AtomMapRefreshWithUnchangedOperationCount) 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); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, calculation, representation); const int original_count = symm.nrotk; EXPECT_EQ(original_count, 48); @@ -70,8 +71,9 @@ TEST_F(SymmetryTest, AtomMapRefreshWithUnchangedOperationCount) { 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); + symm.analy_sys(ucell.lat, ucell.atoms, ucell.nat, ucell.ntype, + ucell.iat2it, ucell.iat2ia, ucell.itia2iat, + ofs_running, 1e-6, 1, calculation, representation); EXPECT_EQ(symm.nrotk, original_count); for (int operation = 0; operation < symm.nrotk; ++operation) { diff --git a/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp b/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp index cee79ba3fb6..131d9d21150 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/test/test_symm_rotation.cpp @@ -174,8 +174,8 @@ TEST(SymmetryDensityRestoration, RebuildAfterCellSymmetryAnalysis) cell.a1 = ModuleBase::Vector3(1, 0, 0); cell.a2 = ModuleBase::Vector3(0, 1, 0); cell.a3 = ModuleBase::Vector3(0, 0, 1); - cell.st.iat2it = new int[3]{0, 0, 1}; - cell.st.iat2ia = new int[3]{0, 1, 0}; + cell.iat2it = {0, 0, 1}; + cell.iat2ia = {0, 1, 0}; atoms[0].label = "A"; atoms[0].na = 2; atoms[0].nw = 1; @@ -193,7 +193,9 @@ TEST(SymmetryDensityRestoration, RebuildAfterCellSymmetryAnalysis) std::ofstream log; const int representation[2] = {0, 0}; const std::string calculation = "cell-relax"; - cell.symm.analy_sys(cell.lat, cell.st, cell.atoms, log, 1e-6, 1, calculation, representation); + cell.symm.analy_sys(cell.lat, cell.atoms, cell.nat, cell.ntype, + cell.iat2it, cell.iat2ia, cell.itia2iat, + log, 1e-6, 1, calculation, representation); const int old_operations = cell.symm.nrotk; K_Vectors kv; kv.set_nks(1); @@ -213,23 +215,31 @@ TEST(SymmetryDensityRestoration, RebuildAfterCellSymmetryAnalysis) const ModuleSymmetry::TC period = {2, 2, 2}; const std::vector cells = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; CellSymmetryRotation reused; - reused.find_irreducible_sector(cell.symm, cell.atoms, cell.st, cells, period, cell.lat); - reused.cal_Ms(kv, cell, pv, 1); + reused.find_irred_sector(cell.symm, cell.atoms, cell.nat, cell.ntype, + cell.iat2it, cell.iat2ia, cells, period, cell.lat); + reused.cal_Ms(kv, cell.symm, cell.atoms, cell.latvec, cell.lmax, + cell.nat, cell.iat2it, cell.iat2ia, pv, 1); const auto old_rotations = reused.rotations(); atoms[0].taud = {{0.9, 0.2, 0.3}, {0.2, 0.3, 0.2}}; atoms[1].taud = {{0.05, 0.25, 0.25}}; atoms[0].tau = atoms[0].taud; atoms[1].tau = atoms[1].taud; - cell.symm.analy_sys(cell.lat, cell.st, cell.atoms, log, 1e-6, 1, calculation, representation); + cell.symm.analy_sys(cell.lat, cell.atoms, cell.nat, cell.ntype, + cell.iat2it, cell.iat2ia, cell.itia2iat, + log, 1e-6, 1, calculation, representation); ASSERT_EQ(cell.symm.nrotk, old_operations); reused.reset_symmetry(); - reused.find_irreducible_sector(cell.symm, cell.atoms, cell.st, cells, period, cell.lat); - reused.cal_Ms(kv, cell, pv, 1); + reused.find_irred_sector(cell.symm, cell.atoms, cell.nat, cell.ntype, + cell.iat2it, cell.iat2ia, cells, period, cell.lat); + reused.cal_Ms(kv, cell.symm, cell.atoms, cell.latvec, cell.lmax, + cell.nat, cell.iat2it, cell.iat2ia, pv, 1); CellSymmetryRotation fresh; - fresh.find_irreducible_sector(cell.symm, cell.atoms, cell.st, cells, period, cell.lat); - fresh.cal_Ms(kv, cell, pv, 1); - EXPECT_EQ(reused.get_irreducible_sector(), fresh.get_irreducible_sector()); + fresh.find_irred_sector(cell.symm, cell.atoms, cell.nat, cell.ntype, + cell.iat2it, cell.iat2ia, cells, period, cell.lat); + fresh.cal_Ms(kv, cell.symm, cell.atoms, cell.latvec, cell.lmax, + cell.nat, cell.iat2it, cell.iat2ia, pv, 1); + EXPECT_EQ(reused.get_irred_sector(), fresh.get_irred_sector()); EXPECT_EQ(reused.rotations(), fresh.rotations()); EXPECT_NE(old_rotations, fresh.rotations()); for (int atom = 0; atom < cell.nat; ++atom) @@ -249,7 +259,7 @@ TEST(SymmetryDensityRestoration, ReturnLatticePreservesBoundaryRepresentatives) { ModuleSymmetry::Symmetry symmetry; symmetry.epsilon = 3.2e-5; - ModuleSymmetry::Irreducible_Sector sector; + ModuleSymmetry::Irred_Sector sector; const ModuleBase::Matrix3 reflection(-1, 0, 0, 0, 1, 0, 0, 0, 1); const ModuleBase::Vector3 translation(0.0, 0.0, 0.0); const ModuleBase::Vector3 source(0.99999, 0.25, 0.5); @@ -274,7 +284,7 @@ TEST(SymmetryDensityRestoration, ReturnLatticePreservesTranslationRepresentative { ModuleSymmetry::Symmetry symmetry; symmetry.epsilon = 3.2e-5; - ModuleSymmetry::Irreducible_Sector sector; + ModuleSymmetry::Irred_Sector sector; const ModuleBase::Matrix3 identity(1, 0, 0, 0, 1, 0, 0, 0, 1); const ModuleBase::Vector3 source(0.25, 0.5, 0.75); const ModuleBase::Vector3 translation(0.99999, 0.0, 0.0); From 97a935d4fb2097cbfed8de6d562b0f2980554f56 Mon Sep 17 00:00:00 2001 From: Critsium Date: Sat, 10 Oct 2026 18:03:55 +0800 Subject: [PATCH 15/15] Fix test_hsr_writer after merging #8068: iat2it is now std::vector #8068 added a test that assigned `new int[1]` to ucell.iat2it, which no longer compiles now that UnitCell owns iat2it as std::vector. Use resize() as the other test in this file already does. Co-Authored-By: Claude Opus 5.5 (1M context) --- source/source_io/module_hs/unittests/test_hsr_writer.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/source_io/module_hs/unittests/test_hsr_writer.cpp b/source/source_io/module_hs/unittests/test_hsr_writer.cpp index df5ff7b68b1..a5f86e7c74e 100644 --- a/source/source_io/module_hs/unittests/test_hsr_writer.cpp +++ b/source/source_io/module_hs/unittests/test_hsr_writer.cpp @@ -506,8 +506,8 @@ TEST(HsrWriterIo, WriteHsrTextCsrCarriesPerSpinFermiFromEfermi) init_unitcell(ucell); // The __MPI path in write_hsr calls set_atomic_trace and gatherParallels, // which require a real atom-to-orbital map; nullptr with nat=0 throws on - // that path. iat2it is released by UnitCell's Statistics member. - ucell.iat2it = new int[1]; + // that path. + ucell.iat2it.resize(1); ucell.iat2it[0] = 0; ucell.set_iat2iwt(1); const int* iat2iwt_ptr = ucell.get_iat2iwt();