From 528ac1869afee3fd0efc65fdafc8c13d49fa0b02 Mon Sep 17 00:00:00 2001 From: hhb2006 <1035516541@qq.com> Date: Sat, 30 May 2026 23:04:27 +0800 Subject: [PATCH 01/15] Optimize neighbor grid search in source/source_cell/module_neighbor/sltk_grid.cpp --- .../source_cell/module_neighbor/sltk_grid.cpp | 111 +++++++++++------- 1 file changed, 69 insertions(+), 42 deletions(-) diff --git a/source/source_cell/module_neighbor/sltk_grid.cpp b/source/source_cell/module_neighbor/sltk_grid.cpp index 9aad1a64e45..b9b84b36657 100644 --- a/source/source_cell/module_neighbor/sltk_grid.cpp +++ b/source/source_cell/module_neighbor/sltk_grid.cpp @@ -5,6 +5,25 @@ #include "source_base/memory.h" #include "source_base/timer.h" +#include +#include + +namespace +{ +int clamp_box_index(const int index, const int size) +{ + if (index < 0) + { + return 0; + } + if (index >= size) + { + return size - 1; + } + return index; +} +} + Grid::Grid(const int& test_grid_in) : test_grid(test_grid_in) { } @@ -40,6 +59,12 @@ void Grid::Check_Expand_Condition(const UnitCell& ucell) if (!pbc) { + glayerX = 1; + glayerY = 1; + glayerZ = 1; + glayerX_minus = 0; + glayerY_minus = 0; + glayerZ_minus = 0; return; } @@ -170,12 +195,19 @@ void Grid::setMemberVariables(std::ofstream& ofs_in, // output data to ofs 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; + // calculate min & max value - for (int ix = -glayerX_minus; ix < glayerX; ix++) + for (int ix = ix_begin; ix < ix_end; ix++) { - for (int iy = -glayerY_minus; iy < glayerY; iy++) + for (int iy = iy_begin; iy < iy_end; iy++) { - for (int iz = -glayerZ_minus; iz < glayerZ; iz++) + for (int iz = iz_begin; iz < iz_end; iz++) { for (int i = 0; i < ucell.ntype; i++) { @@ -202,25 +234,9 @@ void Grid::setMemberVariables(std::ofstream& ofs_in, // output data to ofs this->box_edge_length = sradius + 0.1; // To avoid edge cases, the size of the box is slightly increased. -/* warning box algorithm - this->box_nx = std::ceil((this->x_max - this->x_min) / box_edge_length) + 1; - this->box_ny = std::ceil((this->y_max - this->y_min) / box_edge_length) + 1; - this->box_nz = std::ceil((this->z_max - this->z_min) / box_edge_length) + 1; - ModuleBase::GlobalFunc::OUT(ofs_in, "BoxNumber", box_nx, box_ny, box_nz); - - atoms_in_box.resize(this->box_nx); - for (int i = 0; i < this->box_nx; i++) - { - atoms_in_box[i].resize(this->box_ny); - for (int j = 0; j < this->box_ny; j++) - { - atoms_in_box[i][j].resize(this->box_nz); - } - } - */ - this->box_nx = glayerX + glayerX_minus; - this->box_ny = glayerY + glayerY_minus; - this->box_nz = glayerZ + glayerZ_minus; + this->box_nx = std::max(1, static_cast(std::floor((this->x_max - this->x_min) / box_edge_length)) + 1); + this->box_ny = std::max(1, static_cast(std::floor((this->y_max - this->y_min) / box_edge_length)) + 1); + 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); atoms_in_box.resize(this->box_nx); @@ -232,11 +248,11 @@ void Grid::setMemberVariables(std::ofstream& ofs_in, // output data to ofs atoms_in_box[i][j].resize(this->box_nz); } } - for (int ix = -glayerX_minus; ix < glayerX; ix++) + for (int ix = ix_begin; ix < ix_end; ix++) { - for (int iy = -glayerY_minus; iy < glayerY; iy++) + for (int iy = iy_begin; iy < iy_end; iy++) { - for (int iz = -glayerZ_minus; iz < glayerZ; iz++) + for (int iz = iz_begin; iz < iz_end; iz++) { for (int i = 0; i < ucell.ntype; i++) { @@ -247,10 +263,10 @@ void Grid::setMemberVariables(std::ofstream& ofs_in, // output data to ofs double z = ucell.atoms[i].tau[j].z + vec1[2] * ix + vec2[2] * iy + vec3[2] * iz; FAtom atom(x, y, z, i, j, ix, iy, iz); int box_i_x, box_i_y, box_i_z; - //this->getBox(box_i_x, box_i_y, box_i_z, x, y, z); - box_i_x = ix + glayerX_minus; - box_i_y = iy + glayerY_minus; - box_i_z = iz + glayerZ_minus; + this->getBox(box_i_x, box_i_y, box_i_z, x, y, z); + box_i_x = clamp_box_index(box_i_x, this->box_nx); + box_i_y = clamp_box_index(box_i_y, this->box_ny); + box_i_z = clamp_box_index(box_i_z, this->box_nz); this->atoms_in_box[box_i_x][box_i_y][box_i_z].push_back(atom); } } @@ -269,17 +285,19 @@ void Grid::Construct_Adjacent(const UnitCell& ucell) { ModuleBase::timer::start("Grid", "constru_adj"); - for (int i_type = 0; i_type < ucell.ntype; i_type++) + for (int i_type = 0; i_type < ucell.ntype; i_type++) { +#ifdef _OPENMP +#pragma omp parallel for schedule(dynamic) +#endif for (int j_atom = 0; j_atom < ucell.atoms[i_type].na; j_atom++) { - FAtom atom(ucell.atoms[i_type].tau[j_atom].x, - ucell.atoms[i_type].tau[j_atom].y, - ucell.atoms[i_type].tau[j_atom].z, - i_type, - j_atom, - 0, 0 ,0); + ucell.atoms[i_type].tau[j_atom].y, + ucell.atoms[i_type].tau[j_atom].z, + i_type, + j_atom, + 0, 0, 0); this->Construct_Adjacent_near_box(atom); } @@ -289,17 +307,27 @@ void Grid::Construct_Adjacent(const UnitCell& ucell) void Grid::Construct_Adjacent_near_box(const FAtom& fatom) { - ModuleBase::timer::start("Grid", "adj_near_box"); int box_i_x=0; int box_i_y=0; int box_i_z=0; this->getBox(box_i_x, box_i_y, box_i_z, fatom.x, fatom.y, fatom.z); - - for (int box_i_x_adj = 0; box_i_x_adj < glayerX + glayerX_minus; box_i_x_adj++) + box_i_x = clamp_box_index(box_i_x, this->box_nx); + box_i_y = clamp_box_index(box_i_y, this->box_ny); + box_i_z = clamp_box_index(box_i_z, this->box_nz); + + const int search_layer = std::max(1, static_cast(std::ceil(this->sradius / this->box_edge_length))); + const int box_x_begin = std::max(0, box_i_x - search_layer); + const int box_x_end = std::min(this->box_nx - 1, box_i_x + search_layer); + const int box_y_begin = std::max(0, box_i_y - search_layer); + const int box_y_end = std::min(this->box_ny - 1, box_i_y + search_layer); + const int box_z_begin = std::max(0, box_i_z - search_layer); + const int box_z_end = std::min(this->box_nz - 1, box_i_z + search_layer); + + for (int box_i_x_adj = box_x_begin; box_i_x_adj <= box_x_end; box_i_x_adj++) { - for (int box_i_y_adj = 0; box_i_y_adj < glayerY + glayerY_minus; box_i_y_adj++) + for (int box_i_y_adj = box_y_begin; box_i_y_adj <= box_y_end; box_i_y_adj++) { - for (int box_i_z_adj = 0; box_i_z_adj < glayerZ + glayerZ_minus; box_i_z_adj++) + for (int box_i_z_adj = box_z_begin; box_i_z_adj <= box_z_end; box_i_z_adj++) { for (auto &fatom2 : this->atoms_in_box[box_i_x_adj][box_i_y_adj][box_i_z_adj]) { @@ -308,7 +336,6 @@ void Grid::Construct_Adjacent_near_box(const FAtom& fatom) } } } - ModuleBase::timer::end("Grid", "adj_near_box"); } void Grid::Construct_Adjacent_final(const FAtom& fatom1, From 95a4a0d44bf8f7197c323f9e94c9768c25d51904 Mon Sep 17 00:00:00 2001 From: Frostmoss <2400011223@stu.pku.edu.cn> Date: Thu, 4 Jun 2026 16:24:04 +0800 Subject: [PATCH 02/15] Add MPI domain decomposition interface --- .../source_cell/module_neighbor/mpi_domain.h | 103 ++++++++++++++++++ 1 file changed, 103 insertions(+) create mode 100644 source/source_cell/module_neighbor/mpi_domain.h diff --git a/source/source_cell/module_neighbor/mpi_domain.h b/source/source_cell/module_neighbor/mpi_domain.h new file mode 100644 index 00000000000..a1205b2d3b8 --- /dev/null +++ b/source/source_cell/module_neighbor/mpi_domain.h @@ -0,0 +1,103 @@ +#ifndef MODULE_NEIGHBOR_MPI_DOMAIN_H +#define MODULE_NEIGHBOR_MPI_DOMAIN_H + +#include +#include + +#ifdef __MPI +#include +#else +// Keep the interface buildable in serial configurations. Runtime MPI behavior is +// implemented only when ABACUS is configured with ENABLE_MPI. +typedef int MPI_Comm; +const MPI_Comm MPI_COMM_WORLD = 0; +const MPI_Comm MPI_COMM_NULL = -1; +#endif + +namespace ModuleNeighbor +{ + +struct DomainBounds +{ + std::array lower; + std::array upper; +}; + +struct MpiAtomRecord +{ + double x; + double y; + double z; + int global_index; + std::array pbc_shift; + bool is_ghost; + + MpiAtomRecord(); + MpiAtomRecord(const double x_in, + const double y_in, + const double z_in, + const int global_index_in, + const std::array& pbc_shift_in = std::array{{0, 0, 0}}, + const bool is_ghost_in = false); +}; + +class MpiDomain +{ + public: + MpiDomain(); + MpiDomain(const MpiDomain&) = delete; + MpiDomain& operator=(const MpiDomain&) = delete; + MpiDomain(MpiDomain&& other); + MpiDomain& operator=(MpiDomain&& other); + ~MpiDomain(); + + void initialize(MPI_Comm parent_comm, + const std::array& global_lower, + const std::array& global_upper, + const double ghost_cutoff, + const bool pbc); + + bool initialized() const; + int rank() const; + int size() const; + MPI_Comm cart_comm() const; + const std::array& dims() const; + const std::array& coords() const; + const std::array& periods() const; + const DomainBounds& local_bounds() const; + double ghost_cutoff() const; + + bool owns(const double x, const double y, const double z) const; + std::vector select_local_atoms(const std::vector& atoms) const; + + // Baseline ghost exchange: each rank publishes its owned atoms, then every + // rank keeps the remote periodic images that touch its local ghost shell. + // This is intentionally simple and correctness-first; neighbor-only MPI + // exchange can replace the internals without changing this interface. + std::vector exchange_ghost_atoms(const std::vector& local_atoms) const; + + private: + void reset(); + void compute_local_bounds(); + bool inside_local(const std::array& point) const; + bool inside_expanded(const std::array& point) const; + bool append_ghost_images(const MpiAtomRecord& atom, std::vector& ghosts) const; + + MPI_Comm cart_comm_; + bool owns_comm_; + bool initialized_; + int rank_; + int size_; + std::array dims_; + std::array coords_; + std::array periods_; + std::array global_lower_; + std::array global_upper_; + std::array global_length_; + DomainBounds local_bounds_; + double ghost_cutoff_; +}; + +} // namespace ModuleNeighbor + +#endif From 16afe77caabfee7a553813ba8cb2ff8ef755eea6 Mon Sep 17 00:00:00 2001 From: Frostmoss <2400011223@stu.pku.edu.cn> Date: Thu, 4 Jun 2026 16:25:15 +0800 Subject: [PATCH 03/15] Implement MPI domain decomposition --- .../module_neighbor/mpi_domain.cpp | 378 ++++++++++++++++++ 1 file changed, 378 insertions(+) create mode 100644 source/source_cell/module_neighbor/mpi_domain.cpp diff --git a/source/source_cell/module_neighbor/mpi_domain.cpp b/source/source_cell/module_neighbor/mpi_domain.cpp new file mode 100644 index 00000000000..0db5c35e908 --- /dev/null +++ b/source/source_cell/module_neighbor/mpi_domain.cpp @@ -0,0 +1,378 @@ +#include "mpi_domain.h" + +#include +#include +#include + +namespace ModuleNeighbor +{ +namespace +{ +const int kPackedAtomFields = 7; + +std::array make_point(const double x, const double y, const double z) +{ + return std::array{{x, y, z}}; +} + +bool in_half_open_range(const double value, + const double lower, + const double upper, + const bool include_upper) +{ + if (include_upper) + { + return value >= lower && value <= upper; + } + return value >= lower && value < upper; +} +} // namespace + +MpiAtomRecord::MpiAtomRecord() + : x(0.0), y(0.0), z(0.0), global_index(-1), pbc_shift{{0, 0, 0}}, is_ghost(false) +{ +} + +MpiAtomRecord::MpiAtomRecord(const double x_in, + const double y_in, + const double z_in, + const int global_index_in, + const std::array& pbc_shift_in, + const bool is_ghost_in) + : x(x_in), y(y_in), z(z_in), global_index(global_index_in), pbc_shift(pbc_shift_in), is_ghost(is_ghost_in) +{ +} + +MpiDomain::MpiDomain() + : cart_comm_(MPI_COMM_NULL), + owns_comm_(false), + initialized_(false), + rank_(0), + size_(1), + dims_{{1, 1, 1}}, + coords_{{0, 0, 0}}, + periods_{{0, 0, 0}}, + global_lower_{{0.0, 0.0, 0.0}}, + global_upper_{{0.0, 0.0, 0.0}}, + global_length_{{0.0, 0.0, 0.0}}, + local_bounds_{{{0.0, 0.0, 0.0}}, {{0.0, 0.0, 0.0}}}, + ghost_cutoff_(0.0) +{ +} + +MpiDomain::MpiDomain(MpiDomain&& other) : MpiDomain() +{ + *this = std::move(other); +} + +MpiDomain& MpiDomain::operator=(MpiDomain&& other) +{ + if (this != &other) + { + reset(); + cart_comm_ = other.cart_comm_; + owns_comm_ = other.owns_comm_; + initialized_ = other.initialized_; + rank_ = other.rank_; + size_ = other.size_; + dims_ = other.dims_; + coords_ = other.coords_; + periods_ = other.periods_; + global_lower_ = other.global_lower_; + global_upper_ = other.global_upper_; + global_length_ = other.global_length_; + local_bounds_ = other.local_bounds_; + ghost_cutoff_ = other.ghost_cutoff_; + + other.cart_comm_ = MPI_COMM_NULL; + other.owns_comm_ = false; + other.initialized_ = false; + } + return *this; +} + +MpiDomain::~MpiDomain() +{ + reset(); +} + +void MpiDomain::reset() +{ +#ifdef __MPI + int mpi_finalized = 0; + MPI_Finalized(&mpi_finalized); + if (!mpi_finalized && owns_comm_ && cart_comm_ != MPI_COMM_NULL) + { + MPI_Comm_free(&cart_comm_); + } +#endif + cart_comm_ = MPI_COMM_NULL; + owns_comm_ = false; + initialized_ = false; +} + +void MpiDomain::initialize(MPI_Comm parent_comm, + const std::array& global_lower, + const std::array& global_upper, + const double ghost_cutoff, + const bool pbc) +{ + reset(); + if (ghost_cutoff < 0.0) + { + throw std::invalid_argument("MpiDomain ghost cutoff must be non-negative."); + } + for (int axis = 0; axis < 3; ++axis) + { + if (!(global_upper[axis] > global_lower[axis])) + { + throw std::invalid_argument("MpiDomain global bounds must have positive length."); + } + } + + global_lower_ = global_lower; + global_upper_ = global_upper; + ghost_cutoff_ = ghost_cutoff; + periods_ = pbc ? std::array{{1, 1, 1}} : std::array{{0, 0, 0}}; + for (int axis = 0; axis < 3; ++axis) + { + global_length_[axis] = global_upper_[axis] - global_lower_[axis]; + } + +#ifdef __MPI + MPI_Comm_size(parent_comm, &size_); + MPI_Comm_rank(parent_comm, &rank_); + dims_ = std::array{{0, 0, 0}}; + MPI_Dims_create(size_, 3, dims_.data()); + MPI_Cart_create(parent_comm, 3, dims_.data(), periods_.data(), 0, &cart_comm_); + if (cart_comm_ == MPI_COMM_NULL) + { + throw std::runtime_error("MPI_Cart_create returned MPI_COMM_NULL."); + } + owns_comm_ = true; + MPI_Comm_rank(cart_comm_, &rank_); + MPI_Cart_coords(cart_comm_, rank_, 3, coords_.data()); +#else + (void)parent_comm; + size_ = 1; + rank_ = 0; + dims_ = std::array{{1, 1, 1}}; + coords_ = std::array{{0, 0, 0}}; + cart_comm_ = MPI_COMM_WORLD; + owns_comm_ = false; +#endif + + compute_local_bounds(); + initialized_ = true; +} + +void MpiDomain::compute_local_bounds() +{ + for (int axis = 0; axis < 3; ++axis) + { + const double width = global_length_[axis] / static_cast(dims_[axis]); + local_bounds_.lower[axis] = global_lower_[axis] + width * static_cast(coords_[axis]); + local_bounds_.upper[axis] = coords_[axis] == dims_[axis] - 1 + ? global_upper_[axis] + : global_lower_[axis] + width * static_cast(coords_[axis] + 1); + } +} + +bool MpiDomain::initialized() const +{ + return initialized_; +} + +int MpiDomain::rank() const +{ + return rank_; +} + +int MpiDomain::size() const +{ + return size_; +} + +MPI_Comm MpiDomain::cart_comm() const +{ + return cart_comm_; +} + +const std::array& MpiDomain::dims() const +{ + return dims_; +} + +const std::array& MpiDomain::coords() const +{ + return coords_; +} + +const std::array& MpiDomain::periods() const +{ + return periods_; +} + +const DomainBounds& MpiDomain::local_bounds() const +{ + return local_bounds_; +} + +double MpiDomain::ghost_cutoff() const +{ + return ghost_cutoff_; +} + +bool MpiDomain::inside_local(const std::array& point) const +{ + for (int axis = 0; axis < 3; ++axis) + { + const bool include_upper = coords_[axis] == dims_[axis] - 1; + if (!in_half_open_range(point[axis], local_bounds_.lower[axis], local_bounds_.upper[axis], include_upper)) + { + return false; + } + } + return true; +} + +bool MpiDomain::inside_expanded(const std::array& point) const +{ + for (int axis = 0; axis < 3; ++axis) + { + if (point[axis] < local_bounds_.lower[axis] - ghost_cutoff_ + || point[axis] > local_bounds_.upper[axis] + ghost_cutoff_) + { + return false; + } + } + return true; +} + +bool MpiDomain::owns(const double x, const double y, const double z) const +{ + return inside_local(make_point(x, y, z)); +} + +std::vector MpiDomain::select_local_atoms(const std::vector& atoms) const +{ + std::vector local_indices; + for (int i = 0; i < static_cast(atoms.size()); ++i) + { + if (owns(atoms[i].x, atoms[i].y, atoms[i].z)) + { + local_indices.push_back(i); + } + } + return local_indices; +} + +bool MpiDomain::append_ghost_images(const MpiAtomRecord& atom, std::vector& ghosts) const +{ + bool appended = false; + const int sx_begin = periods_[0] ? -1 : 0; + const int sx_end = periods_[0] ? 1 : 0; + const int sy_begin = periods_[1] ? -1 : 0; + const int sy_end = periods_[1] ? 1 : 0; + const int sz_begin = periods_[2] ? -1 : 0; + const int sz_end = periods_[2] ? 1 : 0; + + for (int sx = sx_begin; sx <= sx_end; ++sx) + { + for (int sy = sy_begin; sy <= sy_end; ++sy) + { + for (int sz = sz_begin; sz <= sz_end; ++sz) + { + MpiAtomRecord image = atom; + image.x += static_cast(sx) * global_length_[0]; + image.y += static_cast(sy) * global_length_[1]; + image.z += static_cast(sz) * global_length_[2]; + image.pbc_shift[0] += sx; + image.pbc_shift[1] += sy; + image.pbc_shift[2] += sz; + image.is_ghost = true; + + const std::array point = make_point(image.x, image.y, image.z); + if (inside_expanded(point) && !inside_local(point)) + { + ghosts.push_back(image); + appended = true; + } + } + } + } + return appended; +} + +std::vector MpiDomain::exchange_ghost_atoms(const std::vector& local_atoms) const +{ + if (!initialized_) + { + throw std::runtime_error("MpiDomain must be initialized before ghost exchange."); + } + +#ifdef __MPI + const int send_count = static_cast(local_atoms.size()) * kPackedAtomFields; + std::vector send_buffer(send_count, 0.0); + for (int i = 0; i < static_cast(local_atoms.size()); ++i) + { + const int offset = i * kPackedAtomFields; + send_buffer[offset + 0] = local_atoms[i].x; + send_buffer[offset + 1] = local_atoms[i].y; + send_buffer[offset + 2] = local_atoms[i].z; + send_buffer[offset + 3] = static_cast(local_atoms[i].global_index); + send_buffer[offset + 4] = static_cast(local_atoms[i].pbc_shift[0]); + send_buffer[offset + 5] = static_cast(local_atoms[i].pbc_shift[1]); + send_buffer[offset + 6] = static_cast(local_atoms[i].pbc_shift[2]); + } + + std::vector counts(size_, 0); + MPI_Allgather(&send_count, 1, MPI_INT, counts.data(), 1, MPI_INT, cart_comm_); + + std::vector displacements(size_, 0); + int total_count = 0; + for (int i = 0; i < size_; ++i) + { + displacements[i] = total_count; + total_count += counts[i]; + } + + std::vector recv_buffer(total_count, 0.0); + MPI_Allgatherv(send_buffer.empty() ? 0 : send_buffer.data(), + send_count, + MPI_DOUBLE, + recv_buffer.empty() ? 0 : recv_buffer.data(), + counts.data(), + displacements.data(), + MPI_DOUBLE, + cart_comm_); + + std::vector ghosts; + for (int owner_rank = 0; owner_rank < size_; ++owner_rank) + { + if (owner_rank == rank_) + { + continue; + } + for (int offset = displacements[owner_rank]; offset < displacements[owner_rank] + counts[owner_rank]; + offset += kPackedAtomFields) + { + MpiAtomRecord atom(recv_buffer[offset + 0], + recv_buffer[offset + 1], + recv_buffer[offset + 2], + static_cast(recv_buffer[offset + 3]), + std::array{{static_cast(recv_buffer[offset + 4]), + static_cast(recv_buffer[offset + 5]), + static_cast(recv_buffer[offset + 6])}}, + false); + append_ghost_images(atom, ghosts); + } + } + return ghosts; +#else + (void)local_atoms; + return std::vector(); +#endif +} + +} // namespace ModuleNeighbor From 9ff6fe6e106d43aa5371487f5f9a29ac88a4c204 Mon Sep 17 00:00:00 2001 From: Frostmoss <2400011223@stu.pku.edu.cn> Date: Thu, 4 Jun 2026 16:25:48 +0800 Subject: [PATCH 04/15] Add MPI domain decomposition tests --- .../module_neighbor/test/mpi_domain_test.cpp | 114 ++++++++++++++++++ 1 file changed, 114 insertions(+) create mode 100644 source/source_cell/module_neighbor/test/mpi_domain_test.cpp diff --git a/source/source_cell/module_neighbor/test/mpi_domain_test.cpp b/source/source_cell/module_neighbor/test/mpi_domain_test.cpp new file mode 100644 index 00000000000..3774c994500 --- /dev/null +++ b/source/source_cell/module_neighbor/test/mpi_domain_test.cpp @@ -0,0 +1,114 @@ +#include "../mpi_domain.h" + +#include "gtest/gtest.h" + +#ifdef __MPI +#include +#endif + +#include +#include + +namespace +{ +class ScopedMpiInit +{ + public: + ScopedMpiInit() : initialized_before_(false), should_finalize_(false) + { +#ifdef __MPI + int initialized = 0; + MPI_Initialized(&initialized); + initialized_before_ = initialized != 0; + if (!initialized_before_) + { + int argc = 0; + char** argv = nullptr; + MPI_Init(&argc, &argv); + should_finalize_ = true; + } +#endif + } + + ~ScopedMpiInit() + { +#ifdef __MPI + if (should_finalize_) + { + MPI_Finalize(); + } +#endif + } + + private: + bool initialized_before_; + bool should_finalize_; +}; +} // namespace + +TEST(MpiDomainTest, CreatesCartesianDomainAndOwnsLocalPoint) +{ + ScopedMpiInit mpi; + + ModuleNeighbor::MpiDomain domain; + domain.initialize(MPI_COMM_WORLD, + std::array{{0.0, 0.0, 0.0}}, + std::array{{8.0, 4.0, 2.0}}, + 0.5, + true); + + EXPECT_TRUE(domain.initialized()); + EXPECT_GE(domain.size(), 1); + EXPECT_EQ(domain.periods()[0], 1); + EXPECT_EQ(domain.periods()[1], 1); + EXPECT_EQ(domain.periods()[2], 1); + EXPECT_GT(domain.local_bounds().upper[0], domain.local_bounds().lower[0]); + EXPECT_GT(domain.local_bounds().upper[1], domain.local_bounds().lower[1]); + EXPECT_GT(domain.local_bounds().upper[2], domain.local_bounds().lower[2]); + + const double x_mid = 0.5 * (domain.local_bounds().lower[0] + domain.local_bounds().upper[0]); + const double y_mid = 0.5 * (domain.local_bounds().lower[1] + domain.local_bounds().upper[1]); + const double z_mid = 0.5 * (domain.local_bounds().lower[2] + domain.local_bounds().upper[2]); + EXPECT_TRUE(domain.owns(x_mid, y_mid, z_mid)); +} + +TEST(MpiDomainTest, SelectsLocalAtomsAndCountsGhostCandidates) +{ + ScopedMpiInit mpi; + + const double cutoff = 0.5; + ModuleNeighbor::MpiDomain domain; + domain.initialize(MPI_COMM_WORLD, + std::array{{0.0, 0.0, 0.0}}, + std::array{{8.0, 4.0, 2.0}}, + cutoff, + true); + + const double x_local = std::max(domain.local_bounds().lower[0], domain.local_bounds().upper[0] - 0.25 * cutoff); + const double y_local = 0.5 * (domain.local_bounds().lower[1] + domain.local_bounds().upper[1]); + const double z_local = 0.5 * (domain.local_bounds().lower[2] + domain.local_bounds().upper[2]); + + std::vector atoms; + atoms.push_back(ModuleNeighbor::MpiAtomRecord(x_local, y_local, z_local, domain.rank())); + + const std::vector local_indices = domain.select_local_atoms(atoms); + EXPECT_EQ(local_indices.size(), 1); + + const std::vector ghosts = domain.exchange_ghost_atoms(atoms); + +#ifdef __MPI + int local_ghost_count = static_cast(ghosts.size()); + int global_ghost_count = 0; + MPI_Allreduce(&local_ghost_count, &global_ghost_count, 1, MPI_INT, MPI_SUM, domain.cart_comm()); + if (domain.size() == 1) + { + EXPECT_EQ(global_ghost_count, 0); + } + else + { + EXPECT_GT(global_ghost_count, 0); + } +#else + EXPECT_TRUE(ghosts.empty()); +#endif +} From 9704c05f057f7c3ac1236c1ae91f8e38e529d46f Mon Sep 17 00:00:00 2001 From: Frostmoss <2400011223@stu.pku.edu.cn> Date: Thu, 4 Jun 2026 16:26:04 +0800 Subject: [PATCH 05/15] Build MPI domain decomposition component --- source/source_cell/module_neighbor/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/source/source_cell/module_neighbor/CMakeLists.txt b/source/source_cell/module_neighbor/CMakeLists.txt index 9b5fea4d147..78d7f553a51 100644 --- a/source/source_cell/module_neighbor/CMakeLists.txt +++ b/source/source_cell/module_neighbor/CMakeLists.txt @@ -5,6 +5,7 @@ add_library( sltk_atom_arrange.cpp sltk_grid.cpp sltk_grid_driver.cpp + mpi_domain.cpp ) if(ENABLE_COVERAGE) From 9a9602d7c1c89d1c0c70f4d3db454c2a4fc24d43 Mon Sep 17 00:00:00 2001 From: Frostmoss <2400011223@stu.pku.edu.cn> Date: Thu, 4 Jun 2026 16:26:17 +0800 Subject: [PATCH 06/15] Add MPI domain decomposition test target --- source/source_cell/module_neighbor/test/CMakeLists.txt | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/source/source_cell/module_neighbor/test/CMakeLists.txt b/source/source_cell/module_neighbor/test/CMakeLists.txt index 3d891aa20c3..d3c7d02f65f 100644 --- a/source/source_cell/module_neighbor/test/CMakeLists.txt +++ b/source/source_cell/module_neighbor/test/CMakeLists.txt @@ -22,4 +22,10 @@ AddTest( SOURCES sltk_atom_arrange_test.cpp ../sltk_atom_arrange.cpp ../sltk_grid_driver.cpp ../sltk_grid.cpp ../sltk_atom.cpp ../../../source_io/module_output/output.cpp -) \ No newline at end of file +) + +AddTest( + TARGET MODULE_CELL_NEIGHBOR_mpi_domain + LIBS ${math_libs} + SOURCES mpi_domain_test.cpp ../mpi_domain.cpp +) From d686797fe8a843782ce826e40347060a2e35c93d Mon Sep 17 00:00:00 2001 From: Frostmoss <2400011223@stu.pku.edu.cn> Date: Thu, 4 Jun 2026 16:27:55 +0800 Subject: [PATCH 07/15] Fix MPI domain C++11 includes and initialization --- source/source_cell/module_neighbor/mpi_domain.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/source/source_cell/module_neighbor/mpi_domain.cpp b/source/source_cell/module_neighbor/mpi_domain.cpp index 0db5c35e908..f3f535caa6f 100644 --- a/source/source_cell/module_neighbor/mpi_domain.cpp +++ b/source/source_cell/module_neighbor/mpi_domain.cpp @@ -3,6 +3,7 @@ #include #include #include +#include namespace ModuleNeighbor { @@ -55,9 +56,11 @@ MpiDomain::MpiDomain() global_lower_{{0.0, 0.0, 0.0}}, global_upper_{{0.0, 0.0, 0.0}}, global_length_{{0.0, 0.0, 0.0}}, - local_bounds_{{{0.0, 0.0, 0.0}}, {{0.0, 0.0, 0.0}}}, + local_bounds_(), ghost_cutoff_(0.0) { + local_bounds_.lower = std::array{{0.0, 0.0, 0.0}}; + local_bounds_.upper = std::array{{0.0, 0.0, 0.0}}; } MpiDomain::MpiDomain(MpiDomain&& other) : MpiDomain() From eea6eeae4e08dcb962404001445532e6530c7b34 Mon Sep 17 00:00:00 2001 From: Frostmoss <2400011223@stu.pku.edu.cn> Date: Thu, 4 Jun 2026 16:28:22 +0800 Subject: [PATCH 08/15] Fix MPI domain test includes --- source/source_cell/module_neighbor/test/mpi_domain_test.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/source/source_cell/module_neighbor/test/mpi_domain_test.cpp b/source/source_cell/module_neighbor/test/mpi_domain_test.cpp index 3774c994500..db8b11fafc1 100644 --- a/source/source_cell/module_neighbor/test/mpi_domain_test.cpp +++ b/source/source_cell/module_neighbor/test/mpi_domain_test.cpp @@ -6,6 +6,7 @@ #include #endif +#include #include #include From ddd25ed3048a5359e78ba583161d24d489b157b0 Mon Sep 17 00:00:00 2001 From: Frostmoss <2400011223@stu.pku.edu.cn> Date: Thu, 4 Jun 2026 16:30:02 +0800 Subject: [PATCH 09/15] Guard MPI communicator cleanup --- source/source_cell/module_neighbor/mpi_domain.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/source/source_cell/module_neighbor/mpi_domain.cpp b/source/source_cell/module_neighbor/mpi_domain.cpp index f3f535caa6f..48bef1a0cd8 100644 --- a/source/source_cell/module_neighbor/mpi_domain.cpp +++ b/source/source_cell/module_neighbor/mpi_domain.cpp @@ -102,11 +102,14 @@ MpiDomain::~MpiDomain() void MpiDomain::reset() { #ifdef __MPI - int mpi_finalized = 0; - MPI_Finalized(&mpi_finalized); - if (!mpi_finalized && owns_comm_ && cart_comm_ != MPI_COMM_NULL) + if (owns_comm_ && cart_comm_ != MPI_COMM_NULL) { - MPI_Comm_free(&cart_comm_); + int mpi_finalized = 0; + MPI_Finalized(&mpi_finalized); + if (!mpi_finalized) + { + MPI_Comm_free(&cart_comm_); + } } #endif cart_comm_ = MPI_COMM_NULL; From b7de4f915f215b63ee7b101946d640d08c60526c Mon Sep 17 00:00:00 2001 From: Frostmoss <2400011223@stu.pku.edu.cn> Date: Thu, 4 Jun 2026 16:36:15 +0800 Subject: [PATCH 10/15] Keep MPI initialized across domain tests --- .../module_neighbor/test/mpi_domain_test.cpp | 42 +++++-------------- 1 file changed, 11 insertions(+), 31 deletions(-) diff --git a/source/source_cell/module_neighbor/test/mpi_domain_test.cpp b/source/source_cell/module_neighbor/test/mpi_domain_test.cpp index db8b11fafc1..a6964eff773 100644 --- a/source/source_cell/module_neighbor/test/mpi_domain_test.cpp +++ b/source/source_cell/module_neighbor/test/mpi_domain_test.cpp @@ -12,44 +12,24 @@ namespace { -class ScopedMpiInit +void ensure_mpi_initialized() { - public: - ScopedMpiInit() : initialized_before_(false), should_finalize_(false) - { #ifdef __MPI - int initialized = 0; - MPI_Initialized(&initialized); - initialized_before_ = initialized != 0; - if (!initialized_before_) - { - int argc = 0; - char** argv = nullptr; - MPI_Init(&argc, &argv); - should_finalize_ = true; - } -#endif - } - - ~ScopedMpiInit() + int initialized = 0; + MPI_Initialized(&initialized); + if (!initialized) { -#ifdef __MPI - if (should_finalize_) - { - MPI_Finalize(); - } -#endif + int argc = 0; + char** argv = nullptr; + MPI_Init(&argc, &argv); } - - private: - bool initialized_before_; - bool should_finalize_; -}; +#endif +} } // namespace TEST(MpiDomainTest, CreatesCartesianDomainAndOwnsLocalPoint) { - ScopedMpiInit mpi; + ensure_mpi_initialized(); ModuleNeighbor::MpiDomain domain; domain.initialize(MPI_COMM_WORLD, @@ -75,7 +55,7 @@ TEST(MpiDomainTest, CreatesCartesianDomainAndOwnsLocalPoint) TEST(MpiDomainTest, SelectsLocalAtomsAndCountsGhostCandidates) { - ScopedMpiInit mpi; + ensure_mpi_initialized(); const double cutoff = 0.5; ModuleNeighbor::MpiDomain domain; From 891db22f6cee2279a5aaa09609a2774501390ecf Mon Sep 17 00:00:00 2001 From: Yanglijf Date: Fri, 5 Jun 2026 09:55:08 +0800 Subject: [PATCH 11/15] Add AtomPack grid storage and MPI neighbor prototype tests --- .../module_neighbor/CMakeLists.txt | 7 +- .../source_cell/module_neighbor/atom_pack.cpp | 304 ++++++++++++++++++ .../source_cell/module_neighbor/atom_pack.h | 87 +++++ .../module_neighbor/test/CMakeLists.txt | 13 + .../module_neighbor/test/atom_pack_test.cpp | 217 +++++++++++++ .../module_neighbor/test/mpi_domain_test.cpp | 26 ++ .../test/mpi_neighbor_proto_test.cpp | 116 +++++++ .../module_neighbor/test/sltk_grid_test.cpp | 56 ++++ 8 files changed, 822 insertions(+), 4 deletions(-) create mode 100644 source/source_cell/module_neighbor/atom_pack.cpp create mode 100644 source/source_cell/module_neighbor/atom_pack.h create mode 100644 source/source_cell/module_neighbor/test/atom_pack_test.cpp create mode 100644 source/source_cell/module_neighbor/test/mpi_neighbor_proto_test.cpp diff --git a/source/source_cell/module_neighbor/CMakeLists.txt b/source/source_cell/module_neighbor/CMakeLists.txt index 78d7f553a51..65a1dd7330c 100644 --- a/source/source_cell/module_neighbor/CMakeLists.txt +++ b/source/source_cell/module_neighbor/CMakeLists.txt @@ -5,6 +5,7 @@ add_library( sltk_atom_arrange.cpp sltk_grid.cpp sltk_grid_driver.cpp + atom_pack.cpp mpi_domain.cpp ) @@ -13,7 +14,5 @@ if(ENABLE_COVERAGE) endif() if(BUILD_TESTING) - if(ENABLE_MPI) - add_subdirectory(test) - endif() -endif() \ No newline at end of file + add_subdirectory(test) +endif() diff --git a/source/source_cell/module_neighbor/atom_pack.cpp b/source/source_cell/module_neighbor/atom_pack.cpp new file mode 100644 index 00000000000..17feddc2819 --- /dev/null +++ b/source/source_cell/module_neighbor/atom_pack.cpp @@ -0,0 +1,304 @@ +#include "atom_pack.h" + +#include +#include +#include +#include + +namespace ModuleNeighbor +{ +namespace +{ +struct ExpandLayers +{ + int x_plus; + int x_minus; + int y_plus; + int y_minus; + int z_plus; + int z_minus; +}; + +ExpandLayers compute_expand_layers(const UnitCell& ucell, const double radius_lat0, const bool pbc) +{ + if (radius_lat0 < 0.0) + { + throw std::invalid_argument("AtomPack radius must be non-negative."); + } + if (!pbc) + { + // Non-periodic search should not manufacture image atoms. The loop bounds + // below use [-minus, plus), so {plus=1, minus=0} keeps only shift (0,0,0). + return ExpandLayers{1, 0, 1, 0, 1, 0}; + } + + // Use the same projection-based layer estimate as Grid::Check_Expand_Condition(). + // The three cross products measure the cell-face areas used to convert a + // Cartesian cutoff into lattice-vector image counts for skewed cells. + const double a23_1 = ucell.latvec.e22 * ucell.latvec.e33 - ucell.latvec.e23 * ucell.latvec.e32; + const double a23_2 = ucell.latvec.e21 * ucell.latvec.e33 - ucell.latvec.e23 * ucell.latvec.e31; + const double a23_3 = ucell.latvec.e21 * ucell.latvec.e32 - ucell.latvec.e22 * ucell.latvec.e31; + const double a23_norm = std::sqrt(a23_1 * a23_1 + a23_2 * a23_2 + a23_3 * a23_3); + const int extend_x = static_cast( + std::ceil(a23_norm * radius_lat0 / ucell.omega * ucell.lat0 * ucell.lat0 * ucell.lat0)); + + const double a31_1 = ucell.latvec.e32 * ucell.latvec.e13 - ucell.latvec.e33 * ucell.latvec.e12; + const double a31_2 = ucell.latvec.e31 * ucell.latvec.e13 - ucell.latvec.e33 * ucell.latvec.e11; + const double a31_3 = ucell.latvec.e31 * ucell.latvec.e12 - ucell.latvec.e32 * ucell.latvec.e11; + const double a31_norm = std::sqrt(a31_1 * a31_1 + a31_2 * a31_2 + a31_3 * a31_3); + const int extend_y = static_cast( + std::ceil(a31_norm * radius_lat0 / ucell.omega * ucell.lat0 * ucell.lat0 * ucell.lat0)); + + const double a12_1 = ucell.latvec.e12 * ucell.latvec.e23 - ucell.latvec.e13 * ucell.latvec.e22; + const double a12_2 = ucell.latvec.e11 * ucell.latvec.e23 - ucell.latvec.e13 * ucell.latvec.e21; + const double a12_3 = ucell.latvec.e11 * ucell.latvec.e22 - ucell.latvec.e12 * ucell.latvec.e21; + const double a12_norm = std::sqrt(a12_1 * a12_1 + a12_2 * a12_2 + a12_3 * a12_3); + const int extend_z = static_cast( + std::ceil(a12_norm * radius_lat0 / ucell.omega * ucell.lat0 * ucell.lat0 * ucell.lat0)); + + // Grid expands images with loops over [-minus, plus). Keeping plus one layer + // larger preserves the current asymmetric loop convention and test comparability. + return ExpandLayers{extend_x + 1, extend_x, extend_y + 1, extend_y, extend_z + 1, extend_z}; +} + +int count_images(const ExpandLayers& layers) +{ + return (layers.x_plus + layers.x_minus) * (layers.y_plus + layers.y_minus) + * (layers.z_plus + layers.z_minus); +} + +int clamp_index(const int index, const int size) +{ + if (index < 0) + { + return 0; + } + if (index >= size) + { + return size - 1; + } + return index; +} +} // namespace + +void AtomPack::clear() +{ + x.clear(); + y.clear(); + z.clear(); + type.clear(); + natom.clear(); + cell_x.clear(); + cell_y.clear(); + cell_z.clear(); + global_index.clear(); + is_ghost.clear(); +} + +int AtomPack::size() const +{ + return static_cast(x.size()); +} + +bool AtomPack::empty() const +{ + return x.empty(); +} + +void AtomPack::reserve(const int count) +{ + x.reserve(count); + y.reserve(count); + z.reserve(count); + type.reserve(count); + natom.reserve(count); + cell_x.reserve(count); + cell_y.reserve(count); + cell_z.reserve(count); + global_index.reserve(count); + is_ghost.reserve(count); +} + +void AtomPack::append_atom(const double x_in, + const double y_in, + const double z_in, + const int type_in, + const int natom_in, + const int cell_x_in, + const int cell_y_in, + const int cell_z_in, + const int global_index_in, + const bool is_ghost_in) +{ + x.push_back(x_in); + y.push_back(y_in); + z.push_back(z_in); + type.push_back(type_in); + natom.push_back(natom_in); + cell_x.push_back(cell_x_in); + cell_y.push_back(cell_y_in); + cell_z.push_back(cell_z_in); + global_index.push_back(global_index_in); + is_ghost.push_back(is_ghost_in); +} + +void AtomPack::append_mpi_record(const MpiAtomRecord& record, const int type_in, const int natom_in) +{ + append_atom(record.x, + record.y, + record.z, + type_in, + natom_in, + record.pbc_shift[0], + record.pbc_shift[1], + record.pbc_shift[2], + record.global_index, + record.is_ghost); +} + +void GridStorage::clear() +{ + atoms_in_box.clear(); + box_offset.clear(); + box_count.clear(); + box_nx = 0; + box_ny = 0; + box_nz = 0; + x_min = 0.0; + y_min = 0.0; + z_min = 0.0; + x_max = 0.0; + y_max = 0.0; + z_max = 0.0; + box_edge_length = 0.0; +} + +int GridStorage::box_size() const +{ + return box_nx * box_ny * box_nz; +} + +int GridStorage::get_box_id(const int bx, const int by, const int bz) const +{ + if (box_nx <= 0 || box_ny <= 0 || box_nz <= 0) + { + throw std::runtime_error("GridStorage has not been initialized."); + } + const int bx_safe = clamp_index(bx, box_nx); + const int by_safe = clamp_index(by, box_ny); + const int bz_safe = clamp_index(bz, box_nz); + return (bx_safe * box_ny + by_safe) * box_nz + bz_safe; +} + +int GridStorage::get_box_id_from_coord(const double x, const double y, const double z) const +{ + if (box_edge_length <= 0.0) + { + throw std::runtime_error("GridStorage box edge length must be positive."); + } + const int bx = static_cast(std::floor((x - x_min) / box_edge_length)); + const int by = static_cast(std::floor((y - y_min) / box_edge_length)); + const int bz = static_cast(std::floor((z - z_min) / box_edge_length)); + return get_box_id(bx, by, bz); +} + +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); + + // UnitCell::tau is already stored in lattice-constant units after check_dtau(). + // A periodic image is therefore obtained by adding integer multiples of latvec. + 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); + + AtomPack pack; + pack.reserve(ucell.nat * count_images(layers)); + + int global_index_base = 0; + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + // Periodic images are search candidates only. They keep the original + // (type, natom, global_index) identity so neighbor results can be mapped + // back to the physical atom after distance checks. + for (int ix = -layers.x_minus; ix < layers.x_plus; ++ix) + { + for (int iy = -layers.y_minus; iy < layers.y_plus; ++iy) + { + for (int iz = -layers.z_minus; iz < layers.z_plus; ++iz) + { + const double x = ucell.atoms[it].tau[ia].x + vec1[0] * ix + vec2[0] * iy + vec3[0] * iz; + const double y = ucell.atoms[it].tau[ia].y + vec1[1] * ix + vec2[1] * iy + vec3[1] * iz; + const double z = ucell.atoms[it].tau[ia].z + vec1[2] * ix + vec2[2] * iy + vec3[2] * iz; + pack.append_atom(x, y, z, it, ia, ix, iy, iz, global_index_base + ia, false); + } + } + } + } + global_index_base += ucell.atoms[it].na; + } + + return pack; +} + +GridStorage build_grid_storage_from_atom_pack(const AtomPack& pack, const double box_edge_length) +{ + if (box_edge_length <= 0.0) + { + throw std::invalid_argument("GridStorage box edge length must be positive."); + } + if (pack.empty()) + { + throw std::invalid_argument("GridStorage cannot be built from an empty AtomPack."); + } + + GridStorage storage; + storage.box_edge_length = box_edge_length; + storage.x_min = storage.x_max = pack.x[0]; + storage.y_min = storage.y_max = pack.y[0]; + storage.z_min = storage.z_max = pack.z[0]; + + for (int i = 1; i < pack.size(); ++i) + { + storage.x_min = std::min(storage.x_min, pack.x[i]); + storage.x_max = std::max(storage.x_max, pack.x[i]); + storage.y_min = std::min(storage.y_min, pack.y[i]); + storage.y_max = std::max(storage.y_max, pack.y[i]); + storage.z_min = std::min(storage.z_min, pack.z[i]); + storage.z_max = std::max(storage.z_max, pack.z[i]); + } + + storage.box_nx = std::max(1, static_cast(std::floor((storage.x_max - storage.x_min) / box_edge_length)) + 1); + storage.box_ny = std::max(1, static_cast(std::floor((storage.y_max - storage.y_min) / box_edge_length)) + 1); + storage.box_nz = std::max(1, static_cast(std::floor((storage.z_max - storage.z_min) / box_edge_length)) + 1); + + const int nbox = storage.box_size(); + storage.box_count.assign(nbox, 0); + // First pass: count atoms per box. The count array is later reused as the + // stable box length table, so it must not be modified during insertion. + for (int i = 0; i < pack.size(); ++i) + { + ++storage.box_count[storage.get_box_id_from_coord(pack.x[i], pack.y[i], pack.z[i])]; + } + + // Prefix sum converts counts into disjoint ranges in atoms_in_box. + // box_offset[0] stays 0; box_offset[b] is the sum of all previous counts. + storage.box_offset.assign(nbox, 0); + std::partial_sum(storage.box_count.begin(), storage.box_count.end() - 1, storage.box_offset.begin() + 1); + + storage.atoms_in_box.assign(pack.size(), -1); + std::vector next_offset = storage.box_offset; + // Second pass: write AtomPack indices into their box ranges. next_offset is + // a scratch cursor so box_offset remains the public range table. + for (int i = 0; i < pack.size(); ++i) + { + const int box_id = storage.get_box_id_from_coord(pack.x[i], pack.y[i], pack.z[i]); + storage.atoms_in_box[next_offset[box_id]++] = i; + } + + return storage; +} + +} // namespace ModuleNeighbor diff --git a/source/source_cell/module_neighbor/atom_pack.h b/source/source_cell/module_neighbor/atom_pack.h new file mode 100644 index 00000000000..c2db5daddbe --- /dev/null +++ b/source/source_cell/module_neighbor/atom_pack.h @@ -0,0 +1,87 @@ +#ifndef MODULE_NEIGHBOR_ATOM_PACK_H +#define MODULE_NEIGHBOR_ATOM_PACK_H + +#include "mpi_domain.h" +#include "source_cell/unitcell.h" + +#include + +namespace ModuleNeighbor +{ + +struct AtomPack +{ + // A compact SoA container for neighbor-search input atoms. + // All vectors are kept at the same length and index i refers to one atom record. + // This layout is intended to replace pointer-heavy temporary records in later + // cell-list and MPI paths. + std::vector x; + std::vector y; + std::vector z; + std::vector type; + std::vector natom; + + // Periodic image shift relative to the origin cell. The origin atom uses + // (0, 0, 0); ghost atoms reuse the shift carried by MpiAtomRecord. + std::vector cell_x; + std::vector cell_y; + std::vector cell_z; + + // global_index identifies the original atom before PBC image expansion. + // is_ghost separates atoms imported from neighboring MPI domains from local atoms. + std::vector global_index; + std::vector is_ghost; + + void clear(); + int size() const; + bool empty() const; + void reserve(const int count); + void append_atom(const double x_in, + const double y_in, + const double z_in, + const int type_in, + const int natom_in, + const int cell_x_in, + const int cell_y_in, + const int cell_z_in, + const int global_index_in, + const bool is_ghost_in); + // Convert one MpiAtomRecord into the same layout used by local atoms, so local + // and ghost atoms can be traversed by one later neighbor-search kernel. + void append_mpi_record(const MpiAtomRecord& record, const int type_in = -1, const int natom_in = -1); +}; + +struct GridStorage +{ + // Flat box storage for AtomPack indices. Atoms in box b are stored in + // atoms_in_box[box_offset[b], box_offset[b] + box_count[b]). + std::vector atoms_in_box; + std::vector box_offset; + std::vector box_count; + + int box_nx = 0; + int box_ny = 0; + int box_nz = 0; + double x_min = 0.0; + double y_min = 0.0; + double z_min = 0.0; + double x_max = 0.0; + double y_max = 0.0; + double z_max = 0.0; + double box_edge_length = 0.0; + + void clear(); + int box_size() const; + int get_box_id(const int bx, const int by, const int bz) const; + int get_box_id_from_coord(const double x, const double y, const double z) 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); + +} // namespace ModuleNeighbor + +#endif diff --git a/source/source_cell/module_neighbor/test/CMakeLists.txt b/source/source_cell/module_neighbor/test/CMakeLists.txt index d3c7d02f65f..c4d055ddcb0 100644 --- a/source/source_cell/module_neighbor/test/CMakeLists.txt +++ b/source/source_cell/module_neighbor/test/CMakeLists.txt @@ -29,3 +29,16 @@ AddTest( LIBS ${math_libs} SOURCES mpi_domain_test.cpp ../mpi_domain.cpp ) + +AddTest( + TARGET MODULE_CELL_NEIGHBOR_atom_pack + LIBS parameter ${math_libs} base device cell_info + SOURCES atom_pack_test.cpp ../atom_pack.cpp ../mpi_domain.cpp ../sltk_grid.cpp ../sltk_atom.cpp + ../../../source_io/module_output/output.cpp +) + +AddTest( + TARGET MODULE_CELL_NEIGHBOR_mpi_neighbor_proto + LIBS ${math_libs} + SOURCES mpi_neighbor_proto_test.cpp ../atom_pack.cpp ../mpi_domain.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 new file mode 100644 index 00000000000..a3471a572da --- /dev/null +++ b/source/source_cell/module_neighbor/test/atom_pack_test.cpp @@ -0,0 +1,217 @@ +#include "../atom_pack.h" +#include "../sltk_grid.h" + +#include "gtest/gtest.h" +#include "prepare_unitcell.h" + +#include "source_cell/read_stru.h" + +#include +#include +#include + +#ifdef __LCAO +InfoNonlocal::InfoNonlocal() +{ +} +InfoNonlocal::~InfoNonlocal() +{ +} +LCAO_Orbitals::LCAO_Orbitals() +{ +} +LCAO_Orbitals::~LCAO_Orbitals() +{ +} +#endif +Magnetism::Magnetism() +{ + this->tot_mag = 0.0; + this->abs_mag = 0.0; + this->start_mag = nullptr; +} +Magnetism::~Magnetism() +{ + delete[] this->start_mag; +} + +namespace +{ +int count_origin_cell_atoms(const ModuleNeighbor::AtomPack& pack) +{ + int count = 0; + for (int i = 0; i < pack.size(); ++i) + { + if (pack.cell_x[i] == 0 && pack.cell_y[i] == 0 && pack.cell_z[i] == 0) + { + ++count; + } + } + return count; +} + +std::vector collect_grid_box_counts(const Grid& grid) +{ + std::vector counts; + counts.reserve(grid.box_nx * grid.box_ny * grid.box_nz); + for (int bx = 0; bx < grid.box_nx; ++bx) + { + for (int by = 0; by < grid.box_ny; ++by) + { + for (int bz = 0; bz < grid.box_nz; ++bz) + { + counts.push_back(static_cast(grid.atoms_in_box[bx][by][bz].size())); + } + } + } + return counts; +} +} // namespace + +TEST(AtomPackTest, BuildsNonPeriodicPackWithoutImages) +{ + UcellTestPrepare utp = UcellTestLib["Si"]; + UnitCell* ucell = utp.SetUcellInfo(); + unitcell::check_dtau(ucell->atoms, ucell->ntype, ucell->lat0, ucell->latvec); + + const ModuleNeighbor::AtomPack pack = ModuleNeighbor::build_atom_pack_from_unitcell(*ucell, 0.5, false); + + EXPECT_EQ(pack.size(), ucell->nat); + EXPECT_EQ(count_origin_cell_atoms(pack), ucell->nat); + EXPECT_FALSE(pack.empty()); + for (int i = 0; i < pack.size(); ++i) + { + EXPECT_FALSE(pack.is_ghost[i]); + EXPECT_EQ(pack.cell_x[i], 0); + EXPECT_EQ(pack.cell_y[i], 0); + EXPECT_EQ(pack.cell_z[i], 0); + } + + delete ucell; +} + +TEST(AtomPackTest, BuildsPeriodicPackWithImages) +{ + UcellTestPrepare utp = UcellTestLib["Si"]; + UnitCell* ucell = utp.SetUcellInfo(); + unitcell::check_dtau(ucell->atoms, ucell->ntype, ucell->lat0, ucell->latvec); + + const ModuleNeighbor::AtomPack pack = ModuleNeighbor::build_atom_pack_from_unitcell(*ucell, 0.5, true); + + EXPECT_GT(pack.size(), ucell->nat); + EXPECT_EQ(count_origin_cell_atoms(pack), ucell->nat); + EXPECT_EQ(pack.type.size(), pack.x.size()); + EXPECT_EQ(pack.global_index.size(), pack.x.size()); + + bool saw_shifted_image = false; + for (int i = 0; i < pack.size(); ++i) + { + saw_shifted_image = saw_shifted_image || pack.cell_x[i] != 0 || pack.cell_y[i] != 0 || pack.cell_z[i] != 0; + EXPECT_FALSE(pack.is_ghost[i]); + EXPECT_GE(pack.global_index[i], 0); + EXPECT_LT(pack.global_index[i], ucell->nat); + } + EXPECT_TRUE(saw_shifted_image); + + delete ucell; +} + +TEST(AtomPackTest, AppendsMpiGhostRecord) +{ + ModuleNeighbor::AtomPack pack; + pack.append_atom(0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, false); + + const ModuleNeighbor::MpiAtomRecord ghost(1.0, 2.0, 3.0, 7, std::array{{1, 0, -1}}, true); + pack.append_mpi_record(ghost); + + ASSERT_EQ(pack.size(), 2); + EXPECT_TRUE(pack.is_ghost[1]); + EXPECT_EQ(pack.global_index[1], 7); + EXPECT_EQ(pack.cell_x[1], 1); + EXPECT_EQ(pack.cell_y[1], 0); + EXPECT_EQ(pack.cell_z[1], -1); + EXPECT_DOUBLE_EQ(pack.x[1], 1.0); + EXPECT_DOUBLE_EQ(pack.y[1], 2.0); + EXPECT_DOUBLE_EQ(pack.z[1], 3.0); +} + +TEST(AtomPackTest, BuildsFlatGridStorageWithoutLossOrDuplication) +{ + 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); + + ASSERT_EQ(storage.atoms_in_box.size(), static_cast(pack.size())); + ASSERT_EQ(storage.box_offset.size(), static_cast(storage.box_size())); + ASSERT_EQ(storage.box_count.size(), static_cast(storage.box_size())); + EXPECT_GT(storage.box_nx, 0); + EXPECT_GT(storage.box_ny, 0); + EXPECT_GT(storage.box_nz, 0); + + std::vector visited(pack.size(), false); + for (int box_id = 0; box_id < storage.box_size(); ++box_id) + { + const int begin = storage.box_offset[box_id]; + const int end = begin + storage.box_count[box_id]; + ASSERT_GE(begin, 0); + ASSERT_GE(end, begin); + ASSERT_LE(end, static_cast(storage.atoms_in_box.size())); + + for (int offset = begin; offset < end; ++offset) + { + const int atom_index = storage.atoms_in_box[offset]; + ASSERT_GE(atom_index, 0); + ASSERT_LT(atom_index, pack.size()); + EXPECT_FALSE(visited[atom_index]); + visited[atom_index] = true; + EXPECT_EQ(storage.get_box_id_from_coord(pack.x[atom_index], pack.y[atom_index], pack.z[atom_index]), box_id); + } + } + + EXPECT_TRUE(std::all_of(visited.begin(), visited.end(), [](const bool value) { return value; })); + + delete ucell; +} + +TEST(AtomPackTest, GridStorageClampsOutOfRangeCoordinates) +{ + ModuleNeighbor::AtomPack pack; + pack.append_atom(0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, false); + pack.append_atom(1.0, 1.0, 1.0, 0, 1, 0, 0, 0, 1, false); + + const ModuleNeighbor::GridStorage storage = ModuleNeighbor::build_grid_storage_from_atom_pack(pack, 0.6); + + EXPECT_EQ(storage.get_box_id(-100, -100, -100), 0); + EXPECT_EQ(storage.get_box_id_from_coord(storage.x_min - 100.0, storage.y_min - 100.0, storage.z_min - 100.0), 0); + EXPECT_EQ(storage.get_box_id(100, 100, 100), storage.box_size() - 1); + EXPECT_EQ(storage.get_box_id_from_coord(storage.x_max + 100.0, storage.y_max + 100.0, storage.z_max + 100.0), + storage.box_size() - 1); +} + +TEST(AtomPackTest, FlatGridStorageMatchesLegacyGridBoxCountsWithoutPbc) +{ + 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); + + std::ofstream ofs("atom_pack_grid_compare.out"); + Grid grid(0); + grid.init(ofs, *ucell, radius, false); + ofs.close(); + + EXPECT_EQ(storage.box_nx, grid.box_nx); + EXPECT_EQ(storage.box_ny, grid.box_ny); + EXPECT_EQ(storage.box_nz, grid.box_nz); + EXPECT_EQ(storage.box_count, collect_grid_box_counts(grid)); + + remove("atom_pack_grid_compare.out"); + delete ucell; +} diff --git a/source/source_cell/module_neighbor/test/mpi_domain_test.cpp b/source/source_cell/module_neighbor/test/mpi_domain_test.cpp index a6964eff773..de82e2edfe4 100644 --- a/source/source_cell/module_neighbor/test/mpi_domain_test.cpp +++ b/source/source_cell/module_neighbor/test/mpi_domain_test.cpp @@ -93,3 +93,29 @@ TEST(MpiDomainTest, SelectsLocalAtomsAndCountsGhostCandidates) EXPECT_TRUE(ghosts.empty()); #endif } + +int main(int argc, char** argv) +{ + ::testing::InitGoogleTest(&argc, argv); +#ifdef __MPI + // These tests call MPI collectives inside the test body. Manage MPI here + // instead of relying on gtest_main, otherwise mpirun can report an abnormal + // exit even when all gtest assertions pass. + int initialized = 0; + MPI_Initialized(&initialized); + if (!initialized) + { + MPI_Init(&argc, &argv); + } +#endif + const int result = RUN_ALL_TESTS(); +#ifdef __MPI + int finalized = 0; + MPI_Finalized(&finalized); + if (!finalized) + { + MPI_Finalize(); + } +#endif + return result; +} diff --git a/source/source_cell/module_neighbor/test/mpi_neighbor_proto_test.cpp b/source/source_cell/module_neighbor/test/mpi_neighbor_proto_test.cpp new file mode 100644 index 00000000000..d4bb7e34b06 --- /dev/null +++ b/source/source_cell/module_neighbor/test/mpi_neighbor_proto_test.cpp @@ -0,0 +1,116 @@ +#include "../atom_pack.h" +#include "../mpi_domain.h" + +#include "gtest/gtest.h" + +#ifdef __MPI +#include +#endif + +#include +#include +#include + +namespace +{ +void ensure_mpi_initialized() +{ +#ifdef __MPI + int initialized = 0; + MPI_Initialized(&initialized); + if (!initialized) + { + int argc = 0; + char** argv = nullptr; + MPI_Init(&argc, &argv); + } +#endif +} + +int count_ghost_atoms(const ModuleNeighbor::AtomPack& pack) +{ + int count = 0; + for (int i = 0; i < pack.size(); ++i) + { + if (pack.is_ghost[i]) + { + ++count; + } + } + return count; +} +} // namespace + +TEST(MpiNeighborPrototypeTest, AppendsGhostsToAtomPack) +{ + ensure_mpi_initialized(); + + const double cutoff = 0.75; + ModuleNeighbor::MpiDomain domain; + domain.initialize(MPI_COMM_WORLD, + std::array{{0.0, 0.0, 0.0}}, + std::array{{8.0, 4.0, 2.0}}, + cutoff, + true); + + const double x_local = std::max(domain.local_bounds().lower[0], domain.local_bounds().upper[0] - 0.25 * cutoff); + const double y_local = 0.5 * (domain.local_bounds().lower[1] + domain.local_bounds().upper[1]); + const double z_local = 0.5 * (domain.local_bounds().lower[2] + domain.local_bounds().upper[2]); + + std::vector local_records; + local_records.push_back(ModuleNeighbor::MpiAtomRecord(x_local, y_local, z_local, domain.rank())); + + ModuleNeighbor::AtomPack pack; + pack.append_mpi_record(local_records.front(), 0, domain.rank()); + + const std::vector ghosts = domain.exchange_ghost_atoms(local_records); + for (const ModuleNeighbor::MpiAtomRecord& ghost: ghosts) + { + pack.append_mpi_record(ghost); + } + + EXPECT_EQ(pack.size(), static_cast(ghosts.size()) + 1); + EXPECT_EQ(count_ghost_atoms(pack), static_cast(ghosts.size())); + +#ifdef __MPI + int local_ghost_count = static_cast(ghosts.size()); + int global_ghost_count = 0; + MPI_Allreduce(&local_ghost_count, &global_ghost_count, 1, MPI_INT, MPI_SUM, domain.cart_comm()); + if (domain.size() == 1) + { + EXPECT_EQ(global_ghost_count, 0); + } + else + { + EXPECT_GT(global_ghost_count, 0); + } +#else + EXPECT_TRUE(ghosts.empty()); +#endif +} + +int main(int argc, char** argv) +{ + ::testing::InitGoogleTest(&argc, argv); +#ifdef __MPI + // These tests call MPI collectives inside the test body. Manage MPI here + // instead of relying on gtest_main, otherwise mpirun can report an abnormal + // exit even when all gtest assertions pass. + int initialized = 0; + MPI_Initialized(&initialized); + if (!initialized) + { + MPI_Init(&argc, &argv); + } +#endif + const int result = RUN_ALL_TESTS(); +#ifdef __MPI + int finalized = 0; + MPI_Finalized(&finalized); + if (!finalized) + { + MPI_Finalize(); + } +#endif + return result; +} 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 044feafc2de..4cd27cd8772 100644 --- a/source/source_cell/module_neighbor/test/sltk_grid_test.cpp +++ b/source/source_cell/module_neighbor/test/sltk_grid_test.cpp @@ -1,6 +1,14 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" +#include +#include +#include + +#ifdef _OPENMP +#include +#endif + #define private public #include "../sltk_grid.h" #include "prepare_unitcell.h" @@ -74,6 +82,28 @@ class SltkGridTest : public testing::Test using SltkGridDeathTest = SltkGridTest; +namespace +{ +using NeighborKey = std::tuple; + +std::vector collect_neighbor_keys(const Grid& grid) +{ + std::vector keys; + for (const auto& type_neighbors: grid.all_adj_info) + { + for (const auto& atom_neighbors: type_neighbors) + { + for (const FAtom* atom: atom_neighbors) + { + keys.push_back(NeighborKey(atom->type, atom->natom, atom->cell_x, atom->cell_y, atom->cell_z)); + } + } + } + std::sort(keys.begin(), keys.end()); + return keys; +} +} // namespace + TEST_F(SltkGridTest, Init) { ofs.open("test.out"); @@ -122,6 +152,32 @@ TEST_F(SltkGridTest, InitSmall) remove("test.out"); } +TEST_F(SltkGridTest, OpenMPThreadCountKeepsNeighborSet) +{ + unitcell::check_dtau(ucell->atoms, ucell->ntype, ucell->lat0, ucell->latvec); + + std::ofstream ofs_one("test_one.out"); + Grid grid_one(PARAM.input.test_grid); +#ifdef _OPENMP + omp_set_num_threads(1); +#endif + grid_one.init(ofs_one, *ucell, radius, pbc); + ofs_one.close(); + + std::ofstream ofs_many("test_many.out"); + Grid grid_many(PARAM.input.test_grid); +#ifdef _OPENMP + omp_set_num_threads(4); +#endif + grid_many.init(ofs_many, *ucell, radius, pbc); + ofs_many.close(); + + EXPECT_EQ(collect_neighbor_keys(grid_many), collect_neighbor_keys(grid_one)); + + remove("test_one.out"); + remove("test_many.out"); +} + /* // This test cannot pass because setAtomLinkArray() is unsuccessful // if expand_flag is false From 5ebdfae3487e827877978638db3a8265482a3744 Mon Sep 17 00:00:00 2001 From: hhb2006 <1035516541@qq.com> Date: Sat, 6 Jun 2026 11:58:48 +0800 Subject: [PATCH 12/15] Improve neighbor grid pair comparison test --- .../source_cell/module_neighbor/sltk_grid.h | 2 +- .../module_neighbor/test/sltk_grid_test.cpp | 22 +++++++++++++------ 2 files changed, 16 insertions(+), 8 deletions(-) diff --git a/source/source_cell/module_neighbor/sltk_grid.h b/source/source_cell/module_neighbor/sltk_grid.h index b5c3e836295..8faf75fa2ea 100644 --- a/source/source_cell/module_neighbor/sltk_grid.h +++ b/source/source_cell/module_neighbor/sltk_grid.h @@ -39,7 +39,7 @@ class Grid double z_max=0.0; // The algorithm for searching neighboring atoms uses a "box" partitioning method. - // Each box has an edge length of sradius, and the number of boxes in each direction is recorded here. + // Each box has an edge length of sradius + 0.1, and the number of boxes in each direction is recorded here. double box_edge_length=0.0; int box_nx=0; int box_ny=0; 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 4cd27cd8772..65843000981 100644 --- a/source/source_cell/module_neighbor/test/sltk_grid_test.cpp +++ b/source/source_cell/module_neighbor/test/sltk_grid_test.cpp @@ -84,18 +84,26 @@ using SltkGridDeathTest = SltkGridTest; namespace { -using NeighborKey = std::tuple; +using NeighborPairKey = std::tuple; -std::vector collect_neighbor_keys(const Grid& grid) +std::vector collect_neighbor_pair_keys(const Grid& grid) { - std::vector keys; - for (const auto& type_neighbors: grid.all_adj_info) + std::vector keys; + for (int center_type = 0; center_type < static_cast(grid.all_adj_info.size()); ++center_type) { - for (const auto& atom_neighbors: type_neighbors) + 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) { - keys.push_back(NeighborKey(atom->type, atom->natom, atom->cell_x, atom->cell_y, atom->cell_z)); + keys.push_back(NeighborPairKey(center_type, + center_natom, + atom->type, + atom->natom, + atom->cell_x, + atom->cell_y, + atom->cell_z)); } } } @@ -172,7 +180,7 @@ TEST_F(SltkGridTest, OpenMPThreadCountKeepsNeighborSet) grid_many.init(ofs_many, *ucell, radius, pbc); ofs_many.close(); - EXPECT_EQ(collect_neighbor_keys(grid_many), collect_neighbor_keys(grid_one)); + EXPECT_EQ(collect_neighbor_pair_keys(grid_many), collect_neighbor_pair_keys(grid_one)); remove("test_one.out"); remove("test_many.out"); From 74a799b3c9b66aad9ad43de7476e8cb03afa6383 Mon Sep 17 00:00:00 2001 From: Yanglijf Date: Sat, 6 Jun 2026 19:31:34 +0800 Subject: [PATCH 13/15] 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 14/15] 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); } } } From e470c332058c97cf92ef0ec21490efd39c23c10d Mon Sep 17 00:00:00 2001 From: Yanglijf Date: Sun, 7 Jun 2026 01:13:40 +0800 Subject: [PATCH 15/15] Link AtomPack in MD tests --- source/source_md/test/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/source/source_md/test/CMakeLists.txt b/source/source_md/test/CMakeLists.txt index 296e9a10b7c..5a9a5e881ec 100644 --- a/source/source_md/test/CMakeLists.txt +++ b/source/source_md/test/CMakeLists.txt @@ -43,6 +43,7 @@ list(APPEND depend_files ../../source_base/math_integral.cpp ../../source_cell/module_neighbor/sltk_atom_arrange.cpp ../../source_cell/module_neighbor/sltk_atom.cpp + ../../source_cell/module_neighbor/atom_pack.cpp ../../source_cell/module_neighbor/sltk_grid.cpp ../../source_cell/module_neighbor/sltk_grid_driver.cpp ../../source_io/module_output/output.cpp