From 74a799b3c9b66aad9ad43de7476e8cb03afa6383 Mon Sep 17 00:00:00 2001 From: Yanglijf Date: Sat, 6 Jun 2026 19:31:34 +0800 Subject: [PATCH 1/2] Add AtomPack default query path and half-domain neighbor prototype --- .../source_cell/module_neighbor/atom_pack.cpp | 260 ++++++++++++++++++ .../source_cell/module_neighbor/atom_pack.h | 30 ++ .../source_cell/module_neighbor/sltk_grid.cpp | 54 +++- .../source_cell/module_neighbor/sltk_grid.h | 28 +- .../module_neighbor/sltk_grid_driver.cpp | 52 +++- .../module_neighbor/sltk_grid_driver.h | 16 ++ .../module_neighbor/test/CMakeLists.txt | 4 +- .../module_neighbor/test/atom_pack_test.cpp | 163 +++++++++++ .../test/sltk_atom_arrange_test.cpp | 153 +++++++++++ .../module_neighbor/test/sltk_grid_test.cpp | 45 +++ 10 files changed, 798 insertions(+), 7 deletions(-) diff --git a/source/source_cell/module_neighbor/atom_pack.cpp b/source/source_cell/module_neighbor/atom_pack.cpp index 17feddc2819..93874c8dffb 100644 --- a/source/source_cell/module_neighbor/atom_pack.cpp +++ b/source/source_cell/module_neighbor/atom_pack.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -79,6 +80,113 @@ int clamp_index(const int index, const int size) } return index; } + +void box_id_to_indices(const GridStorage& storage, const int box_id, int& bx, int& by, int& bz) +{ + const int yz_size = storage.box_ny * storage.box_nz; + bx = box_id / yz_size; + const int yz_id = box_id % yz_size; + by = yz_id / storage.box_nz; + bz = yz_id % storage.box_nz; +} + +bool is_center_atom(const AtomPack& pack, const int index) +{ + return !pack.is_ghost[index] && pack.cell_x[index] == 0 && pack.cell_y[index] == 0 && pack.cell_z[index] == 0; +} + +using AtomImageKey = std::tuple; + +AtomImageKey make_atom_image_key(const AtomPack& pack, const int index) +{ + return std::make_tuple(pack.type[index], pack.natom[index], pack.cell_x[index], pack.cell_y[index], pack.cell_z[index]); +} + +std::map build_atom_image_index(const AtomPack& pack) +{ + std::map index_by_image; + for (int i = 0; i < pack.size(); ++i) + { + if (!pack.is_ghost[i]) + { + index_by_image[make_atom_image_key(pack, i)] = i; + } + } + return index_by_image; +} + +int find_atom_image_index(const std::map& index_by_image, + const int type, + const int natom, + const int cell_x, + const int cell_y, + const int cell_z) +{ + const auto iter = index_by_image.find(std::make_tuple(type, natom, cell_x, cell_y, cell_z)); + return iter == index_by_image.end() ? -1 : iter->second; +} + +void validate_neighbor_search_input(const GridStorage& storage, const double radius) +{ + if (radius < 0.0) + { + throw std::invalid_argument("Neighbor search radius must be non-negative."); + } + if (storage.box_edge_length <= 0.0 || storage.box_size() <= 0) + { + throw std::runtime_error("GridStorage has not been initialized."); + } +} + +NeighborPair make_neighbor_pair(const AtomPack& pack, const int center, const int candidate) +{ + return NeighborPair{pack.type[center], + pack.natom[center], + pack.type[candidate], + pack.natom[candidate], + pack.cell_x[candidate], + pack.cell_y[candidate], + pack.cell_z[candidate], + center, + candidate}; +} + +NeighborPair make_restored_reverse_pair(const AtomPack& pack, + const NeighborPair& pair, + const std::map& index_by_image) +{ + const int reverse_center = find_atom_image_index(index_by_image, pair.neighbor_type, pair.neighbor_natom, 0, 0, 0); + const int reverse_neighbor = find_atom_image_index(index_by_image, + pair.center_type, + pair.center_natom, + -pair.cell_x, + -pair.cell_y, + -pair.cell_z); + + return NeighborPair{pair.neighbor_type, + pair.neighbor_natom, + pair.center_type, + pair.center_natom, + -pair.cell_x, + -pair.cell_y, + -pair.cell_z, + reverse_center, + reverse_neighbor}; +} + +bool is_within_radius(const AtomPack& pack, const int center, const int candidate, const double radius2) +{ + const double dx = pack.x[center] - pack.x[candidate]; + const double dy = pack.y[center] - pack.y[candidate]; + const double dz = pack.z[center] - pack.z[candidate]; + const double dr = dx * dx + dy * dy + dz * dz; + return dr != 0.0 && dr <= radius2; +} + +bool is_half_domain_offset(const int dbx, const int dby, const int dbz) +{ + return dbx > 0 || (dbx == 0 && dby > 0) || (dbx == 0 && dby == 0 && dbz >= 0); +} } // namespace void AtomPack::clear() @@ -202,6 +310,21 @@ int GridStorage::get_box_id_from_coord(const double x, const double y, const dou return get_box_id(bx, by, bz); } +std::tuple NeighborPair::key() const +{ + return std::make_tuple(center_type, center_natom, neighbor_type, neighbor_natom, cell_x, cell_y, cell_z); +} + +bool NeighborPair::operator<(const NeighborPair& rhs) const +{ + return key() < rhs.key(); +} + +bool NeighborPair::operator==(const NeighborPair& rhs) const +{ + return key() == rhs.key(); +} + AtomPack build_atom_pack_from_unitcell(const UnitCell& ucell, const double radius_lat0, const bool pbc) { const ExpandLayers layers = compute_expand_layers(ucell, radius_lat0, pbc); @@ -301,4 +424,141 @@ GridStorage build_grid_storage_from_atom_pack(const AtomPack& pack, const double return storage; } +std::vector build_neighbor_pairs_27(const AtomPack& pack, + const GridStorage& storage, + const double radius) +{ + validate_neighbor_search_input(storage, radius); + + std::vector pairs; + const double radius2 = radius * radius; + const int search_layer = std::max(1, static_cast(std::ceil(radius / storage.box_edge_length))); + + for (int center = 0; center < pack.size(); ++center) + { + if (!is_center_atom(pack, center)) + { + continue; + } + + int center_bx = 0; + int center_by = 0; + int center_bz = 0; + box_id_to_indices(storage, + storage.get_box_id_from_coord(pack.x[center], pack.y[center], pack.z[center]), + center_bx, + center_by, + center_bz); + + const int bx_begin = std::max(0, center_bx - search_layer); + const int bx_end = std::min(storage.box_nx - 1, center_bx + search_layer); + const int by_begin = std::max(0, center_by - search_layer); + const int by_end = std::min(storage.box_ny - 1, center_by + search_layer); + const int bz_begin = std::max(0, center_bz - search_layer); + const int bz_end = std::min(storage.box_nz - 1, center_bz + search_layer); + + for (int bx = bx_begin; bx <= bx_end; ++bx) + { + for (int by = by_begin; by <= by_end; ++by) + { + for (int bz = bz_begin; bz <= bz_end; ++bz) + { + const int box_id = storage.get_box_id(bx, by, bz); + const int begin = storage.box_offset[box_id]; + const int end = begin + storage.box_count[box_id]; + for (int offset = begin; offset < end; ++offset) + { + const int candidate = storage.atoms_in_box[offset]; + if (is_within_radius(pack, center, candidate, radius2)) + { + pairs.push_back(make_neighbor_pair(pack, center, candidate)); + } + } + } + } + } + } + + std::sort(pairs.begin(), pairs.end()); + return pairs; +} + +std::vector build_neighbor_pairs_14(const AtomPack& pack, + const GridStorage& storage, + const double radius) +{ + validate_neighbor_search_input(storage, radius); + + std::vector pairs; + const double radius2 = radius * radius; + const int search_layer = std::max(1, static_cast(std::ceil(radius / storage.box_edge_length))); + const std::map index_by_image = build_atom_image_index(pack); + + for (int center = 0; center < pack.size(); ++center) + { + if (!is_center_atom(pack, center)) + { + continue; + } + + int center_bx = 0; + int center_by = 0; + int center_bz = 0; + box_id_to_indices(storage, + storage.get_box_id_from_coord(pack.x[center], pack.y[center], pack.z[center]), + center_bx, + center_by, + center_bz); + + for (int dbx = -search_layer; dbx <= search_layer; ++dbx) + { + for (int dby = -search_layer; dby <= search_layer; ++dby) + { + for (int dbz = -search_layer; dbz <= search_layer; ++dbz) + { + if (!is_half_domain_offset(dbx, dby, dbz)) + { + continue; + } + + const int bx = center_bx + dbx; + const int by = center_by + dby; + const int bz = center_bz + dbz; + if (bx < 0 || bx >= storage.box_nx || by < 0 || by >= storage.box_ny || bz < 0 + || bz >= storage.box_nz) + { + continue; + } + + const int box_id = storage.get_box_id(bx, by, bz); + const int begin = storage.box_offset[box_id]; + const int end = begin + storage.box_count[box_id]; + for (int offset = begin; offset < end; ++offset) + { + const int candidate = storage.atoms_in_box[offset]; + if (!is_within_radius(pack, center, candidate, radius2)) + { + continue; + } + + const NeighborPair forward = make_neighbor_pair(pack, center, candidate); + const NeighborPair reverse = make_restored_reverse_pair(pack, forward, index_by_image); + if (forward.key() == reverse.key()) + { + continue; + } + + pairs.push_back(forward); + pairs.push_back(reverse); + } + } + } + } + } + + std::sort(pairs.begin(), pairs.end()); + pairs.erase(std::unique(pairs.begin(), pairs.end()), pairs.end()); + return pairs; +} + } // namespace ModuleNeighbor diff --git a/source/source_cell/module_neighbor/atom_pack.h b/source/source_cell/module_neighbor/atom_pack.h index c2db5daddbe..d3ccbc7484f 100644 --- a/source/source_cell/module_neighbor/atom_pack.h +++ b/source/source_cell/module_neighbor/atom_pack.h @@ -4,6 +4,7 @@ #include "mpi_domain.h" #include "source_cell/unitcell.h" +#include #include namespace ModuleNeighbor @@ -76,12 +77,41 @@ struct GridStorage int get_box_id_from_coord(const double x, const double y, const double z) const; }; +struct NeighborPair +{ + int center_type = 0; + int center_natom = 0; + int neighbor_type = 0; + int neighbor_natom = 0; + int cell_x = 0; + int cell_y = 0; + int cell_z = 0; + int center_index = -1; + int neighbor_index = -1; + + std::tuple key() const; + bool operator<(const NeighborPair& rhs) const; + bool operator==(const NeighborPair& rhs) const; +}; + // Build a flat atom pack from UnitCell. When pbc is true, the image layers follow // Grid::Check_Expand_Condition() so this helper remains comparable with Grid. AtomPack build_atom_pack_from_unitcell(const UnitCell& ucell, const double radius_lat0, const bool pbc); GridStorage build_grid_storage_from_atom_pack(const AtomPack& pack, const double box_edge_length); +std::vector build_neighbor_pairs_27(const AtomPack& pack, + const GridStorage& storage, + const double radius); + +// Build the same directed neighbor-pair result as build_neighbor_pairs_27(), but +// traverse only one half of the box-neighbor domain and restore the opposite +// direction explicitly. This keeps the result directly comparable with the +// current 27-direction baseline. +std::vector build_neighbor_pairs_14(const AtomPack& pack, + const GridStorage& storage, + const double radius); + } // namespace ModuleNeighbor #endif diff --git a/source/source_cell/module_neighbor/sltk_grid.cpp b/source/source_cell/module_neighbor/sltk_grid.cpp index b9b84b36657..cd784ea2b0c 100644 --- a/source/source_cell/module_neighbor/sltk_grid.cpp +++ b/source/source_cell/module_neighbor/sltk_grid.cpp @@ -33,7 +33,11 @@ Grid::~Grid() this->clear_atoms(); } -void Grid::init(std::ofstream& ofs_in, const UnitCell& ucell, const double radius_in, const bool boundary) +void Grid::init(std::ofstream& ofs_in, + const UnitCell& ucell, + const double radius_in, + const bool boundary, + const NeighborBuildMode build_mode) { ModuleBase::TITLE("Grid", "init"); ModuleBase::timer::start("Grid", "init"); @@ -49,7 +53,11 @@ void Grid::init(std::ofstream& ofs_in, const UnitCell& ucell, const double radiu ModuleBase::GlobalFunc::OUT(ofs_in, "Min number of cells", glayerX_minus, glayerY_minus, glayerZ_minus); this->setMemberVariables(ofs_in, ucell); - this->Construct_Adjacent(ucell); + this->Build_AtomPack_Search_Path(ucell); + if (build_mode == NeighborBuildMode::AtomPackAndLegacy) + { + this->Build_Legacy_Search_Path(ucell); + } ModuleBase::timer::end("Grid", "init"); } @@ -239,6 +247,24 @@ void Grid::setMemberVariables(std::ofstream& ofs_in, // output data to ofs this->box_nz = std::max(1, static_cast(std::floor((this->z_max - this->z_min) / box_edge_length)) + 1); ModuleBase::GlobalFunc::OUT(ofs_in, "Number of needed cells", box_nx, box_ny, box_nz); +} + +void Grid::Build_Legacy_Search_Path(const UnitCell& ucell) +{ + atoms_in_box.clear(); + all_adj_info.clear(); + + ModuleBase::Vector3 vec1(ucell.latvec.e11, ucell.latvec.e12, ucell.latvec.e13); + ModuleBase::Vector3 vec2(ucell.latvec.e21, ucell.latvec.e22, ucell.latvec.e23); + ModuleBase::Vector3 vec3(ucell.latvec.e31, ucell.latvec.e32, ucell.latvec.e33); + + const int ix_begin = this->pbc ? -glayerX_minus : 0; + const int ix_end = this->pbc ? glayerX : 1; + const int iy_begin = this->pbc ? -glayerY_minus : 0; + const int iy_end = this->pbc ? glayerY : 1; + const int iz_begin = this->pbc ? -glayerZ_minus : 0; + const int iz_end = this->pbc ? glayerZ : 1; + atoms_in_box.resize(this->box_nx); for (int i = 0; i < this->box_nx; i++) { @@ -273,12 +299,14 @@ void Grid::setMemberVariables(std::ofstream& ofs_in, // output data to ofs } } } - + this->all_adj_info.resize(ucell.ntype); for (int i = 0; i < ucell.ntype; i++) { this->all_adj_info[i].resize(ucell.atoms[i].na); } + + this->Construct_Adjacent(ucell); } void Grid::Construct_Adjacent(const UnitCell& ucell) @@ -358,3 +386,23 @@ void Grid::Construct_Adjacent_final(const FAtom& fatom1, all_adj_info[fatom1.type][fatom1.natom].push_back(fatom2); } } + +void Grid::Build_AtomPack_Search_Path(const UnitCell& ucell) +{ + atom_pack = ModuleNeighbor::build_atom_pack_from_unitcell(ucell, sradius, pbc); + grid_storage = ModuleNeighbor::build_grid_storage_from_atom_pack(atom_pack, box_edge_length); + neighbor_pairs_27 = ModuleNeighbor::build_neighbor_pairs_27(atom_pack, grid_storage, sradius); + + neighbor_pair_indices.clear(); + neighbor_pair_indices.resize(ucell.ntype); + for (int it = 0; it < ucell.ntype; ++it) + { + neighbor_pair_indices[it].resize(ucell.atoms[it].na); + } + + for (int pair_index = 0; pair_index < static_cast(neighbor_pairs_27.size()); ++pair_index) + { + const ModuleNeighbor::NeighborPair& pair = neighbor_pairs_27[pair_index]; + neighbor_pair_indices[pair.center_type][pair.center_natom].push_back(pair_index); + } +} diff --git a/source/source_cell/module_neighbor/sltk_grid.h b/source/source_cell/module_neighbor/sltk_grid.h index 8faf75fa2ea..a39b9238923 100644 --- a/source/source_cell/module_neighbor/sltk_grid.h +++ b/source/source_cell/module_neighbor/sltk_grid.h @@ -1,6 +1,7 @@ #ifndef GRID_H #define GRID_H +#include "atom_pack.h" #include "source_cell/unitcell.h" #include "sltk_atom.h" #include "sltk_util.h" @@ -15,6 +16,12 @@ typedef std::vector AtomMap; class Grid { public: + enum class NeighborBuildMode + { + AtomPackOnly, + AtomPackAndLegacy + }; + // Constructors and destructor // Grid is Global class,so init it with constant number Grid() : test_grid(0){}; @@ -23,7 +30,11 @@ class Grid Grid& operator=(Grid&&) = default; - void init(std::ofstream& ofs, const UnitCell& ucell, const double radius_in, const bool boundary = true); + void init(std::ofstream& ofs, + const UnitCell& ucell, + const double radius_in, + const bool boundary = true, + const NeighborBuildMode build_mode = NeighborBuildMode::AtomPackAndLegacy); // Data bool pbc=false; // When pbc is set to false, periodic boundary conditions are explicitly ignored. @@ -56,6 +67,15 @@ class Grid // Stores the adjacent information of atoms. [ntype][natom][adj list] std::vector >> all_adj_info; + + // Phase 2.1 flat search path. Grid_Driver::Find_atom() uses this + // integer-indexed route by default; the legacy FAtom* route above can still + // be built as a regression baseline. + ModuleNeighbor::AtomPack atom_pack; + ModuleNeighbor::GridStorage grid_storage; + std::vector neighbor_pairs_27; + std::vector>> neighbor_pair_indices; + void clear_atoms() { // we have to clear the all_adj_info @@ -63,6 +83,10 @@ class Grid all_adj_info.clear(); atoms_in_box.clear(); + atom_pack.clear(); + grid_storage.clear(); + neighbor_pairs_27.clear(); + neighbor_pair_indices.clear(); } void clear_adj_info() { @@ -102,6 +126,8 @@ class Grid void Construct_Adjacent(const UnitCell& ucell); void Construct_Adjacent_near_box(const FAtom& fatom); void Construct_Adjacent_final(const FAtom& fatom1, FAtom* fatom2); + void Build_AtomPack_Search_Path(const UnitCell& ucell); + void Build_Legacy_Search_Path(const UnitCell& ucell); void Check_Expand_Condition(const UnitCell& ucell); int glayerX=0; diff --git a/source/source_cell/module_neighbor/sltk_grid_driver.cpp b/source/source_cell/module_neighbor/sltk_grid_driver.cpp index 7986a37eebf..902ee0c42ff 100644 --- a/source/source_cell/module_neighbor/sltk_grid_driver.cpp +++ b/source/source_cell/module_neighbor/sltk_grid_driver.cpp @@ -28,11 +28,27 @@ void Grid_Driver::Find_atom(const UnitCell& ucell, AdjacentAtomInfo* adjs) const { ModuleBase::timer::start("Grid_Driver", "Find_atom"); + this->Find_atom_from_atom_pack(ucell, ntype, nnumber, adjs); + ModuleBase::timer::end("Grid_Driver", "Find_atom"); + return; +} + +void Grid_Driver::Find_atom_from_legacy(const UnitCell& ucell, + const int ntype, + const int nnumber, + AdjacentAtomInfo* adjs) const +{ // std::cout << "lenght in Find atom = " << atomlink[offset].fatom.getAdjacentSet()->getLength() << std::endl; // store result in member adj_info when parameter adjs is NULL AdjacentAtomInfo* local_adjs = adjs == nullptr ? &this->adj_info : adjs; local_adjs->clear(); + if (ntype < 0 || ntype >= static_cast(all_adj_info.size()) + || nnumber < 0 || nnumber >= static_cast(all_adj_info[ntype].size())) + { + throw std::runtime_error("Legacy Grid neighbor path is not built for this atom."); + } + const std::vector& all_atom = all_adj_info[ntype][nnumber]; for (const FAtom* atom: all_atom) @@ -51,9 +67,43 @@ void Grid_Driver::Find_atom(const UnitCell& ucell, local_adjs->natom.push_back(nnumber); local_adjs->box.push_back(ModuleBase::Vector3(0, 0, 0)); local_adjs->adjacent_tau.push_back(ModuleBase::Vector3(ucell.atoms[ntype].tau[nnumber].x, ucell.atoms[ntype].tau[nnumber].y, ucell.atoms[ntype].tau[nnumber].z)); - ModuleBase::timer::end("Grid_Driver", "Find_atom"); return; } + +void Grid_Driver::Find_atom_from_atom_pack(const UnitCell& ucell, + const int ntype, + const int nnumber, + AdjacentAtomInfo* adjs) const +{ + AdjacentAtomInfo* local_adjs = adjs == nullptr ? &this->adj_info : adjs; + local_adjs->clear(); + + const std::vector& pair_indices = neighbor_pair_indices[ntype][nnumber]; + for (const int pair_index: pair_indices) + { + const ModuleNeighbor::NeighborPair& pair = neighbor_pairs_27[pair_index]; + const int atom_index = pair.neighbor_index; + if (atom_index < 0 || atom_index >= atom_pack.size()) + { + throw std::runtime_error("AtomPack neighbor index is out of range."); + } + + local_adjs->ntype.push_back(pair.neighbor_type); + local_adjs->natom.push_back(pair.neighbor_natom); + local_adjs->box.push_back(ModuleBase::Vector3(pair.cell_x, pair.cell_y, pair.cell_z)); + local_adjs->adjacent_tau.push_back( + ModuleBase::Vector3(atom_pack.x[atom_index], atom_pack.y[atom_index], atom_pack.z[atom_index])); + local_adjs->adj_num++; + } + + local_adjs->ntype.push_back(ntype); + local_adjs->natom.push_back(nnumber); + local_adjs->box.push_back(ModuleBase::Vector3(0, 0, 0)); + local_adjs->adjacent_tau.push_back(ModuleBase::Vector3(ucell.atoms[ntype].tau[nnumber].x, + ucell.atoms[ntype].tau[nnumber].y, + ucell.atoms[ntype].tau[nnumber].z)); +} + void Grid_Driver::Find_atom(const UnitCell& ucell, const ModuleBase::Vector3& cartesian_posi, const int& ntype, diff --git a/source/source_cell/module_neighbor/sltk_grid_driver.h b/source/source_cell/module_neighbor/sltk_grid_driver.h index 9576617f59d..6942836d410 100644 --- a/source/source_cell/module_neighbor/sltk_grid_driver.h +++ b/source/source_cell/module_neighbor/sltk_grid_driver.h @@ -70,6 +70,22 @@ class Grid_Driver : public Grid const int nnumber, AdjacentAtomInfo* adjs = nullptr) const; + // Phase 2.1 transition helper: query neighbors through the AtomPack + + // GridStorage path built by Grid::init(). This is now the default + // Find_atom() route; keep the named helper for direct regression checks. + void Find_atom_from_atom_pack(const UnitCell& ucell, + const int ntype, + const int nnumber, + AdjacentAtomInfo* adjs = nullptr) const; + + // Phase 2.1 transition helper: keep the legacy all_adj_info route available + // for correctness comparisons while AtomPack replaces FAtom* in production + // queries. + void Find_atom_from_legacy(const UnitCell& ucell, + const int ntype, + const int nnumber, + AdjacentAtomInfo* adjs = nullptr) const; + // cartesian_posi and ucell is deprecated 20241204 zhanghaochong // this interface is deprecated, please use Find_atom above void Find_atom(const UnitCell& ucell, diff --git a/source/source_cell/module_neighbor/test/CMakeLists.txt b/source/source_cell/module_neighbor/test/CMakeLists.txt index c4d055ddcb0..4b7da8087cc 100644 --- a/source/source_cell/module_neighbor/test/CMakeLists.txt +++ b/source/source_cell/module_neighbor/test/CMakeLists.txt @@ -12,7 +12,7 @@ AddTest( AddTest( TARGET MODULE_CELL_NEIGHBOR_sltk_grid LIBS parameter ${math_libs} base device cell_info - SOURCES sltk_grid_test.cpp ../sltk_grid.cpp ../sltk_atom.cpp + SOURCES sltk_grid_test.cpp ../sltk_grid.cpp ../sltk_atom.cpp ../atom_pack.cpp ../../../source_io/module_output/output.cpp ) @@ -20,7 +20,7 @@ AddTest( TARGET MODULE_CELL_NEIGHBOR_sltk_atom_arrange LIBS parameter ${math_libs} base device cell_info SOURCES sltk_atom_arrange_test.cpp ../sltk_atom_arrange.cpp ../sltk_grid_driver.cpp ../sltk_grid.cpp - ../sltk_atom.cpp + ../sltk_atom.cpp ../atom_pack.cpp ../../../source_io/module_output/output.cpp ) diff --git a/source/source_cell/module_neighbor/test/atom_pack_test.cpp b/source/source_cell/module_neighbor/test/atom_pack_test.cpp index a3471a572da..0595324b825 100644 --- a/source/source_cell/module_neighbor/test/atom_pack_test.cpp +++ b/source/source_cell/module_neighbor/test/atom_pack_test.cpp @@ -66,6 +66,53 @@ std::vector collect_grid_box_counts(const Grid& grid) } return counts; } + +std::vector collect_legacy_neighbor_pairs(const Grid& grid) +{ + std::vector pairs; + for (int center_type = 0; center_type < static_cast(grid.all_adj_info.size()); ++center_type) + { + const auto& type_neighbors = grid.all_adj_info[center_type]; + for (int center_natom = 0; center_natom < static_cast(type_neighbors.size()); ++center_natom) + { + const auto& atom_neighbors = type_neighbors[center_natom]; + for (const FAtom* atom: atom_neighbors) + { + pairs.push_back(ModuleNeighbor::NeighborPair{center_type, + center_natom, + atom->type, + atom->natom, + atom->cell_x, + atom->cell_y, + atom->cell_z}); + } + } + } + std::sort(pairs.begin(), pairs.end()); + return pairs; +} + +void expect_neighbor_pair_indices_match_pack(const std::vector& pairs, + const ModuleNeighbor::AtomPack& pack) +{ + for (const ModuleNeighbor::NeighborPair& pair: pairs) + { + ASSERT_GE(pair.center_index, 0); + ASSERT_LT(pair.center_index, pack.size()); + ASSERT_GE(pair.neighbor_index, 0); + ASSERT_LT(pair.neighbor_index, pack.size()); + EXPECT_EQ(pack.type[pair.center_index], pair.center_type); + EXPECT_EQ(pack.natom[pair.center_index], pair.center_natom); + EXPECT_EQ(pack.cell_x[pair.center_index], 0); + EXPECT_EQ(pack.cell_y[pair.center_index], 0); + EXPECT_EQ(pack.cell_z[pair.center_index], 0); + EXPECT_EQ(pack.type[pair.neighbor_index], pair.neighbor_type); + EXPECT_EQ(pack.natom[pair.neighbor_index], pair.neighbor_natom); + EXPECT_EQ(pack.cell_x[pair.neighbor_index], pair.cell_x); + EXPECT_EQ(pack.cell_y[pair.neighbor_index], pair.cell_y); + EXPECT_EQ(pack.cell_z[pair.neighbor_index], pair.cell_z); + } +} } // namespace TEST(AtomPackTest, BuildsNonPeriodicPackWithoutImages) @@ -215,3 +262,119 @@ TEST(AtomPackTest, FlatGridStorageMatchesLegacyGridBoxCountsWithoutPbc) remove("atom_pack_grid_compare.out"); delete ucell; } + +TEST(AtomPackTest, FlatGridSearchMatchesLegacyGridPairsWithoutPbc) +{ + UcellTestPrepare utp = UcellTestLib["Si"]; + UnitCell* ucell = utp.SetUcellInfo(); + unitcell::check_dtau(ucell->atoms, ucell->ntype, ucell->lat0, ucell->latvec); + + const double radius = 0.5; + const ModuleNeighbor::AtomPack pack = ModuleNeighbor::build_atom_pack_from_unitcell(*ucell, radius, false); + const ModuleNeighbor::GridStorage storage = ModuleNeighbor::build_grid_storage_from_atom_pack(pack, radius + 0.1); + const std::vector flat_pairs + = ModuleNeighbor::build_neighbor_pairs_27(pack, storage, radius); + + std::ofstream ofs("atom_pack_neighbor_compare_non_pbc.out"); + Grid grid(0); + grid.init(ofs, *ucell, radius, false); + ofs.close(); + const std::vector legacy_pairs = collect_legacy_neighbor_pairs(grid); + + EXPECT_EQ(flat_pairs, legacy_pairs); + EXPECT_EQ(std::adjacent_find(flat_pairs.begin(), flat_pairs.end()), flat_pairs.end()); + + remove("atom_pack_neighbor_compare_non_pbc.out"); + delete ucell; +} + +TEST(AtomPackTest, FlatGridSearchMatchesLegacyGridPairsWithPbc) +{ + UcellTestPrepare utp = UcellTestLib["Si"]; + UnitCell* ucell = utp.SetUcellInfo(); + unitcell::check_dtau(ucell->atoms, ucell->ntype, ucell->lat0, ucell->latvec); + + const double radius = 0.5; + const ModuleNeighbor::AtomPack pack = ModuleNeighbor::build_atom_pack_from_unitcell(*ucell, radius, true); + const ModuleNeighbor::GridStorage storage = ModuleNeighbor::build_grid_storage_from_atom_pack(pack, radius + 0.1); + const std::vector flat_pairs + = ModuleNeighbor::build_neighbor_pairs_27(pack, storage, radius); + + std::ofstream ofs("atom_pack_neighbor_compare_pbc.out"); + Grid grid(0); + grid.init(ofs, *ucell, radius, true); + ofs.close(); + const std::vector legacy_pairs = collect_legacy_neighbor_pairs(grid); + + EXPECT_EQ(flat_pairs, legacy_pairs); + EXPECT_EQ(std::adjacent_find(flat_pairs.begin(), flat_pairs.end()), flat_pairs.end()); + + remove("atom_pack_neighbor_compare_pbc.out"); + delete ucell; +} + +TEST(AtomPackTest, HalfDomainSearchRestoresFullPairsWithoutPbc) +{ + UcellTestPrepare utp = UcellTestLib["Si"]; + UnitCell* ucell = utp.SetUcellInfo(); + unitcell::check_dtau(ucell->atoms, ucell->ntype, ucell->lat0, ucell->latvec); + + const double radius = 0.5; + const ModuleNeighbor::AtomPack pack = ModuleNeighbor::build_atom_pack_from_unitcell(*ucell, radius, false); + const ModuleNeighbor::GridStorage storage = ModuleNeighbor::build_grid_storage_from_atom_pack(pack, radius + 0.1); + const std::vector full_pairs + = ModuleNeighbor::build_neighbor_pairs_27(pack, storage, radius); + const std::vector restored_pairs + = ModuleNeighbor::build_neighbor_pairs_14(pack, storage, radius); + + EXPECT_EQ(restored_pairs, full_pairs); + EXPECT_EQ(std::adjacent_find(restored_pairs.begin(), restored_pairs.end()), restored_pairs.end()); + expect_neighbor_pair_indices_match_pack(restored_pairs, pack); + + delete ucell; +} + +TEST(AtomPackTest, HalfDomainSearchRestoresFullPairsWithPbc) +{ + UcellTestPrepare utp = UcellTestLib["Si"]; + UnitCell* ucell = utp.SetUcellInfo(); + unitcell::check_dtau(ucell->atoms, ucell->ntype, ucell->lat0, ucell->latvec); + + const double radius = 0.5; + const ModuleNeighbor::AtomPack pack = ModuleNeighbor::build_atom_pack_from_unitcell(*ucell, radius, true); + const ModuleNeighbor::GridStorage storage = ModuleNeighbor::build_grid_storage_from_atom_pack(pack, radius + 0.1); + const std::vector full_pairs + = ModuleNeighbor::build_neighbor_pairs_27(pack, storage, radius); + const std::vector restored_pairs + = ModuleNeighbor::build_neighbor_pairs_14(pack, storage, radius); + + EXPECT_EQ(restored_pairs, full_pairs); + EXPECT_EQ(std::adjacent_find(restored_pairs.begin(), restored_pairs.end()), restored_pairs.end()); + expect_neighbor_pair_indices_match_pack(restored_pairs, pack); + + delete ucell; +} + +TEST(AtomPackTest, FlatGridSearchHandlesSingleAtomWithoutSelfNeighbor) +{ + ModuleNeighbor::AtomPack pack; + pack.append_atom(0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, false); + + const ModuleNeighbor::GridStorage storage = ModuleNeighbor::build_grid_storage_from_atom_pack(pack, 0.6); + const std::vector pairs + = ModuleNeighbor::build_neighbor_pairs_27(pack, storage, 0.5); + + EXPECT_TRUE(pairs.empty()); +} + +TEST(AtomPackTest, HalfDomainSearchHandlesSingleAtomWithoutSelfNeighbor) +{ + ModuleNeighbor::AtomPack pack; + pack.append_atom(0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, false); + + const ModuleNeighbor::GridStorage storage = ModuleNeighbor::build_grid_storage_from_atom_pack(pack, 0.6); + const std::vector pairs + = ModuleNeighbor::build_neighbor_pairs_14(pack, storage, 0.5); + + EXPECT_TRUE(pairs.empty()); +} diff --git a/source/source_cell/module_neighbor/test/sltk_atom_arrange_test.cpp b/source/source_cell/module_neighbor/test/sltk_atom_arrange_test.cpp index 329554d62c3..632a7db6ea9 100644 --- a/source/source_cell/module_neighbor/test/sltk_atom_arrange_test.cpp +++ b/source/source_cell/module_neighbor/test/sltk_atom_arrange_test.cpp @@ -3,8 +3,11 @@ #define private public #include "source_io/module_parameter/parameter.h" #undef private +#include #include #include +#include +#include #include "gmock/gmock.h" #include "gtest/gtest.h" @@ -35,6 +38,140 @@ Magnetism::~Magnetism() delete[] this->start_mag; } +namespace +{ +using AdjacentKey = std::tuple; + +std::vector collect_adjacent_keys(const AdjacentAtomInfo& adjs) +{ + std::vector keys; + for (int i = 0; i < static_cast(adjs.ntype.size()); ++i) + { + keys.push_back(AdjacentKey(adjs.ntype[i], + adjs.natom[i], + adjs.box[i].x, + adjs.box[i].y, + adjs.box[i].z, + adjs.adjacent_tau[i].x, + adjs.adjacent_tau[i].y, + adjs.adjacent_tau[i].z)); + } + std::sort(keys.begin(), keys.end()); + return keys; +} + +void expect_self_is_last(const AdjacentAtomInfo& adjs, + const UnitCell& ucell, + const int ntype, + const int nnumber) +{ + ASSERT_FALSE(adjs.ntype.empty()); + const int last = static_cast(adjs.ntype.size()) - 1; + EXPECT_EQ(adjs.ntype[last], ntype); + EXPECT_EQ(adjs.natom[last], nnumber); + EXPECT_EQ(adjs.box[last].x, 0); + EXPECT_EQ(adjs.box[last].y, 0); + EXPECT_EQ(adjs.box[last].z, 0); + EXPECT_DOUBLE_EQ(adjs.adjacent_tau[last].x, ucell.atoms[ntype].tau[nnumber].x); + EXPECT_DOUBLE_EQ(adjs.adjacent_tau[last].y, ucell.atoms[ntype].tau[nnumber].y); + EXPECT_DOUBLE_EQ(adjs.adjacent_tau[last].z, ucell.atoms[ntype].tau[nnumber].z); +} + +void expect_grid_driver_default_path_matches_legacy(const UnitCell& ucell, const double radius, const bool pbc) +{ + std::ofstream ofs("grid_driver_atom_pack_compare.out"); + Grid_Driver grid_d(PARAM.input.test_deconstructor, PARAM.input.test_grid); + grid_d.init(ofs, ucell, radius, pbc); + ofs.close(); + + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + SCOPED_TRACE("pbc=" + std::to_string(pbc) + ", radius=" + std::to_string(radius) + ", type=" + + std::to_string(it) + ", atom=" + std::to_string(ia)); + + AdjacentAtomInfo default_adjs; + AdjacentAtomInfo legacy_adjs; + AdjacentAtomInfo atom_pack_adjs; + AdjacentAtomInfo deprecated_adjs; + grid_d.Find_atom(ucell, it, ia, &default_adjs); + grid_d.Find_atom_from_legacy(ucell, it, ia, &legacy_adjs); + grid_d.Find_atom_from_atom_pack(ucell, it, ia, &atom_pack_adjs); + grid_d.Find_atom(ucell, ucell.atoms[it].tau[ia], it, ia, &deprecated_adjs); + + EXPECT_EQ(default_adjs.adj_num, legacy_adjs.adj_num); + EXPECT_EQ(default_adjs.ntype.size(), legacy_adjs.ntype.size()); + EXPECT_EQ(default_adjs.natom.size(), legacy_adjs.natom.size()); + EXPECT_EQ(default_adjs.box.size(), legacy_adjs.box.size()); + EXPECT_EQ(default_adjs.adjacent_tau.size(), legacy_adjs.adjacent_tau.size()); + EXPECT_EQ(collect_adjacent_keys(default_adjs), collect_adjacent_keys(legacy_adjs)); + EXPECT_EQ(atom_pack_adjs.adj_num, legacy_adjs.adj_num); + EXPECT_EQ(atom_pack_adjs.ntype.size(), legacy_adjs.ntype.size()); + EXPECT_EQ(atom_pack_adjs.natom.size(), legacy_adjs.natom.size()); + EXPECT_EQ(atom_pack_adjs.box.size(), legacy_adjs.box.size()); + EXPECT_EQ(atom_pack_adjs.adjacent_tau.size(), legacy_adjs.adjacent_tau.size()); + EXPECT_EQ(collect_adjacent_keys(atom_pack_adjs), collect_adjacent_keys(legacy_adjs)); + EXPECT_EQ(deprecated_adjs.adj_num, legacy_adjs.adj_num); + EXPECT_EQ(deprecated_adjs.ntype.size(), legacy_adjs.ntype.size()); + EXPECT_EQ(deprecated_adjs.natom.size(), legacy_adjs.natom.size()); + EXPECT_EQ(deprecated_adjs.box.size(), legacy_adjs.box.size()); + EXPECT_EQ(deprecated_adjs.adjacent_tau.size(), legacy_adjs.adjacent_tau.size()); + EXPECT_EQ(collect_adjacent_keys(deprecated_adjs), collect_adjacent_keys(legacy_adjs)); + expect_self_is_last(default_adjs, ucell, it, ia); + expect_self_is_last(legacy_adjs, ucell, it, ia); + expect_self_is_last(atom_pack_adjs, ucell, it, ia); + expect_self_is_last(deprecated_adjs, ucell, it, ia); + } + } + + remove("grid_driver_atom_pack_compare.out"); +} + +void expect_grid_driver_atom_pack_only_path(const UnitCell& ucell, const double radius, const bool pbc) +{ + std::ofstream ofs("grid_driver_atom_pack_only.out"); + Grid_Driver grid_d(PARAM.input.test_deconstructor, PARAM.input.test_grid); + grid_d.init(ofs, ucell, radius, pbc, Grid::NeighborBuildMode::AtomPackOnly); + ofs.close(); + + EXPECT_TRUE(grid_d.atoms_in_box.empty()); + EXPECT_TRUE(grid_d.all_adj_info.empty()); + EXPECT_FALSE(grid_d.atom_pack.empty()); + EXPECT_FALSE(grid_d.neighbor_pair_indices.empty()); + + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + SCOPED_TRACE("AtomPackOnly pbc=" + std::to_string(pbc) + ", radius=" + std::to_string(radius) + + ", type=" + std::to_string(it) + ", atom=" + std::to_string(ia)); + + AdjacentAtomInfo default_adjs; + AdjacentAtomInfo atom_pack_adjs; + AdjacentAtomInfo deprecated_adjs; + grid_d.Find_atom(ucell, it, ia, &default_adjs); + grid_d.Find_atom_from_atom_pack(ucell, it, ia, &atom_pack_adjs); + grid_d.Find_atom(ucell, ucell.atoms[it].tau[ia], it, ia, &deprecated_adjs); + + EXPECT_EQ(default_adjs.adj_num, atom_pack_adjs.adj_num); + EXPECT_EQ(default_adjs.ntype.size(), atom_pack_adjs.ntype.size()); + EXPECT_EQ(default_adjs.natom.size(), atom_pack_adjs.natom.size()); + EXPECT_EQ(default_adjs.box.size(), atom_pack_adjs.box.size()); + EXPECT_EQ(default_adjs.adjacent_tau.size(), atom_pack_adjs.adjacent_tau.size()); + EXPECT_EQ(collect_adjacent_keys(default_adjs), collect_adjacent_keys(atom_pack_adjs)); + EXPECT_EQ(collect_adjacent_keys(deprecated_adjs), collect_adjacent_keys(atom_pack_adjs)); + expect_self_is_last(default_adjs, ucell, it, ia); + expect_self_is_last(atom_pack_adjs, ucell, it, ia); + expect_self_is_last(deprecated_adjs, ucell, it, ia); + EXPECT_THROW(grid_d.Find_atom_from_legacy(ucell, it, ia, nullptr), std::runtime_error); + } + } + + remove("grid_driver_atom_pack_only.out"); +} +} // namespace + /************************************************ * unit test of atom_arrange ***********************************************/ @@ -153,3 +290,19 @@ TEST_F(SltkAtomArrangeTest, Filteradjs) filter_adjs(is_adjs, adjs); EXPECT_EQ(adjs.adj_num, 0); } + +TEST_F(SltkAtomArrangeTest, GridDriverAtomPackPathMatchesLegacyPath) +{ + unitcell::check_dtau(ucell->atoms, ucell->ntype, ucell->lat0, ucell->latvec); + + expect_grid_driver_default_path_matches_legacy(*ucell, radius, true); + expect_grid_driver_default_path_matches_legacy(*ucell, 0.5, false); +} + +TEST_F(SltkAtomArrangeTest, GridDriverAtomPackOnlyPathWorksWithoutLegacyBuild) +{ + unitcell::check_dtau(ucell->atoms, ucell->ntype, ucell->lat0, ucell->latvec); + + expect_grid_driver_atom_pack_only_path(*ucell, radius, true); + expect_grid_driver_atom_pack_only_path(*ucell, 0.5, false); +} diff --git a/source/source_cell/module_neighbor/test/sltk_grid_test.cpp b/source/source_cell/module_neighbor/test/sltk_grid_test.cpp index 65843000981..6edc6078b5c 100644 --- a/source/source_cell/module_neighbor/test/sltk_grid_test.cpp +++ b/source/source_cell/module_neighbor/test/sltk_grid_test.cpp @@ -110,6 +110,31 @@ std::vector collect_neighbor_pair_keys(const Grid& grid) std::sort(keys.begin(), keys.end()); return keys; } + +std::vector collect_legacy_neighbor_pairs(const Grid& grid) +{ + std::vector pairs; + for (int center_type = 0; center_type < static_cast(grid.all_adj_info.size()); ++center_type) + { + const auto& type_neighbors = grid.all_adj_info[center_type]; + for (int center_natom = 0; center_natom < static_cast(type_neighbors.size()); ++center_natom) + { + const auto& atom_neighbors = type_neighbors[center_natom]; + for (const FAtom* atom: atom_neighbors) + { + pairs.push_back(ModuleNeighbor::NeighborPair{center_type, + center_natom, + atom->type, + atom->natom, + atom->cell_x, + atom->cell_y, + atom->cell_z}); + } + } + } + std::sort(pairs.begin(), pairs.end()); + return pairs; +} } // namespace TEST_F(SltkGridTest, Init) @@ -186,6 +211,26 @@ TEST_F(SltkGridTest, OpenMPThreadCountKeepsNeighborSet) remove("test_many.out"); } +TEST_F(SltkGridTest, AtomPackNeighborPairsMatchLegacyGridAfterInit) +{ + unitcell::check_dtau(ucell->atoms, ucell->ntype, ucell->lat0, ucell->latvec); + + std::ofstream ofs_pbc("test_pair_pbc.out"); + Grid grid_pbc(PARAM.input.test_grid); + grid_pbc.init(ofs_pbc, *ucell, radius, true); + ofs_pbc.close(); + EXPECT_EQ(grid_pbc.neighbor_pairs_27, collect_legacy_neighbor_pairs(grid_pbc)); + + std::ofstream ofs_non_pbc("test_pair_non_pbc.out"); + Grid grid_non_pbc(PARAM.input.test_grid); + grid_non_pbc.init(ofs_non_pbc, *ucell, 0.5, false); + ofs_non_pbc.close(); + EXPECT_EQ(grid_non_pbc.neighbor_pairs_27, collect_legacy_neighbor_pairs(grid_non_pbc)); + + remove("test_pair_pbc.out"); + remove("test_pair_non_pbc.out"); +} + /* // This test cannot pass because setAtomLinkArray() is unsuccessful // if expand_flag is false From c0c8d64e3ee5468aceabbd19045c55ec085cc8f2 Mon Sep 17 00:00:00 2001 From: Yanglijf Date: Sat, 6 Jun 2026 20:16:49 +0800 Subject: [PATCH 2/2] Fix NeighborPair initialization for C++11 builds --- .../source_cell/module_neighbor/atom_pack.cpp | 46 +++++++++++-------- .../source_cell/module_neighbor/atom_pack.h | 6 +++ .../source_cell/module_neighbor/sltk_grid.cpp | 6 +++ .../module_neighbor/sltk_grid_driver.cpp | 10 ++++ .../module_neighbor/test/atom_pack_test.cpp | 16 ++++--- .../module_neighbor/test/sltk_grid_test.cpp | 16 ++++--- 6 files changed, 68 insertions(+), 32 deletions(-) diff --git a/source/source_cell/module_neighbor/atom_pack.cpp b/source/source_cell/module_neighbor/atom_pack.cpp index 93874c8dffb..a20f3bb63e4 100644 --- a/source/source_cell/module_neighbor/atom_pack.cpp +++ b/source/source_cell/module_neighbor/atom_pack.cpp @@ -140,21 +140,29 @@ void validate_neighbor_search_input(const GridStorage& storage, const double rad NeighborPair make_neighbor_pair(const AtomPack& pack, const int center, const int candidate) { - return NeighborPair{pack.type[center], - pack.natom[center], - pack.type[candidate], - pack.natom[candidate], - pack.cell_x[candidate], - pack.cell_y[candidate], - pack.cell_z[candidate], - center, - candidate}; + // Keep explicit assignment for C++11 CI builds: NeighborPair has default + // member initializers and member functions, so brace-list aggregate + // initialization is not accepted by all toolchain variants. + NeighborPair pair; + pair.center_type = pack.type[center]; + pair.center_natom = pack.natom[center]; + pair.neighbor_type = pack.type[candidate]; + pair.neighbor_natom = pack.natom[candidate]; + pair.cell_x = pack.cell_x[candidate]; + pair.cell_y = pack.cell_y[candidate]; + pair.cell_z = pack.cell_z[candidate]; + pair.center_index = center; + pair.neighbor_index = candidate; + return pair; } NeighborPair make_restored_reverse_pair(const AtomPack& pack, const NeighborPair& pair, const std::map& index_by_image) { + // Half-domain search visits only one direction. The restored pair uses the + // neighbor's origin-cell image as the new center and the opposite image shift + // for the original center atom. const int reverse_center = find_atom_image_index(index_by_image, pair.neighbor_type, pair.neighbor_natom, 0, 0, 0); const int reverse_neighbor = find_atom_image_index(index_by_image, pair.center_type, @@ -163,15 +171,17 @@ NeighborPair make_restored_reverse_pair(const AtomPack& pack, -pair.cell_y, -pair.cell_z); - return NeighborPair{pair.neighbor_type, - pair.neighbor_natom, - pair.center_type, - pair.center_natom, - -pair.cell_x, - -pair.cell_y, - -pair.cell_z, - reverse_center, - reverse_neighbor}; + NeighborPair reverse; + reverse.center_type = pair.neighbor_type; + reverse.center_natom = pair.neighbor_natom; + reverse.neighbor_type = pair.center_type; + reverse.neighbor_natom = pair.center_natom; + reverse.cell_x = -pair.cell_x; + reverse.cell_y = -pair.cell_y; + reverse.cell_z = -pair.cell_z; + reverse.center_index = reverse_center; + reverse.neighbor_index = reverse_neighbor; + return reverse; } bool is_within_radius(const AtomPack& pack, const int center, const int candidate, const double radius2) diff --git a/source/source_cell/module_neighbor/atom_pack.h b/source/source_cell/module_neighbor/atom_pack.h index d3ccbc7484f..213e4889412 100644 --- a/source/source_cell/module_neighbor/atom_pack.h +++ b/source/source_cell/module_neighbor/atom_pack.h @@ -79,6 +79,8 @@ struct GridStorage struct NeighborPair { + // Complete directed pair key used for correctness comparisons and stable + // sorting. cell_* is the periodic image shift of the neighbor atom. int center_type = 0; int center_natom = 0; int neighbor_type = 0; @@ -86,6 +88,10 @@ struct NeighborPair int cell_x = 0; int cell_y = 0; int cell_z = 0; + + // Internal AtomPack indices used to recover coordinates when converting the + // pair list back to AdjacentAtomInfo. They are deliberately ignored by + // operator< and operator== so tests compare the physical neighbor relation. int center_index = -1; int neighbor_index = -1; diff --git a/source/source_cell/module_neighbor/sltk_grid.cpp b/source/source_cell/module_neighbor/sltk_grid.cpp index cd784ea2b0c..50a74c23d3e 100644 --- a/source/source_cell/module_neighbor/sltk_grid.cpp +++ b/source/source_cell/module_neighbor/sltk_grid.cpp @@ -389,10 +389,16 @@ void Grid::Construct_Adjacent_final(const FAtom& fatom1, void Grid::Build_AtomPack_Search_Path(const UnitCell& ucell) { + // Build the Phase 2.1 integer-indexed path in the same init() pass as the + // legacy boxes. This keeps both paths comparable while Grid_Driver moves to + // AtomPack as the default query backend. atom_pack = ModuleNeighbor::build_atom_pack_from_unitcell(ucell, sradius, pbc); grid_storage = ModuleNeighbor::build_grid_storage_from_atom_pack(atom_pack, box_edge_length); neighbor_pairs_27 = ModuleNeighbor::build_neighbor_pairs_27(atom_pack, grid_storage, sradius); + // Convert the global pair list to per-center lookup indices so Find_atom() + // can keep its existing single-atom query interface without scanning all + // neighbor pairs each time. neighbor_pair_indices.clear(); neighbor_pair_indices.resize(ucell.ntype); for (int it = 0; it < ucell.ntype; ++it) diff --git a/source/source_cell/module_neighbor/sltk_grid_driver.cpp b/source/source_cell/module_neighbor/sltk_grid_driver.cpp index 902ee0c42ff..fcaffae67d8 100644 --- a/source/source_cell/module_neighbor/sltk_grid_driver.cpp +++ b/source/source_cell/module_neighbor/sltk_grid_driver.cpp @@ -28,6 +28,8 @@ void Grid_Driver::Find_atom(const UnitCell& ucell, AdjacentAtomInfo* adjs) const { ModuleBase::timer::start("Grid_Driver", "Find_atom"); + // The public interface now uses the Phase 2.1 AtomPack path by default. + // Find_atom_from_legacy() is kept only as a regression and fallback route. this->Find_atom_from_atom_pack(ucell, ntype, nnumber, adjs); ModuleBase::timer::end("Grid_Driver", "Find_atom"); return; @@ -77,6 +79,11 @@ void Grid_Driver::Find_atom_from_atom_pack(const UnitCell& ucell, { AdjacentAtomInfo* local_adjs = adjs == nullptr ? &this->adj_info : adjs; local_adjs->clear(); + if (ntype < 0 || ntype >= static_cast(neighbor_pair_indices.size()) + || nnumber < 0 || nnumber >= static_cast(neighbor_pair_indices[ntype].size())) + { + throw std::runtime_error("AtomPack neighbor path is not built for this atom."); + } const std::vector& pair_indices = neighbor_pair_indices[ntype][nnumber]; for (const int pair_index: pair_indices) @@ -96,6 +103,9 @@ void Grid_Driver::Find_atom_from_atom_pack(const UnitCell& ucell, local_adjs->adj_num++; } + // Keep the ABACUS compatibility rule from the legacy path: the center atom + // itself is appended after all real neighbors, and adj_num counts only the + // real neighbor entries. local_adjs->ntype.push_back(ntype); local_adjs->natom.push_back(nnumber); local_adjs->box.push_back(ModuleBase::Vector3(0, 0, 0)); diff --git a/source/source_cell/module_neighbor/test/atom_pack_test.cpp b/source/source_cell/module_neighbor/test/atom_pack_test.cpp index 0595324b825..2221cf92794 100644 --- a/source/source_cell/module_neighbor/test/atom_pack_test.cpp +++ b/source/source_cell/module_neighbor/test/atom_pack_test.cpp @@ -78,13 +78,15 @@ std::vector collect_legacy_neighbor_pairs(const Gr const auto& atom_neighbors = type_neighbors[center_natom]; for (const FAtom* atom: atom_neighbors) { - pairs.push_back(ModuleNeighbor::NeighborPair{center_type, - center_natom, - atom->type, - atom->natom, - atom->cell_x, - atom->cell_y, - atom->cell_z}); + ModuleNeighbor::NeighborPair pair; + pair.center_type = center_type; + pair.center_natom = center_natom; + pair.neighbor_type = atom->type; + pair.neighbor_natom = atom->natom; + pair.cell_x = atom->cell_x; + pair.cell_y = atom->cell_y; + pair.cell_z = atom->cell_z; + pairs.push_back(pair); } } } diff --git a/source/source_cell/module_neighbor/test/sltk_grid_test.cpp b/source/source_cell/module_neighbor/test/sltk_grid_test.cpp index 6edc6078b5c..884e79c2513 100644 --- a/source/source_cell/module_neighbor/test/sltk_grid_test.cpp +++ b/source/source_cell/module_neighbor/test/sltk_grid_test.cpp @@ -122,13 +122,15 @@ std::vector collect_legacy_neighbor_pairs(const Gr const auto& atom_neighbors = type_neighbors[center_natom]; for (const FAtom* atom: atom_neighbors) { - pairs.push_back(ModuleNeighbor::NeighborPair{center_type, - center_natom, - atom->type, - atom->natom, - atom->cell_x, - atom->cell_y, - atom->cell_z}); + ModuleNeighbor::NeighborPair pair; + pair.center_type = center_type; + pair.center_natom = center_natom; + pair.neighbor_type = atom->type; + pair.neighbor_natom = atom->natom; + pair.cell_x = atom->cell_x; + pair.cell_y = atom->cell_y; + pair.cell_z = atom->cell_z; + pairs.push_back(pair); } } }