From 72353ee125efbd24928247f76c95fe67345a7311 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Wed, 7 Oct 2026 10:18:33 -0400 Subject: [PATCH 1/7] fix(exx): refresh hybrid state and Coulomb tables across cell steps Refresh the initial hybrid Hamiltonian and EXX energy using the current geometry on later ionic steps. Retain Eexx alongside Hexxs until both are replaced, instead of resetting only the energy after convergence. Reinitialize convolved orbitals, integral tables and associated RI caches when the cell-dependent Coulomb settings change, including Spencer cutoff. Reuse the existing auxiliary basis and leave fixed-kernel settings unchanged. Add a two-atom diamond DDH cell-relax regression that compares the evaluated second-step structure with a fresh DDH single point using the same FFT grid. Both separate and simultaneous EXX loops pass a predefined 1e-5 eV tolerance: separate_loop=1: 3.7717291e-6 eV; separate_loop=0: 3.7752644e-6 eV. Both ionic SCFs and independent single points converge, and second-step initial EXX energies remain nonzero. Validation: cmake --build build -j6; Python regression for both EXX loops; py_compile; git diff --cached --check; staged governance check (exit 0). No INPUT names, public APIs, or default arguments change, so no user-facing documentation update is needed. Governance's default-parameter warning is an existing INPUT comparison, not a default argument; global migration budget is neutral. Existing diagnostic source edits are excluded. fix(exx): rebuild symmetry after cell-relax analysis Refresh LCAO EXX symmetry immediately after cell-relax reanalyzes the unit-cell symmetry. Clear lattice returns, sectors, Bloch rotations and angular-rotation caches before rebuilding, including invalidation of the RI tensor cache. Keep the extra refresh limited to cell-relax. Add a regression with unchanged point-group size and atoms crossing a cell boundary. The refreshed sectors, lattice returns and M matrices match a fresh object. Validation: - cmake --build build --target abacus_std_para -j2: passed. - cmake --build build_test --target MODULE_RI_EXX_SYMMETRY_rotation -j2: passed; OMP_NUM_THREADS=1 ctest --test-dir build_test -V -R '^MODULE_RI_EXX_SYMMETRY_rotation$': all 13 tests passed. - python3 tests/08_RI/test_exx_lri_interface.py with the built executable, SG15/Orb-v2 data and one MPI rank: both loops passed the 1e-5 eV tolerance (3.7717291e-6 and 3.7752644e-6 eV). - git diff --cached --check and staged agent governance check: passed; the documentation warning needs no parameter docs update because no INPUT option, default or intended user-facing semantics change. No new GlobalV/GlobalC/PARAM dependencies. Pre-existing unrelated edits are excluded. refactor(exx): separate Coulomb refresh from basis initialization Retain original orbitals, auxiliary basis and the shared Gaunt table when cell-dependent Coulomb parameters change. Rebuild only changed kernels' convolved orbitals, integral interpolations and point-value caches. Destroy old kernel objects first because Matrix_Orbs initialization inserts entries without clearing previous tables. Preserve unchanged kernel objects. The regression checks both EXX loop modes against independent singlepoints and requires one full init and two Coulomb refreshes across two ionic steps. No INPUT parameters or intended numerical behavior change; no documentation update is required. This extraction adds no global dependencies. Verification: - cmake --build build --target abacus_std_para -j2: passed. - python3 -m py_compile tests/08_RI/test_exx_lri_interface.py: passed. - OMP_NUM_THREADS=1 python3 tests/08_RI/test_exx_lri_interface.py --binary build/abacus_std_para --pseudo-dir /home/fortneu/work/pp_orbs/SG15_ONCV_v1.0_upf --orbital-dir /home/fortneu/work/pp_orbs/Orb-v2 --work-dir .diagnostics/ddh-cell-audit/coulomb-refactor-loop1-20261008 --mpi-ranks 1 --separate-loop 1: passed, delta 3.7717290979344398e-6 eV. - Same runtime command with work directory coulomb-refactor-loop0-20261008 and --separate-loop 0: passed, delta 3.77526436068365e-6 eV. - Both runtime modes: basis_initializations=1, coulomb_refreshes=2. - Executable version: ABACUS v3.11.0-beta10. - python3 tools/03_code_analysis/agent_governance_check.py --staged: no blockers; documentation warning explained above. - git diff --cached --check: passed. refactor(exx): refresh symmetry without adding ESolver hooks Remove the new ESolver_FP virtual hook and ESolver_KS_LCAO override. Refresh EXX symmetry in the existing LCAO before_scf immediately after its base call, which has completed symmetry analysis and k-point updates. Keep refresh limited to cell-relax with an updated cell and symmetry enabled, before Hamiltonian construction and EXX ionic-step preparation. Verification: git diff --cached --check passed; staged agent governance check passed without blockers. Test/documentation review warnings apply: per explicit user instruction, the build was stopped and regression was not run for this follow-up; the user will validate actual cases. Existing symmetry and EXX regression coverage is retained. No INPUT behavior or user-facing interface changes require documentation updates. refactor(exx): refresh symmetry inside ionic-step preparation Extract symmetry initialization into Exx_LRI_Interface::refresh_symmetry. Use it from exx_before_all_runners and at the beginning of exx_beforescf on subsequent cell-relax ionic steps with runtime symmetry enabled. Pass parallel orbital layout and the refresh condition explicitly. Update both ground-state and TDDFT callers and remove the ESolver refresh block and Exx_NAO refresh forwarding function. Runtime symm_flag is used because INPUT conversion can disable symmetry after parsing the symmetry setting. Verification: all call sites inspected; git diff --cached --check and staged agent governance check passed without blockers. No build or regression run, per the user's instruction to commit directly and validate actual cases later. Existing tests are retained. No INPUT options/defaults change, so parameter documentation updates are unnecessary. The existing symmetry flag read moves from ESolver to EXX setup with no net global budget increase. --- .../module_symmetry/symm_rotation_k.cpp | 11 +++ .../module_symmetry/symm_rotation_k.h | 3 + source/source_esolver/esolver_ks_lcao.cpp | 2 +- .../source_esolver/esolver_ks_lcao_tddft.cpp | 2 +- source/source_lcao/module_ri/exx_lri.h | 14 ++- source/source_lcao/module_ri/exx_lri.hpp | 61 +++++++++--- .../source_lcao/module_ri/exx_lri_interface.h | 10 +- .../module_ri/exx_lri_interface.hpp | 38 +++++++- .../module_exx_symmetry/test/CMakeLists.txt | 3 +- .../test/test_symm_rotation.cpp | 97 +++++++++++++++++++ source/source_lcao/setup_exx.cpp | 8 +- source/source_lcao/setup_exx.h | 2 + 12 files changed, 221 insertions(+), 30 deletions(-) diff --git a/source/source_cell/module_symmetry/symm_rotation_k.cpp b/source/source_cell/module_symmetry/symm_rotation_k.cpp index f645419b05f..71b145a2147 100644 --- a/source/source_cell/module_symmetry/symm_rotation_k.cpp +++ b/source/source_cell/module_symmetry/symm_rotation_k.cpp @@ -10,6 +10,17 @@ namespace ModuleSymmetry { + void Symmetry_rotation_k::reset_symmetry() + { + this->irs_ = Irreducible_Sector(); + this->Ms_.clear(); + this->little_groups_.clear(); + this->spin_U_.clear(); + this->rotmat_Slm_.clear(); + // Invalidate the derived RI::Tensor rotation cache as well. + ++this->rotmat_Slm_version_; + } + std::vector Symmetry_rotation_k::get_bvk_cells(const TC& period) { std::vector cells; diff --git a/source/source_cell/module_symmetry/symm_rotation_k.h b/source/source_cell/module_symmetry/symm_rotation_k.h index 48ee839f4ff..40efdc18ba9 100644 --- a/source/source_cell/module_symmetry/symm_rotation_k.h +++ b/source/source_cell/module_symmetry/symm_rotation_k.h @@ -24,6 +24,9 @@ namespace ModuleSymmetry Symmetry_rotation_k() {}; virtual ~Symmetry_rotation_k() {}; + /// Discard geometry-dependent caches before rebuilding symmetry for a new cell. + void reset_symmetry(); + //-------------------------------------------------------------------------------- // getters const std::map>& get_irreducible_sector()const { return this->irs_.get_irreducible_sector(); } diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index f8c89b038ef..1a473c120de 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -192,7 +192,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) init_deltaspin_lcao(ucell, *this->inp_, &(this->pv), this->kv, this->p_hamilt, this->psi, this->dmat.dm, this->pelec); // 11) set xc type before the first cal of xc in pelec->init_scf, Peize Lin add 2016-12-03 - this->exx_nao.before_scf(ucell, this->kv, orb_, this->p_chgmix, istep, *this->inp_, this->exx_info_); + this->exx_nao.before_scf(ucell, this->kv, orb_, this->pv, this->p_chgmix, istep, *this->inp_, this->exx_info_); // 12) initalize DM(R), which has the same size with Hamiltonian(R) auto* hamilt_lcao = dynamic_cast*>(this->p_hamilt); diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index a31cb15b0a6..e63ada82957 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -239,7 +239,7 @@ void ESolver_KS_LCAO_TDDFT::runner(BaseCell& basecell, const int ist GlobalV::ofs_warning, atomic_rho_cfg_tddft, PARAM.globalv.has_float_data); - this->exx_nao.before_scf(ucell, this->kv, this->orb_, this->p_chgmix, totstep, *this->inp_, this->exx_info_); + this->exx_nao.before_scf(ucell, this->kv, this->orb_, this->pv, this->p_chgmix, totstep, *this->inp_, this->exx_info_); elecstate::init_scf(ucell, this->Pgrid, this->sf.strucFac, diff --git a/source/source_lcao/module_ri/exx_lri.h b/source/source_lcao/module_ri/exx_lri.h index 67eaa321c28..9240ace5299 100644 --- a/source/source_lcao/module_ri/exx_lri.h +++ b/source/source_lcao/module_ri/exx_lri.h @@ -57,6 +57,11 @@ class Exx_LRI using TatomR = std::array; // tmp public: + using CoulombSettings = std::map>>>>; + Exx_LRI(const Exx_Info_RI& info_in) :info(info_in) {} Exx_LRI operator=(const Exx_LRI&) = delete; Exx_LRI operator=(Exx_LRI&&); @@ -67,6 +72,10 @@ class Exx_LRI const K_Vectors &kv_in, const LCAO_Orbitals& orb, const std::vector>>& abfs_in = {}); + /// Refresh changed Coulomb kernels while retaining the initialized basis and Gaunt table. + void refresh_coulomb(const UnitCell& ucell, + const LCAO_Orbitals& orb, + const CoulombSettings& settings); void init_spencer(const MPI_Comm& mpi_comm_in, const UnitCell& ucell, const K_Vectors& kv_in, @@ -128,10 +137,7 @@ class Exx_LRI std::map> exx_objs; //LRI_CV cv; RI::Exx exx_lri; - std::map>>>> coulomb_settings; + CoulombSettings coulomb_settings; void post_process_Hexx( std::map>> &Hexxs_io ) const; double post_process_Eexx(const double& Eexx_in) const; diff --git a/source/source_lcao/module_ri/exx_lri.hpp b/source/source_lcao/module_ri/exx_lri.hpp index 04a134a7748..c89336ba5cb 100644 --- a/source/source_lcao/module_ri/exx_lri.hpp +++ b/source/source_lcao/module_ri/exx_lri.hpp @@ -69,27 +69,62 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, } this->exx_objs.clear(); - this->coulomb_settings = RI_Util::update_coulomb_settings(this->info.coulomb_param, ucell, this->p_kv); - + const auto settings = RI_Util::update_coulomb_settings(this->info.coulomb_param, ucell, this->p_kv); this->MGT = std::make_shared(); - for(const auto &settings_list : this->coulomb_settings) + this->refresh_coulomb(ucell, orb, settings); + + ModuleBase::timer::end("Exx_LRI", "init"); +} + +template +void Exx_LRI::refresh_coulomb(const UnitCell& ucell, + const LCAO_Orbitals& orb, + const CoulombSettings& settings) +{ + ModuleBase::TITLE("Exx_LRI", "refresh_coulomb"); + ModuleBase::timer::start("Exx_LRI", "refresh_coulomb"); + + for (auto object = this->exx_objs.begin(); object != this->exx_objs.end();) { - this->exx_objs[settings_list.first].abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, settings_list.second.second, this->info.ccp_rmesh_times); - this->exx_objs[settings_list.first].cv.set_orbitals(ucell, orb, - this->lcaos, this->abfs, this->exx_objs[settings_list.first].abfs_ccp, - this->info.kmesh_times, this->MGT, settings_list.second.first ); - this->exx_objs[settings_list.first].cv.set_info_ri(&this->info); + if (settings.count(object->first) == 0) + { + object = this->exx_objs.erase(object); + } + else + { + ++object; + } + } + for (const auto& settings_list : settings) + { + const auto previous = this->coulomb_settings.find(settings_list.first); + const bool unchanged = previous != this->coulomb_settings.end() + && previous->second == settings_list.second; + if (unchanged && this->exx_objs.count(settings_list.first) != 0) + { + continue; + } + + // Discard interpolation and point-value caches belonging to the old kernel. + // In particular, reusing Matrix_Orbs maps would retain old entries inserted by init(). + this->exx_objs.erase(settings_list.first); + auto& object = this->exx_objs[settings_list.first]; + object.abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp( + this->abfs, settings_list.second.second, this->info.ccp_rmesh_times); + object.cv.set_orbitals(ucell, orb, this->lcaos, this->abfs, object.abfs_ccp, + this->info.kmesh_times, this->MGT, settings_list.second.first); + object.cv.set_info_ri(&this->info); if (settings_list.first == Conv_Coulomb_Pot_K::Coulomb_Method::Ewald) { const int evq_abfs_Lmax = this->abfs_Lmax_; - this->exx_objs[settings_list.first].evq.init(ucell, orb, - this->mpi_comm, this->p_kv, this->lcaos, this->abfs, - settings_list.second.second, this->MGT, this->info.ccp_rmesh_times, this->info.kmesh_times, - evq_abfs_Lmax); + object.evq.init(ucell, orb, this->mpi_comm, this->p_kv, this->lcaos, this->abfs, + settings_list.second.second, this->MGT, this->info.ccp_rmesh_times, + this->info.kmesh_times, evq_abfs_Lmax); } } + this->coulomb_settings = settings; - ModuleBase::timer::end("Exx_LRI", "init"); + ModuleBase::timer::end("Exx_LRI", "refresh_coulomb"); } template diff --git a/source/source_lcao/module_ri/exx_lri_interface.h b/source/source_lcao/module_ri/exx_lri_interface.h index 888f5ad91a6..8b5bc61cf64 100644 --- a/source/source_lcao/module_ri/exx_lri_interface.h +++ b/source/source_lcao/module_ri/exx_lri_interface.h @@ -89,13 +89,15 @@ class Exx_LRI_Interface void cal_exx_stress(const double& omega, const double& lat0); // Processes in ESolver_KS_LCAO - /// @brief in before_all_runners: set symmetry according to irreducible k-points - /// since k-points are not reduced again after the variation of the cell and exx-symmetry must be consistent with k-points. - /// In the future, we will reduce k-points again during cell-relax, then this setting can be moved to `exx_beforescf`. + /// Initialize EXX symmetry before the first ionic step. void exx_before_all_runners(const K_Vectors& kv, const UnitCell& ucell, const Parallel_2D& pv); + /// Discard old atom mappings, lattice returns and rotation caches before rebuilding. + void refresh_symmetry(const K_Vectors& kv, const UnitCell& ucell, const Parallel_2D& pv); + /// @brief in beforescf: set xc type, opt_orb, do DM mixing - void exx_beforescf(const int istep, const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const LCAO_Orbitals& orb); + void exx_beforescf(const int istep, const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const LCAO_Orbitals& orb, + const Parallel_2D& pv, const bool update_symmetry); /// @brief in eachiterinit: do DM mixing and calculate Hexx when entering 2nd SCF void exx_eachiterinit(const int istep, diff --git a/source/source_lcao/module_ri/exx_lri_interface.hpp b/source/source_lcao/module_ri/exx_lri_interface.hpp index c09592f5fca..9bb4a45492a 100644 --- a/source/source_lcao/module_ri/exx_lri_interface.hpp +++ b/source/source_lcao/module_ri/exx_lri_interface.hpp @@ -85,6 +85,15 @@ void Exx_LRI_Interface::exx_before_all_runners( const Parallel_2D& pv) { ModuleBase::TITLE("Exx_LRI_Interface","exx_before_all_runners"); + this->refresh_symmetry(kv, ucell, pv); +} + +template +void Exx_LRI_Interface::refresh_symmetry( + const K_Vectors& kv, const UnitCell& ucell, const Parallel_2D& pv) +{ + ModuleBase::TITLE("Exx_LRI_Interface", "refresh_symmetry"); + this->symrot_.reset_symmetry(); // initialize the rotation matrix in AO representation this->exx_spacegroup_symmetry = (ModuleSymmetry::Symmetry::symm_flag == 1); if (this->exx_spacegroup_symmetry) @@ -103,12 +112,19 @@ void Exx_LRI_Interface::exx_beforescf(const int istep, const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, - const LCAO_Orbitals& orb) + const LCAO_Orbitals& orb, + const Parallel_2D& pv, + const bool update_symmetry) { ModuleBase::TITLE("Exx_LRI_Interface","exx_beforescf"); #ifdef __MPI if (this->info_global.cal_exx) { + if (update_symmetry) + { + this->refresh_symmetry(kv, ucell, pv); + } + if ((GlobalC::restart.info_load.load_H_finish && !GlobalC::restart.info_load.restart_exx) || (istep > 0) || (PARAM.inp.init_wfc == "file")) @@ -120,6 +136,18 @@ void Exx_LRI_Interface::exx_beforescf(const int istep, XC_Functional::set_xc_first_loop(ucell); } + // Spencer's cutoff depends on the current cell volume. Rebuild the + // convolved orbitals and their integral tables when its parameters + // change: the tables depend on abfs_ccp, not just on the unchanged + // atomic basis. Retain the basis and Gaunt table while invalidating + // interpolation and point-value caches belonging to changed kernels. + const auto coulomb_settings = RI_Util::update_coulomb_settings( + this->exx_ptr->get_info_ri().coulomb_param, ucell, &kv); + if (coulomb_settings != this->exx_ptr->coulomb_settings) + { + this->exx_ptr->refresh_coulomb(ucell, orb, coulomb_settings); + } + this->cal_exx_ions(ucell,PARAM.inp.out_ri_cv); } @@ -162,10 +190,10 @@ void Exx_LRI_Interface::exx_eachiterinit(const int istep, || istep > 0 || PARAM.inp.init_wfc == "file") // non separate loop case || (this->info_global.separate_loop - && PARAM.inp.init_wfc == "file" + && (istep > 0 || PARAM.inp.init_wfc == "file") && this->two_level_step == 0 && iter == 1) - ) // the first iter in separate loop case + ) // refresh the initial hybrid Hamiltonian for the current geometry { bool flag_restart = (iter == 1) ? true : false; @@ -346,7 +374,9 @@ bool Exx_LRI_Interface::exx_after_converge( auto restart_reset = [this]() { // avoid calling restart related procedure in the subsequent ion steps GlobalC::restart.info_load.restart_exx = true; - this->exx_ptr->Eexx = 0; + // Keep the energy paired with Hexxs until cal_exx_elec replaces both. + // Clearing only Eexx makes the next ionic step use an exchange + // Hamiltonian without its corresponding energy contribution. }; // no separate_loop case diff --git a/source/source_lcao/module_ri/module_exx_symmetry/test/CMakeLists.txt b/source/source_lcao/module_ri/module_exx_symmetry/test/CMakeLists.txt index 46fe581b419..015c4fc4a16 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/test/CMakeLists.txt +++ b/source/source_lcao/module_ri/module_exx_symmetry/test/CMakeLists.txt @@ -6,5 +6,6 @@ AddTest( LIBS base device symmetry neighbor parameter SOURCES symm_rotation_test.cpp test_symm_rotation.cpp ../symm_rotation.cpp ../symm_rot_out.cpp ../../../../source_basis/module_ao/parallel_orbitals.cpp + ../../../../source_cell/parallel_kpoints.cpp -) \ No newline at end of file +) 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 64b8ad020ca..a0faf2f17ac 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 @@ -147,3 +147,100 @@ TEST(SymmetryDensityRestoration, SpinorDensity) { check_little_group_restoration(4); } + +namespace +{ +class CellSymmetryRotation : public ModuleSymmetry::Symmetry_rotation +{ + public: + const std::vector>>& rotations() const + { + return this->Ms_; + } +}; +} + +TEST(SymmetryDensityRestoration, RebuildAfterCellSymmetryAnalysis) +{ + // A common translation crosses the cell boundary for only part of the basis. + // The point group stays unchanged, but the lattice returns and Bloch phases change. + UnitCell cell; + std::vector atoms(2); + cell.atoms = atoms.data(); + cell.ntype = 2; + cell.nat = 3; + cell.lmax = 0; + cell.latvec = ModuleBase::Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1); + 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}; + atoms[0].label = "A"; + atoms[0].na = 2; + atoms[0].nw = 1; + atoms[0].iw2l = {0}; + atoms[0].stapos_wf = 0; + atoms[0].taud = {{0.1, 0.2, 0.3}, {0.4, 0.3, 0.2}}; + atoms[0].tau = atoms[0].taud; + atoms[1].label = "B"; + atoms[1].na = 1; + atoms[1].nw = 1; + atoms[1].iw2l = {0}; + atoms[1].stapos_wf = 2; + atoms[1].taud = {{0.25, 0.25, 0.25}}; + atoms[1].tau = atoms[1].taud; + 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); + const int old_operations = cell.symm.nrotk; + K_Vectors kv; + kv.set_nks(1); + kv.set_nkstot(1); + kv.set_nkstot_nospin(2); + int total_kpoints = 1; + kv.para_k.kinfo(total_kpoints, 1, 0, 0, 1, 1); + kv.ik2iktot = {0}; + kv.kvec_d = {{0.25, 0, 0}}; + kv.kstars = {{{0, {0.25, 0, 0}}}}; + for (int operation = 0; operation < cell.symm.nrotk; ++operation) + { + kv.kstars[0][operation] = kv.kvec_d[0] * cell.symm.kgmatrix[operation]; + } + Parallel_2D pv; + pv.init(3, 3, 1, MPI_COMM_WORLD); + 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); + 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); + 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); + 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()); + EXPECT_EQ(reused.rotations(), fresh.rotations()); + EXPECT_NE(old_rotations, fresh.rotations()); + for (int atom = 0; atom < cell.nat; ++atom) + { + for (int operation = 0; operation < cell.symm.nrotk; ++operation) + { + const auto actual = reused.get_return_lattice(atom, operation); + const auto expected = fresh.get_return_lattice(atom, operation); + EXPECT_EQ(actual.x, expected.x); + EXPECT_EQ(actual.y, expected.y); + EXPECT_EQ(actual.z, expected.z); + } + } +} diff --git a/source/source_lcao/setup_exx.cpp b/source/source_lcao/setup_exx.cpp index c1371d33b96..0c17b3d2d9b 100644 --- a/source/source_lcao/setup_exx.cpp +++ b/source/source_lcao/setup_exx.cpp @@ -97,6 +97,7 @@ void Exx_NAO::before_scf( const UnitCell &ucell, // unitcell const K_Vectors &kv, const LCAO_Orbitals &orb, // orbital info + const Parallel_Orbitals& pv, Charge_Mixing* p_chgmix, const int istep, const Input_para& inp, @@ -105,13 +106,16 @@ void Exx_NAO::before_scf( #ifdef __EXX if (inp.calculation != "nscf") { + const bool update_symmetry = inp.calculation == "cell-relax" + && istep > 0 + && ModuleSymmetry::Symmetry::symm_flag == 1; if (exx_info.info_ri.real_number) { - this->exd->exx_beforescf(istep, kv, *p_chgmix, ucell, orb); + this->exd->exx_beforescf(istep, kv, *p_chgmix, ucell, orb, pv, update_symmetry); } else { - this->exc->exx_beforescf(istep, kv, *p_chgmix, ucell, orb); + this->exc->exx_beforescf(istep, kv, *p_chgmix, ucell, orb, pv, update_symmetry); } } else diff --git a/source/source_lcao/setup_exx.h b/source/source_lcao/setup_exx.h index 11f1fb52d96..caf1f4c8d0e 100644 --- a/source/source_lcao/setup_exx.h +++ b/source/source_lcao/setup_exx.h @@ -43,11 +43,13 @@ class Exx_NAO const UnitCell &ucell, // unitcell const K_Vectors &kv, const LCAO_Orbitals &orb, // orbital info + const Parallel_Orbitals& pv, Charge_Mixing* p_chgmix, const int istep, const Input_para& inp, Exx_Info& exx_info); + }; From 77654992cc1b3b3352cd906e9d72e398c09b2d4e Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Thu, 8 Oct 2026 13:28:55 -0400 Subject: [PATCH 2/7] fix(symmetry): preserve coordinate representatives in lattice returns Compute the integer return vector directly from r*g+t-r_mapped. Wrapping and tolerance-based snapping of fractional coordinates discarded integer shifts near cell boundaries, while RI integrals retained the original representatives. This corrupted Bloch phases and real-space EXX restoration, even for a fresh singlepoint when cell-relax enlarged symmetry_prec. Add tests for near-boundary atoms, changes of atom representatives, and translation representatives. Update the existing negative-coordinate test: z=-0.2 must remain -0.2, so its z return is zero rather than the old wrapped value -1. Also verify reconstruction of the original mapped coordinate. No INPUT options or defaults change; parameter documentation is unnecessary. Verification: - cmake --build build --target MODULE_RI_EXX_SYMMETRY_rotation -j2: passed. - OMP_NUM_THREADS=1 ctest --test-dir build -V -R '^MODULE_RI_EXX_SYMMETRY_rotation$': all 15 tests passed outside sandbox. - Initial test run exposed the obsolete wrapped-coordinate expectation; the intentional semantic update and coordinate-reconstruction check are described above. No tolerance was relaxed. - git diff --cached --check and staged agent governance: passed, no blockers; documentation warning explained above. - HPC isolated intervention on 72353ee12, 4 MPI ranks x 48 threads, 8-atom diamond step-two geometry, symmetry_prec=3.2e-5: first EXX energy changed from -29976.5115 to -1239.36097 eV; next two iterations -1239.36111 eV. Native EXX relative non-Hermiticity at k=1 decreased from 0.467132 to 6.25019e-5. Base Hamiltonian remained Hermitian to about 2e-16. Diagnostic full build passed, executable version v3.11.0-beta10. Full hybrid SCF and cell-relax convergence were not verified. --- .../module_symmetry/irreducible_sector.cpp | 29 +++----------- .../module_symmetry/irreducible_sector.h | 2 +- .../test/symm_rotation_test.cpp | 8 +++- .../test/test_symm_rotation.cpp | 40 +++++++++++++++++++ 4 files changed, 53 insertions(+), 26 deletions(-) diff --git a/source/source_cell/module_symmetry/irreducible_sector.cpp b/source/source_cell/module_symmetry/irreducible_sector.cpp index c6521c7ab52..2b9ca66d7bf 100644 --- a/source/source_cell/module_symmetry/irreducible_sector.cpp +++ b/source/source_cell/module_symmetry/irreducible_sector.cpp @@ -68,33 +68,14 @@ namespace ModuleSymmetry const ModuleBase::Matrix3& gmatd, const TCdouble gtransd, const TCdouble& posd_a1, const TCdouble& posd_a2)const { - // auto restrict_center = [&symm](const TCdouble& v) -> TCdouble { - // // in [-0.5, 0.5) - // TCdouble vr; - // vr.x = fmod(v.x + 100.5 + 0.5 * symm.epsilon, 1) - 0.5 - 0.5 * symm.epsilon; - // vr.y = fmod(v.y + 100.5 + 0.5 * symm.epsilon, 1) - 0.5 - 0.5 * symm.epsilon; - // vr.z = fmod(v.z + 100.5 + 0.5 * symm.epsilon, 1) - 0.5 - 0.5 * symm.epsilon; - // if (std::abs(vr.x) < symm.epsilon) vr.x = 0.0; - // if (std::abs(vr.y) < symm.epsilon) vr.y = 0.0; - // if (std::abs(vr.z) < symm.epsilon) vr.z = 0.0; - // return vr; - // }; - auto restrict_center = [&symm](const TCdouble& v) -> TCdouble { - // in [0,1) - TCdouble vr; - vr.x = fmod(v.x + 100 + symm.epsilon, 1) - symm.epsilon; - vr.y = fmod(v.y + 100 + symm.epsilon, 1) - symm.epsilon; - vr.z = fmod(v.z + 100 + symm.epsilon, 1) - symm.epsilon; - if (std::abs(vr.x) < symm.epsilon) vr.x = 0.0; - if (std::abs(vr.y) < symm.epsilon) vr.y = 0.0; - if (std::abs(vr.z) < symm.epsilon) vr.z = 0.0; - return vr; - }; + // Keep the representatives used by the orbital integrals. Wrapping an + // atom near a cell boundary without retaining its integer shift changes + // the Bloch phase and the real-space atom-pair mapping. auto check_integer = [&symm](const double x) -> void { assert(symm.equal(x, std::round(x))); }; - TCdouble rotpos1 = restrict_center(posd_a1) * gmatd + restrict_center(gtransd); // row vector - TCdouble return_lattice_double = rotpos1 - restrict_center(posd_a2); + TCdouble rotpos1 = posd_a1 * gmatd + gtransd; // row vector + TCdouble return_lattice_double = rotpos1 - posd_a2; #ifdef __DEBUG check_integer(return_lattice_double.x); check_integer(return_lattice_double.y); diff --git a/source/source_cell/module_symmetry/irreducible_sector.h b/source/source_cell/module_symmetry/irreducible_sector.h index f35be6a9d2f..73a9248130e 100644 --- a/source/source_cell/module_symmetry/irreducible_sector.h +++ b/source/source_cell/module_symmetry/irreducible_sector.h @@ -69,7 +69,7 @@ namespace ModuleSymmetry /// Perfoming {R|t} to atom position r in the R=0 lattice, we get Rr+t, which may get out of R=0 lattice, /// whose image in R=0 lattice is r'=Rr+t-O. This function is to get O for each atom and each symmetry operation. - /// the range of direct position is [-0.5, 0.5). + /// Preserve the input coordinate representatives; do not wrap away their integer lattice shifts. TCdouble get_return_lattice(const Symmetry& symm, const ModuleBase::Matrix3& gmatd, const TCdouble gtransd, const TCdouble& posd_a1, const TCdouble& posd_a2)const; diff --git a/source/source_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..ab93816a077 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 @@ -136,7 +136,13 @@ TEST_F(SymmetryRotationTest, GetReturnLattice) ModuleBase::Vector3 return_lattice = symrot.get_return_lattice(ModuleSymmetry::Symmetry(), gmatd, gtransd, posd_a1, posd_a2); EXPECT_NEAR(return_lattice.x, -1, DOUBLETHRESHOLD); EXPECT_NEAR(return_lattice.y, 0, DOUBLETHRESHOLD); - EXPECT_NEAR(return_lattice.z, -1, DOUBLETHRESHOLD); + // The target keeps z=-0.2 rather than wrapping it to 0.8: no z shift is needed. + EXPECT_NEAR(return_lattice.z, 0, DOUBLETHRESHOLD); + const auto rotated_position = posd_a1 * gmatd + gtransd; + const auto restored_position = rotated_position - return_lattice; + EXPECT_NEAR(restored_position.x, posd_a2.x, DOUBLETHRESHOLD); + EXPECT_NEAR(restored_position.y, posd_a2.y, DOUBLETHRESHOLD); + EXPECT_NEAR(restored_position.z, posd_a2.z, DOUBLETHRESHOLD); } TEST_F(SymmetryRotationTest, SetBlockToMat2d) 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 a0faf2f17ac..cee79ba3fb6 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 @@ -244,3 +244,43 @@ TEST(SymmetryDensityRestoration, RebuildAfterCellSymmetryAnalysis) } } } + +TEST(SymmetryDensityRestoration, ReturnLatticePreservesBoundaryRepresentatives) +{ + ModuleSymmetry::Symmetry symmetry; + symmetry.epsilon = 3.2e-5; + ModuleSymmetry::Irreducible_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); + const ModuleBase::Vector3 mapped(0.00001, 0.25, 0.5); + const auto lattice = sector.get_return_lattice(symmetry, reflection, translation, source, mapped); + EXPECT_DOUBLE_EQ(lattice.x, -1.0); + EXPECT_DOUBLE_EQ(lattice.y, 0.0); + EXPECT_DOUBLE_EQ(lattice.z, 0.0); + + // Changing either atom's representative must change its integer lattice shift. + const ModuleBase::Vector3 shifted_source(1.99999, 0.25, 0.5); + const ModuleBase::Vector3 shifted_mapped(1.00001, 0.25, 0.5); + const auto source_lattice = sector.get_return_lattice( + symmetry, reflection, translation, shifted_source, mapped); + const auto mapped_lattice = sector.get_return_lattice( + symmetry, reflection, translation, source, shifted_mapped); + EXPECT_DOUBLE_EQ(source_lattice.x, -2.0); + EXPECT_DOUBLE_EQ(mapped_lattice.x, -2.0); +} + +TEST(SymmetryDensityRestoration, ReturnLatticePreservesTranslationRepresentative) +{ + ModuleSymmetry::Symmetry symmetry; + symmetry.epsilon = 3.2e-5; + ModuleSymmetry::Irreducible_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); + const ModuleBase::Vector3 mapped(0.24999, 0.5, 0.75); + const auto lattice = sector.get_return_lattice(symmetry, identity, translation, source, mapped); + EXPECT_DOUBLE_EQ(lattice.x, 1.0); + EXPECT_DOUBLE_EQ(lattice.y, 0.0); + EXPECT_DOUBLE_EQ(lattice.z, 0.0); +} From 2482f809e5de1b19086f82ba9059b1101f2945c1 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Thu, 8 Oct 2026 14:59:58 -0400 Subject: [PATCH 3/7] fix(cell): refresh k-star operation indices after cell changes Reindex existing global k-stars in set_after_vc using the newly analyzed symmetry operations, including time reversal. Preserve IBZ representatives, weights, spin and pool ownership so the previous ionic-step density matrix can be used with the refreshed EXX rotations. Pass symmetry explicitly from the existing ESolver_FP caller; no new ESolver functions or global dependencies. Verification: - cmake --build build --target abacus_std_para -j2: passed after removing a pre-existing zero-byte read_atoms_helper.cpp.o from the build cache. - OMP_NUM_THREADS=1 ctest --test-dir build -V -R '^MODULE_CELL_klist_test(_para[14])?$': all three targets passed (37 serial cases, 2 parallel cases with one and four ranks). - git diff --check and working-tree governance checks: no blockers. - No INPUT names, defaults or parsing behavior changed; no INPUT docs needed. - Existing global logger arguments retained with net reference change zero. - HPC cell-relax regression will follow deployment. --- source/source_cell/klist.cpp | 58 ++++++++++++- source/source_cell/klist.h | 10 ++- source/source_cell/test/klist_test.cpp | 90 ++++++++++++++++++++- source/source_cell/test/klist_test_para.cpp | 2 +- source/source_esolver/esolver_fp.cpp | 5 +- 5 files changed, 158 insertions(+), 7 deletions(-) diff --git a/source/source_cell/klist.cpp b/source/source_cell/klist.cpp index c0be2ebc946..47981d578e9 100644 --- a/source/source_cell/klist.cpp +++ b/source/source_cell/klist.cpp @@ -565,10 +565,66 @@ void K_Vectors::reduce_by_symmetry(const UnitCell& ucell, return; } -void K_Vectors::set_after_vc(const ModuleBase::Matrix3& G, std::ofstream& ofs_running) +void K_Vectors::set_after_vc(const ModuleBase::Matrix3& G, + const ModuleSymmetry::Symmetry& symm, + const bool update_symmetry, + std::ofstream& ofs_running) { ofs_running << "\n SETUP K-POINTS" << std::endl; + // Every pool holds the complete stars. Their identity entries retain the + // global IBZ representatives even when kvec_d contains only local points. + // Reassign operation indices without changing the density-matrix ordering. + if (update_symmetry && !this->kstars.empty()) + { + const int operation_capacity = 2 * symm.nrotk + symm.nrotk_anti; + std::vector operations(operation_capacity); + for (int operation = 0; operation < symm.nrotk; ++operation) + { + operations[operation] = symm.kgmatrix[operation]; + } + const int operation_count = KListIO::append_time_reversal_ops(symm, operations, symm.nrotk); + typedef std::map> KStar; + std::vector updated_stars(this->kstars.size()); + for (std::size_t ik = 0; ik < this->kstars.size(); ++ik) + { + const KStar& star = this->kstars[ik]; + const KStar::const_iterator identity = star.find(0); + if (identity == star.end()) + { + ModuleBase::WARNING_QUIT("K_Vectors::set_after_vc", "Missing k-star IBZ representative."); + } + const ModuleBase::Vector3& representative = identity->second; + for (const KStar::value_type& member : star) + { + bool matched = false; + for (int operation = 0; operation < operation_count; ++operation) + { + ModuleBase::Vector3 rotated = member.second * operations[operation]; + ModuleCell::restrict_kpt(rotated, symm.epsilon); + if (symm.equal(rotated.x, representative.x) + && symm.equal(rotated.y, representative.y) + && symm.equal(rotated.z, representative.z)) + { + updated_stars[ik].emplace(operation, member.second); + matched = true; + break; + } + } + if (!matched) + { + ModuleBase::WARNING_QUIT("K_Vectors::set_after_vc", + "New cell symmetry cannot preserve the existing k-star."); + } + } + if (updated_stars[ik].size() != star.size()) + { + ModuleBase::WARNING_QUIT("K_Vectors::set_after_vc", "New cell symmetry changed the k-star size."); + } + } + this->kstars.swap(updated_stars); + } + // set cartesian k vectors. this->kvec_d2c(G); diff --git a/source/source_cell/klist.h b/source/source_cell/klist.h index f02f18a7d18..b5f3c8f68d6 100644 --- a/source/source_cell/klist.h +++ b/source/source_cell/klist.h @@ -167,9 +167,17 @@ class K_Vectors : public ModuleCell::ReciprocalGrid * as up to date. The spin multiplicity is not touched: it was fixed by * set() and never changes during a run. * + * Reindex existing k-stars using the newly analyzed symmetry operations, + * preserving IBZ representatives, weights, and pool/spin ownership. + * * @param G The new reciprocal lattice matrix. + * @param symm The symmetry analyzed for the new cell. + * @param update_symmetry Whether crystal symmetry is enabled. */ - void set_after_vc(const ModuleBase::Matrix3& G, std::ofstream& ofs_running); + void set_after_vc(const ModuleBase::Matrix3& G, + const ModuleSymmetry::Symmetry& symm, + const bool update_symmetry, + std::ofstream& ofs_running); //==================================================================== // Test seam. diff --git a/source/source_cell/test/klist_test.cpp b/source/source_cell/test/klist_test.cpp index a9a546b2196..15565677e20 100644 --- a/source/source_cell/test/klist_test.cpp +++ b/source/source_cell/test/klist_test.cpp @@ -669,7 +669,7 @@ TEST_F(KlistTest, SetAfterVC) kv->kvec_c[0].x = 0; kv->kvec_c[0].y = 0; kv->kvec_c[0].z = 0; - kv->set_after_vc(ucell.G, ofs_running); + kv->set_after_vc(ucell.G, ucell.symm, false, ofs_running); EXPECT_TRUE(kv->kd_done); EXPECT_TRUE(kv->kc_done); @@ -690,7 +690,7 @@ TEST_F(KlistTest, PrintKlists) kv->kvec_c[0].x = 0; kv->kvec_c[0].y = 0; kv->kvec_c[0].z = 0; - kv->set_after_vc(ucell.G, ofs_running); + kv->set_after_vc(ucell.G, ucell.symm, false, ofs_running); EXPECT_TRUE(kv->kd_done); kv->print_klists(ofs_running); ofs_running.close(); @@ -1045,3 +1045,89 @@ TEST_F(KlistTest, IbzKpointCustomWeights) remove("tmp_klist_custom_weights"); } + +TEST_F(KlistTest, SetAfterVCReindexesStarsWithoutChangingIBZ) +{ + kv->set_spin_mult(2); + kv->set_nkstot(2); + kv->set_nks(2); + kv->renew_for_testing(2); + kv->kvec_d = {{0.125, 0.25, 0}, {0.125, 0.25, 0}}; + kv->wk = {0.5, 0.5}; + kv->isk = {0, 1}; + kv->ik2iktot = {0, 1}; + kv->kstars = {{{0, {0.125, 0.25, 0}}, {1, {-0.125, -0.25, 0}}, + {2, {0.25, 0.125, 0}}, {3, {-0.25, -0.125, 0}}}}; + ucell.symm.nrotk = 4; + ucell.symm.nrotk_anti = 0; + ucell.symm.magnetic_nspin4 = false; + ucell.symm.epsilon = 1e-6; + ucell.symm.kgmatrix[0] = ModuleBase::Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1); + ucell.symm.kgmatrix[1] = ModuleBase::Matrix3(-1, 0, 0, 0, -1, 0, 0, 0, -1); + ucell.symm.kgmatrix[2] = ModuleBase::Matrix3(0, -1, 0, -1, 0, 0, 0, 0, -1); + ucell.symm.kgmatrix[3] = ModuleBase::Matrix3(0, 1, 0, 1, 0, 0, 0, 0, 1); + const auto points = kv->kvec_d; + const auto weights = kv->wk; + const auto spins = kv->isk; + const auto ownership = kv->ik2iktot; + kv->set_after_vc(ucell.G, ucell.symm, true, ofs_running); + EXPECT_EQ(kv->kvec_d, points); + EXPECT_EQ(kv->wk, weights); + EXPECT_EQ(kv->isk, spins); + EXPECT_EQ(kv->ik2iktot, ownership); + ASSERT_EQ(kv->kstars.size(), 1); + const auto& star = kv->kstars[0]; + ASSERT_EQ(star.size(), 4); + EXPECT_DOUBLE_EQ(star.at(3).x, 0.25); + EXPECT_DOUBLE_EQ(star.at(2).x, -0.25); + for (const auto& member : star) + { + const auto rotated = member.second * ucell.symm.kgmatrix[member.first]; + EXPECT_NEAR(rotated.x, points[0].x, 1e-12); + EXPECT_NEAR(rotated.y, points[0].y, 1e-12); + EXPECT_NEAR(rotated.z, points[0].z, 1e-12); + } + const auto refreshed = kv->kstars; + kv->set_after_vc(ucell.G, ucell.symm, true, ofs_running); + EXPECT_EQ(kv->kstars, refreshed); +} + +TEST_F(KlistTest, SetAfterVCRejectsIncompatibleStars) +{ + kv->kstars = {{{0, {0.125, 0, 0}}, {1, {-0.125, 0, 0}}}}; + ucell.symm.nrotk = 1; + ucell.symm.nrotk_anti = 0; + ucell.symm.magnetic_nspin4 = true; + ucell.symm.epsilon = 1e-6; + ucell.symm.kgmatrix[0] = ModuleBase::Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1); + EXPECT_EXIT(kv->set_after_vc(ucell.G, ucell.symm, true, ofs_running), + ::testing::ExitedWithCode(1), ""); +} + +TEST_F(KlistTest, SetAfterVCUsesGlobalStarsAndTimeReversal) +{ + // Only the second IBZ point is local, while stars are replicated globally. + kv->set_spin_mult(1); + kv->set_nkstot(2); + kv->set_nks(1); + kv->renew_for_testing(2); + kv->kvec_d = {{0.25, 0, 0}}; + kv->ik2iktot = {1}; + kv->kstars = {{{0, {0.125, 0, 0}}, {7, {-0.125, 0, 0}}}, + {{0, {0.25, 0, 0}}, {7, {-0.25, 0, 0}}}}; + ucell.symm.nrotk = 1; + ucell.symm.nrotk_anti = 0; + ucell.symm.magnetic_nspin4 = false; + ucell.symm.epsilon = 1e-6; + ucell.symm.kgmatrix[0] = ModuleBase::Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1); + const auto original = kv->kstars; + kv->set_after_vc(ucell.G, ucell.symm, false, ofs_running); + EXPECT_EQ(kv->kstars, original); + kv->set_after_vc(ucell.G, ucell.symm, true, ofs_running); + ASSERT_EQ(kv->kstars.size(), 2); + EXPECT_DOUBLE_EQ(kv->kstars[0].at(0).x, 0.125); + EXPECT_DOUBLE_EQ(kv->kstars[0].at(1).x, -0.125); + EXPECT_DOUBLE_EQ(kv->kstars[1].at(0).x, 0.25); + EXPECT_DOUBLE_EQ(kv->kstars[1].at(1).x, -0.25); + EXPECT_EQ(kv->ik2iktot, std::vector{1}); +} diff --git a/source/source_cell/test/klist_test_para.cpp b/source/source_cell/test/klist_test_para.cpp index 0bd7552a853..7f8f76d524c 100644 --- a/source/source_cell/test/klist_test_para.cpp +++ b/source/source_cell/test/klist_test_para.cpp @@ -350,7 +350,7 @@ TEST_F(KlistParaTest, SetAfterVC) } // call set_after_vc here kv->kc_done = false; - kv->set_after_vc(ucell.G, GlobalV::ofs_running); + kv->set_after_vc(ucell.G, symm, true, GlobalV::ofs_running); EXPECT_TRUE(kv->kc_done); EXPECT_TRUE(kv->kd_done); // clear diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 5a661ce4a35..309f1a7bb7f 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -209,7 +209,8 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); // perform symmetry analysis - if (ModuleSymmetry::Symmetry::symm_flag == 1) + const bool update_symmetry = ModuleSymmetry::Symmetry::symm_flag == 1; + if (update_symmetry) { 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, @@ -218,7 +219,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) } // reset k-points - kv.set_after_vc(ucell.G, GlobalV::ofs_running); + kv.set_after_vc(ucell.G, ucell.symm, update_symmetry, GlobalV::ofs_running); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); } From 1115bc761546dbc05c811bae7a491e28c95d9842 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Thu, 8 Oct 2026 23:32:37 -0400 Subject: [PATCH 4/7] fix(symmetry): rebuild atom maps after each symmetry analysis --- .../source_cell/module_symmetry/symmetry.cpp | 7 +-- .../test/symm_test_analysis.cpp | 58 +++++++++++++++++++ 2 files changed, 61 insertions(+), 4 deletions(-) diff --git a/source/source_cell/module_symmetry/symmetry.cpp b/source/source_cell/module_symmetry/symmetry.cpp index c22774aeb02..d7b5f6f158a 100644 --- a/source/source_cell/module_symmetry/symmetry.cpp +++ b/source/source_cell/module_symmetry/symmetry.cpp @@ -11,12 +11,11 @@ bool Symmetry::pricell_loop = true; void Symmetry::set_atom_map(const Atom* atoms) { ModuleBase::TITLE("Symmetry", "set_atom_map"); - if (this->isym_rotiat_.size() == this->nrotk) { - return; - } + // A new analysis can reorder operations without changing their number. + // Rebuild their atom mappings from the current geometry every time. this->isym_rotiat_.resize(this->nrotk); for (int i = 0; i < this->nrotk; ++i) { - this->isym_rotiat_[i].resize(this->nat, -1); + this->isym_rotiat_[i].assign(this->nat, -1); } double* pos = this->newpos; diff --git a/source/source_cell/module_symmetry/test/symm_test_analysis.cpp b/source/source_cell/module_symmetry/test/symm_test_analysis.cpp index 51d6c452601..24dcbbd26d8 100644 --- a/source/source_cell/module_symmetry/test/symm_test_analysis.cpp +++ b/source/source_cell/module_symmetry/test/symm_test_analysis.cpp @@ -39,6 +39,64 @@ SepPot::~SepPot(){} Sep_Cell::Sep_Cell() noexcept {} Sep_Cell::~Sep_Cell() noexcept {} +TEST_F(SymmetryTest, AtomMapRefreshWithUnchangedOperationCount) +{ + stru_ diamond; + diamond.cell = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; + diamond.coordtype = "D"; + atomtype_ carbon; + carbon.atomname = "C"; + carbon.coordinate = {{0.0, 0.0, 0.0}, {0.0, 0.5, 0.5}, + {0.5, 0.0, 0.5}, {0.5, 0.5, 0.0}, + {0.25, 0.25, 0.25}, {0.25, 0.75, 0.75}, + {0.75, 0.25, 0.75}, {0.75, 0.75, 0.25}}; + diamond.all_type.push_back(carbon); + construct_ucell(diamond); + ModuleSymmetry::Symmetry symm; + const int representation[2] = {0, 0}; + const std::string calculation = "cell-relax"; + symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, + 1e-6, 1, calculation, representation); + const int original_count = symm.nrotk; + EXPECT_EQ(original_count, 48); + + // Keep the same crystal symmetry while changing atom representatives. + std::swap(ucell.atoms[0].taud[0], ucell.atoms[0].taud[1]); + ucell.latvec *= 0.98; + ucell.a1 *= 0.98; + ucell.a2 *= 0.98; + ucell.a3 *= 0.98; + for (int atom = 0; atom < ucell.nat; ++atom) + { + ucell.atoms[0].tau[atom] = ucell.atoms[0].taud[atom] * ucell.latvec; + } + symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, + 1e-6, 1, calculation, representation); + EXPECT_EQ(symm.nrotk, original_count); + for (int operation = 0; operation < symm.nrotk; ++operation) + { + for (int atom = 0; atom < ucell.nat; ++atom) + { + const int mapped = symm.get_rotated_atom(operation, atom); + EXPECT_GE(mapped, 0); + EXPECT_LT(mapped, ucell.nat); + if (mapped < 0 || mapped >= ucell.nat) + { + continue; + } + const auto rotated = ucell.atoms[0].taud[atom] * symm.gmatrix[operation]; + const auto residual = rotated + symm.gtrans[operation] - ucell.atoms[0].taud[mapped]; + const double lattice_x = std::round(residual.x); + const double lattice_y = std::round(residual.y); + const double lattice_z = std::round(residual.z); + EXPECT_NEAR(residual.x, lattice_x, 1e-8); + EXPECT_NEAR(residual.y, lattice_y, 1e-8); + EXPECT_NEAR(residual.z, lattice_z, 1e-8); + } + } + ClearUcell(); +} + TEST_F(SymmetryTest, AnalySys) { for (int stru = 0; stru < stru_lib.size(); stru++) From 52a180f8485fb4b38902260db2682ef8dd5ea946 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Fri, 9 Oct 2026 03:39:45 -0400 Subject: [PATCH 5/7] test(exx): update ionic-step HSE energy references Refresh energy references for the current-geometry initial EXX Hamiltonian. Keep the quick tests' inputs, thresholds, outer-loop limits and timing references unchanged. Verified all three cases with Autotest.sh, four MPI ranks and one OpenMP thread; all six energy comparisons passed. No INPUT or interface changes require documentation updates. --- tests/08_RI/md_hse_multik/result.ref | 4 ++-- tests/08_RI/relax_cell_hse_multik/result.ref | 4 ++-- tests/08_RI/relax_hse_gamma/result.ref | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/08_RI/md_hse_multik/result.ref b/tests/08_RI/md_hse_multik/result.ref index ff5a2d05be3..c1d1b8bc493 100644 --- a/tests/08_RI/md_hse_multik/result.ref +++ b/tests/08_RI/md_hse_multik/result.ref @@ -1,3 +1,3 @@ -etotref -467.270790968887 -etotperatomref -155.7569303230 +etotref -467.0419725776744 +etotperatomref -155.6806575259 totaltimeref 3.40 diff --git a/tests/08_RI/relax_cell_hse_multik/result.ref b/tests/08_RI/relax_cell_hse_multik/result.ref index b0a00725fce..cb786101f30 100644 --- a/tests/08_RI/relax_cell_hse_multik/result.ref +++ b/tests/08_RI/relax_cell_hse_multik/result.ref @@ -1,3 +1,3 @@ -etotref -463.4113534713501963 -etotperatomref -154.4704511571 +etotref -463.4160767981939557 +etotperatomref -154.4720255994 totaltimeref 5.46 diff --git a/tests/08_RI/relax_hse_gamma/result.ref b/tests/08_RI/relax_hse_gamma/result.ref index 925d8f9b61d..8787b747da6 100644 --- a/tests/08_RI/relax_hse_gamma/result.ref +++ b/tests/08_RI/relax_hse_gamma/result.ref @@ -1,3 +1,3 @@ -etotref -461.2147597707549 -etotperatomref -153.7382532569 +etotref -461.1879847890294 +etotperatomref -153.7293282630 totaltimeref 3.20 From 6e9e079a99f984d0e8286760e0c028e7a92367f7 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Fri, 9 Oct 2026 10:19:43 -0400 Subject: [PATCH 6/7] fix(exx): print SCF rerun timing in decimal scientific format Use std::scientific to replace the floating-point format rather than adding a scientific flag to an existing fixed flag. This avoids hexadecimal timing output after pressure logging. Verification: cmake --build build --target abacus_std_para -j2 passed. Autotest.sh with four MPI ranks and one OpenMP thread passed relax_hse_gamma and relax_cell_hse_multik (four energy checks). All eight EXX timings matched decimal scientific notation, including those after pressure output. No INPUT or interface changes require documentation updates. --- source/source_lcao/module_ri/exx_lri_interface.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/source_lcao/module_ri/exx_lri_interface.hpp b/source/source_lcao/module_ri/exx_lri_interface.hpp index 9bb4a45492a..13fe79697d5 100644 --- a/source/source_lcao/module_ri/exx_lri_interface.hpp +++ b/source/source_lcao/module_ri/exx_lri_interface.hpp @@ -451,7 +451,7 @@ bool Exx_LRI_Interface::exx_after_converge( timeval t_end; gettimeofday(&t_end, nullptr); std::cout << "and rerun SCF\t" - << std::setprecision(3) << std::setiosflags(std::ios::scientific) + << std::setprecision(3) << std::scientific << (double)(t_end.tv_sec-t_start.tv_sec) + (double)(t_end.tv_usec-t_start.tv_usec)/1000000.0 << std::defaultfloat << " (s)" << std::endl; return false; From 11b7b3cfb8cbc0319298fd429a7bb8ad32c8e1b6 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Fri, 9 Oct 2026 20:07:32 -0400 Subject: [PATCH 7/7] fix: compare k points modulo reciprocal lattice vectors Share periodic equivalence across initial IBZ reduction, rotation consistency checks, initial k-star construction and variable-cell k-star updates. Keep the epsilon-dependent normalization convention while allowing equivalent boundary representatives to match. Add regressions for reciprocal translations and tolerance changes at the k-point boundary. No INPUT defaults or tolerance semantics change; user-facing parameter documentation does not require an update. --- source/source_cell/klist.cpp | 10 +++--- source/source_cell/klist_io.cpp | 4 +-- source/source_cell/klist_io.h | 6 ++-- source/source_cell/reciprocal_grid.cpp | 25 +++++++++++---- source/source_cell/reciprocal_grid.h | 5 +++ source/source_cell/test/klist_test.cpp | 32 +++++++++++++++++++ .../source_cell/test/reciprocal_grid_test.cpp | 15 +++++++++ 7 files changed, 78 insertions(+), 19 deletions(-) diff --git a/source/source_cell/klist.cpp b/source/source_cell/klist.cpp index 47981d578e9..ef3809238c3 100644 --- a/source/source_cell/klist.cpp +++ b/source/source_cell/klist.cpp @@ -7,6 +7,7 @@ #include "source_base/parallel_reduce.h" #include "source_cell/module_symmetry/symmetry.h" + void K_Vectors::set(const UnitCell& ucell, const ModuleSymmetry::Symmetry& symm, const std::string& k_file_name, @@ -546,7 +547,6 @@ void K_Vectors::reduce_by_symmetry(const UnitCell& ucell, nrotkm, kvec_d_ibz, symm.epsilon, - [&symm](double a, double b) { return symm.equal(a, b); }, this->kstars); } @@ -600,11 +600,9 @@ void K_Vectors::set_after_vc(const ModuleBase::Matrix3& G, bool matched = false; for (int operation = 0; operation < operation_count; ++operation) { - ModuleBase::Vector3 rotated = member.second * operations[operation]; - ModuleCell::restrict_kpt(rotated, symm.epsilon); - if (symm.equal(rotated.x, representative.x) - && symm.equal(rotated.y, representative.y) - && symm.equal(rotated.z, representative.z)) + const ModuleBase::Vector3 rotated = member.second * operations[operation]; + // Representatives may retain the boundary convention of an older tolerance. + if (ModuleCell::kpoints_equivalent(rotated, representative, symm.epsilon)) { updated_stars[ik].emplace(operation, member.second); matched = true; diff --git a/source/source_cell/klist_io.cpp b/source/source_cell/klist_io.cpp index f869a0c7baf..9afa9a09eab 100644 --- a/source/source_cell/klist_io.cpp +++ b/source/source_cell/klist_io.cpp @@ -176,7 +176,6 @@ void build_kstars(const std::vector>& kvec_d, const int nrotkm, const std::vector>& kvec_d_ibz, const double epsilon, - const std::function& equal, std::vector>>& kstars) { const int nkstot = static_cast(kvec_d.size()); @@ -194,8 +193,7 @@ void build_kstars(const std::vector>& kvec_d, ModuleCell::restrict_kpt(kvec_rot, epsilon); for (int k = 0; k < nkstot_ibz; ++k) { - if (equal(kvec_rot.x, kvec_d_ibz[k].x) && equal(kvec_rot.y, kvec_d_ibz[k].y) - && equal(kvec_rot.z, kvec_d_ibz[k].z)) + if (ModuleCell::kpoints_equivalent(kvec_rot, kvec_d_ibz[k], epsilon)) { isym = j; exist_number = k; diff --git a/source/source_cell/klist_io.h b/source/source_cell/klist_io.h index c9785b0f85e..2e8cf82411d 100644 --- a/source/source_cell/klist_io.h +++ b/source/source_cell/klist_io.h @@ -5,7 +5,6 @@ #include "source_base/vector3.h" #include -#include #include #include #include @@ -61,15 +60,14 @@ LineK interp_line(std::ifstream& ifk, int nks_special); /// Build the EXX k-stars: for every k point, find the symmetry operation /// (index into `kgmatrix`) that rotates it onto an irreducible k point, and -/// group k points by that IBZ representative. `equal` compares two doubles -/// with the symmetry precision; `epsilon` is the k-restriction tolerance. +/// group k points by that IBZ representative using periodic equivalence. +/// `epsilon` is the symmetry and k-restriction tolerance. /// this-free so the heavy triple loop is isolated and testable. void build_kstars(const std::vector>& kvec_d, const std::vector& kgmatrix, int nrotkm, const std::vector>& kvec_d_ibz, double epsilon, - const std::function& equal, std::vector>>& kstars); /// Append the time-reversal-related k-point symmetry operations into diff --git a/source/source_cell/reciprocal_grid.cpp b/source/source_cell/reciprocal_grid.cpp index 10f5e1a32cc..c185a3bc1b5 100644 --- a/source/source_cell/reciprocal_grid.cpp +++ b/source/source_cell/reciprocal_grid.cpp @@ -6,6 +6,8 @@ */ #include "reciprocal_grid.h" +#include + #include "source_cell/unitcell.h" #include "source_cell/module_symmetry/symmetry.h" #include "source_base/global_function.h" @@ -18,6 +20,19 @@ namespace ModuleCell { +bool kpoints_equivalent(const ModuleBase::Vector3& first, + const ModuleBase::Vector3& second, + const double epsilon) +{ + ModuleBase::Vector3 difference = first - second; + difference.x -= std::round(difference.x); + difference.y -= std::round(difference.y); + difference.z -= std::round(difference.z); + return std::abs(difference.x) < epsilon + && std::abs(difference.y) < epsilon + && std::abs(difference.z) < epsilon; +} + void restrict_kpt(ModuleBase::Vector3& kvec, double epsilon) { // fold into (-0.5, 0.5]; the epsilon shift keeps points sitting on the @@ -329,7 +344,6 @@ void ReciprocalGrid::reduce_ibz(const ModuleBase::Matrix3* rot_ops, std::vector& ibz_index, std::vector& ibz2bz) { - auto equal = [epsilon](double m, double n) { return fabs(m - n) < epsilon; }; // direct coordinates of points in the k-lattice std::vector> kvec_d_k(this->nkstot); @@ -382,15 +396,14 @@ void ReciprocalGrid::reduce_ibz(const ModuleBase::Matrix3* rot_ops, kvec_rot_k = kvec_rot_k * k_lattice * G.Inverse(); // convert to recip lattice restrict_kpt(kvec_rot_k, epsilon); - assert(equal(kvec_rot.x, kvec_rot_k.x)); - assert(equal(kvec_rot.y, kvec_rot_k.y)); - assert(equal(kvec_rot.z, kvec_rot_k.z)); + const bool consistent = kpoints_equivalent(kvec_rot, kvec_rot_k, epsilon); + assert(consistent); + (void)consistent; kvec_rot_k = kvec_rot_k * G * k_lattice.Inverse(); // convert back to k-lattice } for (int k = 0; k < nkstot_ibz; ++k) { - if (equal(kvec_rot.x, kvec_d_ibz[k].x) && equal(kvec_rot.y, kvec_d_ibz[k].y) - && equal(kvec_rot.z, kvec_d_ibz[k].z)) + if (kpoints_equivalent(kvec_rot, kvec_d_ibz[k], epsilon)) { already_exist = true; // find another ibz point, diff --git a/source/source_cell/reciprocal_grid.h b/source/source_cell/reciprocal_grid.h index ced3a0b4ff5..219225d6b8d 100644 --- a/source/source_cell/reciprocal_grid.h +++ b/source/source_cell/reciprocal_grid.h @@ -36,6 +36,11 @@ namespace ModuleCell */ void restrict_kpt(ModuleBase::Vector3& kvec, double epsilon); +/// Compare reciprocal fractional coordinates modulo integer reciprocal vectors. +bool kpoints_equivalent(const ModuleBase::Vector3& first, + const ModuleBase::Vector3& second, + double epsilon); + /** * @brief Abstract base class shared by K_Vectors (electrons) and QList (phonons). * diff --git a/source/source_cell/test/klist_test.cpp b/source/source_cell/test/klist_test.cpp index 15565677e20..62aaef39837 100644 --- a/source/source_cell/test/klist_test.cpp +++ b/source/source_cell/test/klist_test.cpp @@ -1092,6 +1092,38 @@ TEST_F(KlistTest, SetAfterVCReindexesStarsWithoutChangingIBZ) EXPECT_EQ(kv->kstars, refreshed); } +TEST_F(KlistTest, SetAfterVCMatchesPeriodicBoundaryWhenToleranceChanges) +{ + kv->set_spin_mult(1); + kv->set_nkstot(1); + kv->set_nks(1); + kv->renew_for_testing(1); + kv->kvec_d = {{-0.499999, 0.125, -0.499999}}; + kv->wk = {1.0}; + kv->ik2iktot = {0}; + kv->kstars = {{{0, {-0.499999, 0.125, -0.499999}}, + {1, {-0.499999, -0.125, -0.499999}}}}; + ucell.symm.nrotk = 2; + ucell.symm.nrotk_anti = 0; + ucell.symm.magnetic_nspin4 = true; + ucell.symm.kgmatrix[0] = ModuleBase::Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1); + ucell.symm.kgmatrix[1] = ModuleBase::Matrix3(1, 0, 0, 0, -1, 0, 0, 0, 1); + ucell.symm.epsilon = 1e-6; + kv->set_after_vc(ucell.G, ucell.symm, true, ofs_running); + const auto stars = kv->kstars; + const auto points = kv->kvec_d; + const auto weights = kv->wk; + + // The new tolerance folds -0.499999 to 0.500001. Both are the same k point. + ucell.symm.epsilon = 3.2e-5; + kv->set_after_vc(ucell.G, ucell.symm, true, ofs_running); + EXPECT_EQ(kv->kstars, stars); + EXPECT_EQ(kv->kvec_d, points); + EXPECT_EQ(kv->wk, weights); + kv->set_after_vc(ucell.G, ucell.symm, true, ofs_running); + EXPECT_EQ(kv->kstars, stars); +} + TEST_F(KlistTest, SetAfterVCRejectsIncompatibleStars) { kv->kstars = {{{0, {0.125, 0, 0}}, {1, {-0.125, 0, 0}}}}; diff --git a/source/source_cell/test/reciprocal_grid_test.cpp b/source/source_cell/test/reciprocal_grid_test.cpp index 05525749d38..892c6935d9f 100644 --- a/source/source_cell/test/reciprocal_grid_test.cpp +++ b/source/source_cell/test/reciprocal_grid_test.cpp @@ -319,3 +319,18 @@ TEST_F(ReciprocalGridTest, ReduceIbzMpKLattice) } EXPECT_NEAR(sum, 1.0, 1e-12); } + +TEST(ReciprocalGrid, KpointsEquivalentModuloReciprocalVectors) +{ + const ModuleBase::Vector3 first(-0.499999, 0.125, -0.25); + const ModuleBase::Vector3 translated(0.500001, -1.875, 2.75); + const ModuleBase::Vector3 near(0.5000012, -1.875, 2.75); + const ModuleBase::Vector3 different(0.500011, -1.875, 2.75); + const double epsilon = 1e-6; + EXPECT_TRUE(ModuleCell::kpoints_equivalent(first, translated, epsilon)); + EXPECT_TRUE(ModuleCell::kpoints_equivalent(translated, first, epsilon)); + EXPECT_TRUE(ModuleCell::kpoints_equivalent(first, near, epsilon)); + EXPECT_FALSE(ModuleCell::kpoints_equivalent(first, different, epsilon)); + EXPECT_DOUBLE_EQ(first.x, -0.499999); + EXPECT_DOUBLE_EQ(translated.x, 0.500001); +}