diff --git a/source/source_cell/klist.cpp b/source/source_cell/klist.cpp index c0be2ebc946..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); } @@ -565,10 +565,64 @@ void K_Vectors::reduce_by_symmetry(const UnitCell& ucell, return; } -void K_Vectors::set_after_vc(const ModuleBase::Matrix3& G, std::ofstream& ofs_running) +void K_Vectors::set_after_vc(const ModuleBase::Matrix3& G, + const ModuleSymmetry::Symmetry& symm, + const bool update_symmetry, + std::ofstream& ofs_running) { ofs_running << "\n SETUP K-POINTS" << std::endl; + // Every pool holds the complete stars. Their identity entries retain the + // global IBZ representatives even when kvec_d contains only local points. + // Reassign operation indices without changing the density-matrix ordering. + if (update_symmetry && !this->kstars.empty()) + { + const int operation_capacity = 2 * symm.nrotk + symm.nrotk_anti; + std::vector 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) + { + 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; + 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/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/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_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_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++) 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 a9a546b2196..62aaef39837 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,121 @@ 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, 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}}}}; + 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_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); +} 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"); } 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..13fe79697d5 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 @@ -421,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; 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/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 64b8ad020ca..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 @@ -147,3 +147,140 @@ 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); + } + } +} + +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); +} 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); + }; 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