From a18f1675e0a8982084e16875a39bc749475f7705 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Fri, 16 Jan 2026 17:08:08 +0800 Subject: [PATCH 1/8] Optimize install_openmpi.sh script --- toolchain/scripts/stage1/install_openmpi.sh | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/toolchain/scripts/stage1/install_openmpi.sh b/toolchain/scripts/stage1/install_openmpi.sh index 6656fb9a971..6798e8ef178 100755 --- a/toolchain/scripts/stage1/install_openmpi.sh +++ b/toolchain/scripts/stage1/install_openmpi.sh @@ -91,7 +91,12 @@ case "${with_openmpi}" in # else # EXTRA_CONFIGURE_FLAGS="" # fi - ./configure CFLAGS="${CFLAGS}" \ + ./configure \ + CC=gcc \ + CXX=g++ \ + FC=gfortran \ + F77=gfortran \ + CFLAGS="${CFLAGS}" \ --prefix=${pkg_install_dir} \ --libdir="${pkg_install_dir}/lib" \ --with-libevent=internal \ From 47f8a8aaa17c3e1ecbbfa3396223fb24a1941b23 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Wed, 15 Apr 2026 13:59:05 +0800 Subject: [PATCH 2/8] change threshold in 01_PW/208_PW_CG_float --- tests/01_PW/208_PW_CG_float/threshold | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/01_PW/208_PW_CG_float/threshold b/tests/01_PW/208_PW_CG_float/threshold index cc0ad91b672..65774695dee 100644 --- a/tests/01_PW/208_PW_CG_float/threshold +++ b/tests/01_PW/208_PW_CG_float/threshold @@ -1,5 +1,5 @@ # The float type possesses different precision compared to the double type. # This integration aims to test the functionality of the float type # within the plane-wave (pw) basis -threshold 0.00001 +threshold 0.00002 fatal_threshold 1 From 1eb7b268d6b6e16e38a7d8ab5be41a4d1dd638b8 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Sun, 19 Apr 2026 18:37:04 +0800 Subject: [PATCH 3/8] module_neighbor_search --- CMakeLists.txt | 1 + source/Makefile.Objects | 6 + source/source_cell/CMakeLists.txt | 1 + .../module_neighbor_search/CMakeLists.txt | 11 + .../module_neighbor_search/bin_manager.cpp | 190 ++++++++++++++++ .../module_neighbor_search/bin_manager.h | 67 ++++++ .../module_neighbor_search/neighbor_atom.h | 33 +++ .../module_neighbor_search/neighbor_list.h | 95 ++++++++ .../neighbor_search.cpp | 211 ++++++++++++++++++ .../module_neighbor_search/neighbor_search.h | 57 +++++ source/source_esolver/esolver_lj.cpp | 47 +++- 11 files changed, 717 insertions(+), 2 deletions(-) create mode 100644 source/source_cell/module_neighbor_search/CMakeLists.txt create mode 100644 source/source_cell/module_neighbor_search/bin_manager.cpp create mode 100644 source/source_cell/module_neighbor_search/bin_manager.h create mode 100644 source/source_cell/module_neighbor_search/neighbor_atom.h create mode 100644 source/source_cell/module_neighbor_search/neighbor_list.h create mode 100644 source/source_cell/module_neighbor_search/neighbor_search.cpp create mode 100644 source/source_cell/module_neighbor_search/neighbor_search.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 6dc0bf4b5f8..37b8f7e1317 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -762,6 +762,7 @@ target_link_libraries( planewave surchem neighbor + neighbor_search io_input io_basic io_advanced diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 9ad0ebe5094..1e2ba00a7c1 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -16,6 +16,7 @@ VPATH=./src_global:\ ./source_basis/module_ao:\ ./source_basis/module_nao:\ ./source_cell/module_neighbor:\ +./source_cell/module_neighbor_search:\ ./source_cell/module_symmetry:\ ./source_cell:\ ./source_base:\ @@ -86,6 +87,7 @@ ${OBJS_HAMILT}\ ${OBJS_HSOLVER}\ ${OBJS_MD}\ ${OBJS_NEIGHBOR}\ +${OBJS_NEIGHBOR_SEARCH}\ ${OBJS_PSI}\ ${OBJS_PSI_INITIALIZER}\ ${OBJS_PW}\ @@ -411,6 +413,10 @@ OBJS_NEIGHBOR=sltk_atom.o\ sltk_grid.o\ sltk_grid_driver.o\ +OBJS_NEIGHBOR_SEARCH=neighbor_search.o\ + bin_manager.o\ + + OBJS_ORBITAL=ORB_atomic.o\ ORB_atomic_lm.o\ ORB_gaunt_table.o\ diff --git a/source/source_cell/CMakeLists.txt b/source/source_cell/CMakeLists.txt index 39cfa88c22c..69b4b181333 100644 --- a/source/source_cell/CMakeLists.txt +++ b/source/source_cell/CMakeLists.txt @@ -1,5 +1,6 @@ add_subdirectory(module_symmetry) add_subdirectory(module_neighbor) +add_subdirectory(module_neighbor_search) add_library( cell diff --git a/source/source_cell/module_neighbor_search/CMakeLists.txt b/source/source_cell/module_neighbor_search/CMakeLists.txt new file mode 100644 index 00000000000..173d3463cd6 --- /dev/null +++ b/source/source_cell/module_neighbor_search/CMakeLists.txt @@ -0,0 +1,11 @@ +add_library( + neighbor_search + OBJECT + bin_manager.cpp + neighbor_search.cpp +) + +if(ENABLE_COVERAGE) + add_coverage(neighbor_search) +endif() + diff --git a/source/source_cell/module_neighbor_search/bin_manager.cpp b/source/source_cell/module_neighbor_search/bin_manager.cpp new file mode 100644 index 00000000000..66ed75c151d --- /dev/null +++ b/source/source_cell/module_neighbor_search/bin_manager.cpp @@ -0,0 +1,190 @@ +#include +#include +#include +#include "bin_manager.h" + + + +void BinManager::init_bins( + double sr, + const std::vector& inside_atoms, + const std::vector& ghost_atoms +) +{ + sradius = sr; + + x_min = y_min = z_min = std::numeric_limits::max(); + + x_max = y_max = z_max = std::numeric_limits::lowest(); + + auto update_bounds = [&](const std::vector& atoms) + { + for (const auto& atom : atoms) + { + x_min = std::min(x_min, atom.position_x); + x_max = std::max(x_max, atom.position_x); + + y_min = std::min(y_min, atom.position_y); + y_max = std::max(y_max, atom.position_y); + + z_min = std::min(z_min, atom.position_z); + z_max = std::max(z_max, atom.position_z); + } + }; + + update_bounds(inside_atoms); + update_bounds(ghost_atoms); + + bin_sizex = bin_sizey = bin_sizez = sradius; + + nbinx = std::ceil((x_max - x_min) / bin_sizex); + nbiny = std::ceil((y_max - y_min) / bin_sizey); + nbinz = std::ceil((z_max - z_min) / bin_sizez); + + nbinx = std::max(1, nbinx); + nbiny = std::max(1, nbiny); + nbinz = std::max(1, nbinz); + + nbins = nbinx * nbiny * nbinz; + + bins.clear(); + + bins.resize(nbins); + + for (int ix = 0; ix < nbinx; ++ix) + { + for (int iy = 0; iy < nbiny; ++iy) + { + for (int iz = 0; iz < nbinz; ++iz) + { + int idx = ix * nbiny * nbinz + iy * nbinz + iz; + + bins[idx].id_x = ix; + bins[idx].id_y = iy; + bins[idx].id_z = iz; + + bins[idx].atoms.clear(); + } + } + } +} + +void BinManager::do_binning( + const std::vector& inside_atoms, + const std::vector& ghost_atoms +) +{ + auto bin_atom = [&](const NeighborAtom& atom) + { + int ix = std::min( + std::max(int((atom.position_x - x_min) / bin_sizex), 0), + nbinx - 1 + ); + + int iy = std::min( + std::max(int((atom.position_y - y_min) / bin_sizey), 0), + nbiny - 1 + ); + + int iz = std::min( + std::max(int((atom.position_z - z_min) / bin_sizez), 0), + nbinz - 1 + ); + + int idx = ix * nbiny * nbinz + iy * nbinz + iz; + + bins[idx].atoms.push_back(atom); + }; + + for (const auto& atom : inside_atoms) bin_atom(atom); + + for (const auto& atom : ghost_atoms) bin_atom(atom); +} + +void BinManager::build_atom_neighbors( + NeighborList& neighbor_list, + std::vector& atoms +) +{ + assert(atoms.size() == neighbor_list.numneigh.size()); + + double sradius2 = sradius * sradius; + + neighbor_list.reset(); + + for (int i = 0; i < atoms.size(); i++) + { + std::vector neigh_tmp; + + int ix = std::min( + std::max(int((atoms[i].position_x - x_min) / bin_sizex), 0), + nbinx - 1 + ); + + int iy = std::min( + std::max(int((atoms[i].position_y - y_min) / bin_sizey), 0), + nbiny - 1 + ); + + int iz = std::min( + std::max(int((atoms[i].position_z - z_min) / bin_sizez), 0), + nbinz - 1 + ); + + for (int dx = -1; dx <= 1; dx++) + { + for (int dy = -1; dy <= 1; dy++) + { + for (int dz = -1; dz <= 1; dz++) + { + int jx = ix + dx; + int jy = iy + dy; + int jz = iz + dz; + + if (jx < 0 || jx >= nbinx || + jy < 0 || jy >= nbiny || + jz < 0 || jz >= nbinz) + continue; + + int nidx = jx * nbiny * nbinz + jy * nbinz + jz; + + for (auto natom : bins[nidx].atoms) + { + double dx = atoms[i].position_x - natom.position_x; + double dy = atoms[i].position_y - natom.position_y; + double dz = atoms[i].position_z - natom.position_z; + + double dist2 = dx * dx + dy * dy + dz * dz; + + if (dist2 <= sradius2 && dist2 != 0) + { + neigh_tmp.push_back(natom.atom_id); + } + } + } + } + } + + int n = neigh_tmp.size(); + + int* ptr = neighbor_list.allocator.allocate(n); + + for (int k = 0; k < n; k++) + { + ptr[k] = neigh_tmp[k]; + } + + neighbor_list.firstneigh[i] = ptr; + neighbor_list.numneigh[i] = n; + } +} + +void BinManager::clear() +{ + for (auto& bin : bins) + { + bin.atoms.clear(); + } + + bins.clear(); +} \ No newline at end of file diff --git a/source/source_cell/module_neighbor_search/bin_manager.h b/source/source_cell/module_neighbor_search/bin_manager.h new file mode 100644 index 00000000000..bceefd5d76c --- /dev/null +++ b/source/source_cell/module_neighbor_search/bin_manager.h @@ -0,0 +1,67 @@ +#ifndef BIN_MANAGER_H +#define BIN_MANAGER_H + +#include +#include "source_cell/module_neighbor_search/neighbor_atom.h" +#include "source_cell/module_neighbor_search/neighbor_list.h" + + + +class Bin +{ +public: + Bin()=default; + ~Bin()=default; + + int id_x; + int id_y; + int id_z; + + std::vector atoms; +}; + + +class BinManager +{ +public: + + void init_bins( + double sr, + const std::vector& inside_atoms, + const std::vector& ghost_atoms + ); + + + void do_binning( + const std::vector& inside_atoms, + const std::vector& ghost_atoms + ); + + + void build_atom_neighbors( + NeighborList& neighbor_list, + std::vector& atoms + ); + + + void clear(); + + + double sradius; + + double x_min, y_min, z_min; + double x_max, y_max, z_max; + + double bin_sizex; + double bin_sizey; + double bin_sizez; + + int nbinx; + int nbiny; + int nbinz; + int nbins; + + std::vector bins; +}; + +#endif // BIN_MANAGER_H \ No newline at end of file diff --git a/source/source_cell/module_neighbor_search/neighbor_atom.h b/source/source_cell/module_neighbor_search/neighbor_atom.h new file mode 100644 index 00000000000..efc6a4af523 --- /dev/null +++ b/source/source_cell/module_neighbor_search/neighbor_atom.h @@ -0,0 +1,33 @@ +#ifndef NEIGHBOR_ATOM_H +#define NEIGHBOR_ATOM_H + +#include + +class NeighborAtom +{ +public: + double position_x; + double position_y; + double position_z; + int atom_type; + int atom_index; + int atom_id; + bool isghost; + + NeighborAtom(double x, double y, double z, int type, int index, int id) + : position_x(x), position_y(y), position_z(z), + atom_type(type), atom_index(index), atom_id(id) {} +}; + +class InputAtoms +{ +public: + std::vector InputAtom; + double x_low, x_high, y_low, y_high, z_low, z_high; + int n_atoms; + + InputAtoms() + : x_low(0), x_high(0), y_low(0), y_high(0), z_low(0), z_high(0), n_atoms(0) {} +}; + +#endif // NEIGHBOR_ATOM_H \ No newline at end of file diff --git a/source/source_cell/module_neighbor_search/neighbor_list.h b/source/source_cell/module_neighbor_search/neighbor_list.h new file mode 100644 index 00000000000..c54e6ff3f58 --- /dev/null +++ b/source/source_cell/module_neighbor_search/neighbor_list.h @@ -0,0 +1,95 @@ +#pragma once + +#include +#include +#include +#include + +class PageAllocator { +public: + struct Page { + std::vector data; + int capacity; + int offset; + }; + + std::vector pages; + int pgsize; + + static constexpr int DEFAULT_PGSIZE = 1024; + + // Default constructor + PageAllocator() : pgsize(DEFAULT_PGSIZE) { + if (pgsize > 0) new_page(); + } + + PageAllocator(int pgsize_) : pgsize(pgsize_) { + new_page(); + } + + ~PageAllocator() { + // no manual delete[]; vectors clean themselves + } + + PageAllocator(const PageAllocator&) = delete; + PageAllocator& operator=(const PageAllocator&) = delete; + + // Allow move + PageAllocator(PageAllocator&&) = default; + PageAllocator& operator=(PageAllocator&&) = default; + + void new_page() { + Page p; + p.capacity = pgsize; + p.offset = 0; + p.data.resize(pgsize); + pages.push_back(std::move(p)); + } + + int* allocate(int n) { + if (n <= 0) return nullptr; + if (pages.empty()) new_page(); + Page& p = pages.back(); + if (p.offset + n > p.capacity) { + new_page(); + return allocate(n); + } + int* ptr = p.data.data() + p.offset; + p.offset += n; + return ptr; + } + + void reset() { + for (auto& p : pages) { + p.offset = 0; + } + } +}; + +////////////////////////////////////////////////////////////// +// Neighbor List +////////////////////////////////////////////////////////////// + +class NeighborList { +public: + NeighborList() = default; + ~NeighborList() = default; + + int nlocal; + + std::vector numneigh; + std::vector firstneigh; + + PageAllocator allocator; + + void initialize(int n, int pgsize) { + nlocal = n; + allocator = PageAllocator(pgsize); + numneigh.resize(n); + firstneigh.resize(n); + } + + void reset() { + allocator.reset(); + } +}; \ No newline at end of file diff --git a/source/source_cell/module_neighbor_search/neighbor_search.cpp b/source/source_cell/module_neighbor_search/neighbor_search.cpp new file mode 100644 index 00000000000..3b157858769 --- /dev/null +++ b/source/source_cell/module_neighbor_search/neighbor_search.cpp @@ -0,0 +1,211 @@ +#include "source_cell/module_neighbor_search/neighbor_search.h" +#include +#include +#include + + +InputAtoms NeighborSearch::ucell_to_input_atoms(const UnitCell& ucell) +{ + InputAtoms input_atoms; + int atom_count = 0; + + input_atoms.x_low = ucell.atoms[0].tau[0].x; + input_atoms.x_high = ucell.atoms[0].tau[0].x; + input_atoms.y_low = ucell.atoms[0].tau[0].y; + input_atoms.y_high = ucell.atoms[0].tau[0].y; + input_atoms.z_low = ucell.atoms[0].tau[0].z; + input_atoms.z_high = ucell.atoms[0].tau[0].z; + + for (int i = 0; i < ucell.ntype; i++) + { + for (int j = 0; j < ucell.atoms[i].na; j++) + { + NeighborAtom atom( + ucell.atoms[i].tau[j].x, + ucell.atoms[i].tau[j].y, + ucell.atoms[i].tau[j].z, + i, + j, + atom_count + ); + input_atoms.InputAtom.push_back(atom); + + input_atoms.x_low = std::min(input_atoms.x_low, atom.position_x); + input_atoms.x_high = std::max(input_atoms.x_high, atom.position_x); + input_atoms.y_low = std::min(input_atoms.y_low, atom.position_y); + input_atoms.y_high = std::max(input_atoms.y_high, atom.position_y); + input_atoms.z_low = std::min(input_atoms.z_low, atom.position_z); + input_atoms.z_high = std::max(input_atoms.z_high, atom.position_z); + + atom_count++; + } + } + + input_atoms.n_atoms = atom_count; + return input_atoms; +} + +void NeighborSearch::init(const UnitCell& ucell, double sr, int mpi_rank) +{ + search_radius = sr / ucell.lat0; + Check_Expand_Condition(ucell); + setMemberVariables(ucell); + InputAtoms atoms = ucell_to_input_atoms(ucell); + + int mpi_size = 1; + int nx, ny, nz; + decompose(mpi_size, nx, ny, nz); + + z = mpi_rank / (nx * ny); + y = (mpi_rank % (nx * ny)) / ny; + x = mpi_rank % (nx * ny) % ny; + + wide_x = (atoms.x_high - atoms.x_low) / nx; + wide_y = (atoms.y_high - atoms.y_low) / ny; + wide_z = (atoms.z_high - atoms.z_low) / nz; + + for (int i = 0; i < all_atoms.size(); i++) + { + bool in_x = (all_atoms[i].position_x >= atoms.x_low + x * wide_x) && + (all_atoms[i].position_x <= atoms.x_low + (x + 1) * wide_x); + bool in_y = (all_atoms[i].position_y >= atoms.y_low + y * wide_y) && + (all_atoms[i].position_y <= atoms.y_low + (y + 1) * wide_y); + bool in_z = (all_atoms[i].position_z >= atoms.z_low + z * wide_z) && + (all_atoms[i].position_z <= atoms.z_low + (z + 1) * wide_z); + + if (in_x && in_y && in_z) + { + all_atoms[i].isghost = false; + inside_atoms.push_back(all_atoms[i]); + } + else if (distance( + all_atoms[i].position_x, + all_atoms[i].position_y, + all_atoms[i].position_z, + atoms.x_low, + atoms.y_low, + atoms.z_low) <= search_radius * search_radius) + { + all_atoms[i].isghost = true; + ghost_atoms.push_back(all_atoms[i]); + } + } + + neighbor_list.initialize(ucell.nat, 10000000); +} + +void NeighborSearch::build_neighbors() +{ + bin_manager.init_bins(search_radius, inside_atoms, ghost_atoms); + bin_manager.do_binning(inside_atoms, ghost_atoms); + bin_manager.build_atom_neighbors(neighbor_list, inside_atoms); +} + +void NeighborSearch::Check_Expand_Condition(const UnitCell& ucell) +{ + double a23_1 = ucell.latvec.e22 * ucell.latvec.e33 - ucell.latvec.e23 * ucell.latvec.e32; + double a23_2 = ucell.latvec.e21 * ucell.latvec.e33 - ucell.latvec.e23 * ucell.latvec.e31; + double a23_3 = ucell.latvec.e21 * ucell.latvec.e32 - ucell.latvec.e22 * ucell.latvec.e31; + double a23_norm = sqrt(a23_1 * a23_1 + a23_2 * a23_2 + a23_3 * a23_3); + double extend_v = a23_norm * search_radius; + double extend_d1 = extend_v / ucell.omega * ucell.lat0 * ucell.lat0 * ucell.lat0; + int extend_d11 = std::ceil(extend_d1); + + double a31_1 = ucell.latvec.e32 * ucell.latvec.e13 - ucell.latvec.e33 * ucell.latvec.e12; + double a31_2 = ucell.latvec.e31 * ucell.latvec.e13 - ucell.latvec.e33 * ucell.latvec.e11; + double a31_3 = ucell.latvec.e31 * ucell.latvec.e12 - ucell.latvec.e32 * ucell.latvec.e11; + double a31_norm = sqrt(a31_1 * a31_1 + a31_2 * a31_2 + a31_3 * a31_3); + double extend_d2 = a31_norm * search_radius / ucell.omega * ucell.lat0 * ucell.lat0 * ucell.lat0; + int extend_d22 = std::ceil(extend_d2); + + double a12_1 = ucell.latvec.e12 * ucell.latvec.e23 - ucell.latvec.e13 * ucell.latvec.e22; + double a12_2 = ucell.latvec.e11 * ucell.latvec.e23 - ucell.latvec.e13 * ucell.latvec.e21; + double a12_3 = ucell.latvec.e11 * ucell.latvec.e22 - ucell.latvec.e12 * ucell.latvec.e21; + double a12_norm = sqrt(a12_1 * a12_1 + a12_2 * a12_2 + a12_3 * a12_3); + double extend_d3 = a12_norm * search_radius / ucell.omega * ucell.lat0 * ucell.lat0 * ucell.lat0; + int extend_d33 = std::ceil(extend_d3); + + glayerX = extend_d11 + 1; + glayerY = extend_d22 + 1; + glayerZ = extend_d33 + 1; + glayerX_minus = extend_d11; + glayerY_minus = extend_d22; + glayerZ_minus = extend_d33; +} + +void NeighborSearch::setMemberVariables(const UnitCell& ucell) +{ + all_atoms.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); + + int atom_count = 0; + + for (int ix = -glayerX_minus; ix < glayerX; ix++) + { + for (int iy = -glayerY_minus; iy < glayerY; iy++) + { + for (int iz = -glayerZ_minus; iz < glayerZ; iz++) + { + for (int i = 0; i < ucell.ntype; i++) + { + for (int j = 0; j < ucell.atoms[i].na; j++) + { + double x = ucell.atoms[i].tau[j].x + vec1[0] * ix + vec2[0] * iy + vec3[0] * iz; + double y = ucell.atoms[i].tau[j].y + vec1[1] * ix + vec2[1] * iy + vec3[1] * iz; + double z = ucell.atoms[i].tau[j].z + vec1[2] * ix + vec2[2] * iy + vec3[2] * iz; + + NeighborAtom atom(x, y, z, i, j, atom_count); + all_atoms.push_back(atom); + atom_count++; + } + } + } + } + } +} + +double NeighborSearch::distance( + double position_x, + double position_y, + double position_z, + double x_low, + double y_low, + double z_low) +{ + double dx = std::max(0.0, std::max(x_low + x * wide_x - position_x, position_x - (x_low + (x + 1) * wide_x))); + double dy = std::max(0.0, std::max(y_low + y * wide_y - position_y, position_y - (y_low + (y + 1) * wide_y))); + double dz = std::max(0.0, std::max(z_low + z * wide_z - position_z, position_z - (z_low + (z + 1) * wide_z))); + return dx * dx + dy * dy + dz * dz; +} + +void NeighborSearch::decompose(int mpi_size, int &nx, int &ny, int &nz) +{ + nx = 1; + ny = 1; + nz = mpi_size; + + int cube = cbrt(mpi_size); + for (int i = cube; i >= 1; i--) + { + if (mpi_size % i == 0) + { + nx = i; + ny = mpi_size / i; + break; + } + } + + int sq = sqrt(ny); + for (int i = sq; i >= 1; i--) + { + if (ny % i == 0) + { + nz = ny / i; + ny = i; + break; + } + } +} \ No newline at end of file diff --git a/source/source_cell/module_neighbor_search/neighbor_search.h b/source/source_cell/module_neighbor_search/neighbor_search.h new file mode 100644 index 00000000000..1a50d3078ec --- /dev/null +++ b/source/source_cell/module_neighbor_search/neighbor_search.h @@ -0,0 +1,57 @@ +#ifndef NEIGHBOR_SEARCH_H +#define NEIGHBOR_SEARCH_H + + +#include "source_cell/module_neighbor_search/neighbor_atom.h" +#include "source_cell/module_neighbor_search/bin_manager.h" +#include "source_cell/module_neighbor_search/neighbor_list.h" +#include "source_cell/unitcell.h" + +class NeighborSearch +{ +public: + NeighborSearch()=default; + ~NeighborSearch()=default; + + void init(const UnitCell& ucell, double sr, int mpi_rank); + + + void build_neighbors(); + + InputAtoms ucell_to_input_atoms(const UnitCell& ucell); + + void Check_Expand_Condition(const UnitCell& ucell); + + void setMemberVariables(const UnitCell& ucell); + + double distance( + double position_x, + double position_y, + double position_z, + double x_low, + double y_low, + double z_low); + + void decompose(int mpi_size, int& nx, int& ny, int& nz); + + double search_radius; + + + int x, y, z; + double wide_x, wide_y, wide_z; + + // 扩展层 + int glayerX, glayerY, glayerZ; + int glayerX_minus, glayerY_minus, glayerZ_minus; + + // 原子集合 + std::vector all_atoms; + std::vector inside_atoms; + std::vector ghost_atoms; + + // 近邻表 & bin 管理 + NeighborList neighbor_list; + BinManager bin_manager; +}; + +#endif \ No newline at end of file diff --git a/source/source_esolver/esolver_lj.cpp b/source/source_esolver/esolver_lj.cpp index f71698586ea..d38cf3d9211 100644 --- a/source/source_esolver/esolver_lj.cpp +++ b/source/source_esolver/esolver_lj.cpp @@ -4,6 +4,7 @@ #include "source_cell/module_neighbor/sltk_grid_driver.h" #include "source_io/module_output/output_log.h" #include "source_io/module_output/cif_io.h" +#include "source_cell/module_neighbor_search/neighbor_search.h" namespace ModuleESolver @@ -32,7 +33,49 @@ void ESolver_LJ::before_all_runners(UnitCell& ucell, const Input_para& inp) void ESolver_LJ::runner(UnitCell& ucell, const int istep) { - Grid_Driver grid_neigh(PARAM.inp.test_deconstructor, PARAM.inp.test_grid); + NeighborSearch neighbor_search; + neighbor_search.init(ucell, search_radius, 0); + neighbor_search.build_neighbors(); + + double distance = 0.0; + int index = 0; + + // Important! potential, force, virial must be zero per step + lj_potential = 0; + lj_force.zero_out(); + lj_virial.zero_out(); + + ModuleBase::Vector3 tau1, tau2, dtau; + for (int it = 0; it < ucell.ntype; ++it) + { + Atom* atom1 = &ucell.atoms[it]; + for (int ia = 0; ia < atom1->na; ++ia) + { + tau1 = atom1->tau[ia]; + for (int ad = 0; ad < neighbor_search.neighbor_list.numneigh[index]; ++ad) + { + tau2.x = neighbor_search.all_atoms[neighbor_search.neighbor_list.firstneigh[index][ad]].position_x; + tau2.y = neighbor_search.all_atoms[neighbor_search.neighbor_list.firstneigh[index][ad]].position_y; + tau2.z = neighbor_search.all_atoms[neighbor_search.neighbor_list.firstneigh[index][ad]].position_z; + int it2 = neighbor_search.all_atoms[neighbor_search.neighbor_list.firstneigh[index][ad]].atom_type; + dtau = (tau1 - tau2) * ucell.lat0; + distance = dtau.norm(); + if (distance < lj_rcut(it, it2)) + { + lj_potential += LJ_energy(distance, it, it2) - en_shift(it, it2); + ModuleBase::Vector3 f_ij = LJ_force(dtau, it, it2); + lj_force(index, 0) += f_ij.x; + lj_force(index, 1) += f_ij.y; + lj_force(index, 2) += f_ij.z; + LJ_virial(f_ij, dtau); + } + } + index++; + } + } + + + /*Grid_Driver grid_neigh(PARAM.inp.test_deconstructor, PARAM.inp.test_grid); atom_arrange::search(PARAM.globalv.search_pbc, GlobalV::ofs_running, grid_neigh, @@ -74,7 +117,7 @@ void ESolver_LJ::runner(UnitCell& ucell, const int istep) } index++; } - } + }*/ lj_potential /= 2.0; GlobalV::ofs_running << " #TOTAL ENERGY# " << std::setprecision(11) << lj_potential * ModuleBase::Ry_to_eV << " eV" From 888c61871da5f3b03411db2ad5bf922618cf5dbd Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Mon, 20 Apr 2026 16:17:25 +0800 Subject: [PATCH 4/8] module_neighbor_search --- .../module_neighbor_search/bin_manager.cpp | 5 +- .../module_neighbor_search/bin_manager.h | 1 - .../neighbor_search.cpp | 94 ++++++++++--------- .../module_neighbor_search/neighbor_search.h | 10 +- .../unitcell_interface.h | 25 +++++ source/source_cell/unitcell.h | 31 +++++- 6 files changed, 111 insertions(+), 55 deletions(-) create mode 100644 source/source_cell/module_neighbor_search/unitcell_interface.h diff --git a/source/source_cell/module_neighbor_search/bin_manager.cpp b/source/source_cell/module_neighbor_search/bin_manager.cpp index 66ed75c151d..cb4ba57f555 100644 --- a/source/source_cell/module_neighbor_search/bin_manager.cpp +++ b/source/source_cell/module_neighbor_search/bin_manager.cpp @@ -45,7 +45,7 @@ void BinManager::init_bins( nbiny = std::max(1, nbiny); nbinz = std::max(1, nbinz); - nbins = nbinx * nbiny * nbinz; + int nbins = nbinx * nbiny * nbinz; bins.clear(); @@ -164,9 +164,8 @@ void BinManager::build_atom_neighbors( } } } - + int n = neigh_tmp.size(); - int* ptr = neighbor_list.allocator.allocate(n); for (int k = 0; k < n; k++) diff --git a/source/source_cell/module_neighbor_search/bin_manager.h b/source/source_cell/module_neighbor_search/bin_manager.h index bceefd5d76c..a2d050bbc04 100644 --- a/source/source_cell/module_neighbor_search/bin_manager.h +++ b/source/source_cell/module_neighbor_search/bin_manager.h @@ -59,7 +59,6 @@ class BinManager int nbinx; int nbiny; int nbinz; - int nbins; std::vector bins; }; diff --git a/source/source_cell/module_neighbor_search/neighbor_search.cpp b/source/source_cell/module_neighbor_search/neighbor_search.cpp index 3b157858769..b6a4c3336d4 100644 --- a/source/source_cell/module_neighbor_search/neighbor_search.cpp +++ b/source/source_cell/module_neighbor_search/neighbor_search.cpp @@ -4,26 +4,22 @@ #include -InputAtoms NeighborSearch::ucell_to_input_atoms(const UnitCell& ucell) +InputAtoms NeighborSearch::ucell_to_input_atoms(const IAtomProvider& ucell) { InputAtoms input_atoms; int atom_count = 0; - input_atoms.x_low = ucell.atoms[0].tau[0].x; - input_atoms.x_high = ucell.atoms[0].tau[0].x; - input_atoms.y_low = ucell.atoms[0].tau[0].y; - input_atoms.y_high = ucell.atoms[0].tau[0].y; - input_atoms.z_low = ucell.atoms[0].tau[0].z; - input_atoms.z_high = ucell.atoms[0].tau[0].z; + input_atoms.x_low = input_atoms.y_low = input_atoms.z_low = std::numeric_limits::max(); + input_atoms.x_high = input_atoms.y_high = input_atoms.z_high = std::numeric_limits::lowest(); - for (int i = 0; i < ucell.ntype; i++) + for (int i = 0; i < ucell.get_ntype(); i++) { - for (int j = 0; j < ucell.atoms[i].na; j++) + for (int j = 0; j < ucell.get_na(i); j++) { NeighborAtom atom( - ucell.atoms[i].tau[j].x, - ucell.atoms[i].tau[j].y, - ucell.atoms[i].tau[j].z, + ucell.get_tauu(i,j).x, + ucell.get_tauu(i,j).y, + ucell.get_tauu(i,j).z, i, j, atom_count @@ -45,9 +41,9 @@ InputAtoms NeighborSearch::ucell_to_input_atoms(const UnitCell& ucell) return input_atoms; } -void NeighborSearch::init(const UnitCell& ucell, double sr, int mpi_rank) +void NeighborSearch::init(const IAtomProvider& ucell, double sr, int mpi_rank) { - search_radius = sr / ucell.lat0; + search_radius = sr / ucell.get_lat0(); Check_Expand_Condition(ucell); setMemberVariables(ucell); InputAtoms atoms = ucell_to_input_atoms(ucell); @@ -66,14 +62,22 @@ void NeighborSearch::init(const UnitCell& ucell, double sr, int mpi_rank) for (int i = 0; i < all_atoms.size(); i++) { - bool in_x = (all_atoms[i].position_x >= atoms.x_low + x * wide_x) && - (all_atoms[i].position_x <= atoms.x_low + (x + 1) * wide_x); - bool in_y = (all_atoms[i].position_y >= atoms.y_low + y * wide_y) && - (all_atoms[i].position_y <= atoms.y_low + (y + 1) * wide_y); - bool in_z = (all_atoms[i].position_z >= atoms.z_low + z * wide_z) && - (all_atoms[i].position_z <= atoms.z_low + (z + 1) * wide_z); - - if (in_x && in_y && in_z) + + int in_x = std::min( + static_cast(std::floor((all_atoms[i].position_x - atoms.x_low) / wide_x)), + nx - 1 + ); + int in_y = std::min( + static_cast(std::floor((all_atoms[i].position_y - atoms.y_low) / wide_y)), + ny - 1 + ); + int in_z = std::min( + static_cast(std::floor((all_atoms[i].position_z - atoms.z_low) / wide_z)), + nz - 1 + ); + + + if (in_x==x && in_y==y && in_z==z&&all_atoms[i].position_x<=atoms.x_high&&all_atoms[i].position_y<=atoms.y_high&&all_atoms[i].position_z<=atoms.z_high) { all_atoms[i].isghost = false; inside_atoms.push_back(all_atoms[i]); @@ -91,7 +95,7 @@ void NeighborSearch::init(const UnitCell& ucell, double sr, int mpi_rank) } } - neighbor_list.initialize(ucell.nat, 10000000); + neighbor_list.initialize(ucell.get_natom(), 100000000); } void NeighborSearch::build_neighbors() @@ -101,28 +105,28 @@ void NeighborSearch::build_neighbors() bin_manager.build_atom_neighbors(neighbor_list, inside_atoms); } -void NeighborSearch::Check_Expand_Condition(const UnitCell& ucell) +void NeighborSearch::Check_Expand_Condition(const IAtomProvider& ucell) { - double a23_1 = ucell.latvec.e22 * ucell.latvec.e33 - ucell.latvec.e23 * ucell.latvec.e32; - double a23_2 = ucell.latvec.e21 * ucell.latvec.e33 - ucell.latvec.e23 * ucell.latvec.e31; - double a23_3 = ucell.latvec.e21 * ucell.latvec.e32 - ucell.latvec.e22 * ucell.latvec.e31; + double a23_1 = ucell.get_latvec().e22 * ucell.get_latvec().e33 - ucell.get_latvec().e23 * ucell.get_latvec().e32; + double a23_2 = ucell.get_latvec().e21 * ucell.get_latvec().e33 - ucell.get_latvec().e23 * ucell.get_latvec().e31; + double a23_3 = ucell.get_latvec().e21 * ucell.get_latvec().e32 - ucell.get_latvec().e22 * ucell.get_latvec().e31; double a23_norm = sqrt(a23_1 * a23_1 + a23_2 * a23_2 + a23_3 * a23_3); double extend_v = a23_norm * search_radius; - double extend_d1 = extend_v / ucell.omega * ucell.lat0 * ucell.lat0 * ucell.lat0; + double extend_d1 = extend_v / ucell.get_omega() * ucell.get_lat0() * ucell.get_lat0() * ucell.get_lat0(); int extend_d11 = std::ceil(extend_d1); - double a31_1 = ucell.latvec.e32 * ucell.latvec.e13 - ucell.latvec.e33 * ucell.latvec.e12; - double a31_2 = ucell.latvec.e31 * ucell.latvec.e13 - ucell.latvec.e33 * ucell.latvec.e11; - double a31_3 = ucell.latvec.e31 * ucell.latvec.e12 - ucell.latvec.e32 * ucell.latvec.e11; + double a31_1 = ucell.get_latvec().e32 * ucell.get_latvec().e13 - ucell.get_latvec().e33 * ucell.get_latvec().e12; + double a31_2 = ucell.get_latvec().e31 * ucell.get_latvec().e13 - ucell.get_latvec().e33 * ucell.get_latvec().e11; + double a31_3 = ucell.get_latvec().e31 * ucell.get_latvec().e12 - ucell.get_latvec().e32 * ucell.get_latvec().e11; double a31_norm = sqrt(a31_1 * a31_1 + a31_2 * a31_2 + a31_3 * a31_3); - double extend_d2 = a31_norm * search_radius / ucell.omega * ucell.lat0 * ucell.lat0 * ucell.lat0; + double extend_d2 = a31_norm * search_radius / ucell.get_omega() * ucell.get_lat0() * ucell.get_lat0() * ucell.get_lat0(); int extend_d22 = std::ceil(extend_d2); - double a12_1 = ucell.latvec.e12 * ucell.latvec.e23 - ucell.latvec.e13 * ucell.latvec.e22; - double a12_2 = ucell.latvec.e11 * ucell.latvec.e23 - ucell.latvec.e13 * ucell.latvec.e21; - double a12_3 = ucell.latvec.e11 * ucell.latvec.e22 - ucell.latvec.e12 * ucell.latvec.e21; + double a12_1 = ucell.get_latvec().e12 * ucell.get_latvec().e23 - ucell.get_latvec().e13 * ucell.get_latvec().e22; + double a12_2 = ucell.get_latvec().e11 * ucell.get_latvec().e23 - ucell.get_latvec().e13 * ucell.get_latvec().e21; + double a12_3 = ucell.get_latvec().e11 * ucell.get_latvec().e22 - ucell.get_latvec().e12 * ucell.get_latvec().e21; double a12_norm = sqrt(a12_1 * a12_1 + a12_2 * a12_2 + a12_3 * a12_3); - double extend_d3 = a12_norm * search_radius / ucell.omega * ucell.lat0 * ucell.lat0 * ucell.lat0; + double extend_d3 = a12_norm * search_radius / ucell.get_omega() * ucell.get_lat0() * ucell.get_lat0() * ucell.get_lat0(); int extend_d33 = std::ceil(extend_d3); glayerX = extend_d11 + 1; @@ -133,13 +137,13 @@ void NeighborSearch::Check_Expand_Condition(const UnitCell& ucell) glayerZ_minus = extend_d33; } -void NeighborSearch::setMemberVariables(const UnitCell& ucell) +void NeighborSearch::setMemberVariables(const IAtomProvider& ucell) { all_atoms.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); + ModuleBase::Vector3 vec1(ucell.get_latvec().e11, ucell.get_latvec().e12, ucell.get_latvec().e13); + ModuleBase::Vector3 vec2(ucell.get_latvec().e21, ucell.get_latvec().e22, ucell.get_latvec().e23); + ModuleBase::Vector3 vec3(ucell.get_latvec().e31, ucell.get_latvec().e32, ucell.get_latvec().e33); int atom_count = 0; @@ -149,13 +153,13 @@ void NeighborSearch::setMemberVariables(const UnitCell& ucell) { for (int iz = -glayerZ_minus; iz < glayerZ; iz++) { - for (int i = 0; i < ucell.ntype; i++) + for (int i = 0; i < ucell.get_ntype(); i++) { - for (int j = 0; j < ucell.atoms[i].na; j++) + for (int j = 0; j < ucell.get_na(i); j++) { - double x = ucell.atoms[i].tau[j].x + vec1[0] * ix + vec2[0] * iy + vec3[0] * iz; - double y = ucell.atoms[i].tau[j].y + vec1[1] * ix + vec2[1] * iy + vec3[1] * iz; - double z = ucell.atoms[i].tau[j].z + vec1[2] * ix + vec2[2] * iy + vec3[2] * iz; + double x = ucell.get_tauu(i,j).x + vec1[0] * ix + vec2[0] * iy + vec3[0] * iz; + double y = ucell.get_tauu(i,j).y + vec1[1] * ix + vec2[1] * iy + vec3[1] * iz; + double z = ucell.get_tauu(i,j).z + vec1[2] * ix + vec2[2] * iy + vec3[2] * iz; NeighborAtom atom(x, y, z, i, j, atom_count); all_atoms.push_back(atom); diff --git a/source/source_cell/module_neighbor_search/neighbor_search.h b/source/source_cell/module_neighbor_search/neighbor_search.h index 1a50d3078ec..1219345836b 100644 --- a/source/source_cell/module_neighbor_search/neighbor_search.h +++ b/source/source_cell/module_neighbor_search/neighbor_search.h @@ -5,7 +5,7 @@ #include "source_cell/module_neighbor_search/neighbor_atom.h" #include "source_cell/module_neighbor_search/bin_manager.h" #include "source_cell/module_neighbor_search/neighbor_list.h" -#include "source_cell/unitcell.h" +#include "source_cell/module_neighbor_search/unitcell_interface.h" class NeighborSearch { @@ -13,16 +13,16 @@ class NeighborSearch NeighborSearch()=default; ~NeighborSearch()=default; - void init(const UnitCell& ucell, double sr, int mpi_rank); + void init(const IAtomProvider& ucell, double sr, int mpi_rank); void build_neighbors(); - InputAtoms ucell_to_input_atoms(const UnitCell& ucell); + InputAtoms ucell_to_input_atoms(const IAtomProvider& ucell); - void Check_Expand_Condition(const UnitCell& ucell); + void Check_Expand_Condition(const IAtomProvider& ucell); - void setMemberVariables(const UnitCell& ucell); + void setMemberVariables(const IAtomProvider& ucell); double distance( double position_x, diff --git a/source/source_cell/module_neighbor_search/unitcell_interface.h b/source/source_cell/module_neighbor_search/unitcell_interface.h new file mode 100644 index 00000000000..25610d471b3 --- /dev/null +++ b/source/source_cell/module_neighbor_search/unitcell_interface.h @@ -0,0 +1,25 @@ +#ifndef UNITCELL_INTERFACE_H +#define UNITCELL_INTERFACE_H + + +#include "source_base/vector3.h" +#include "source_base/matrix3.h" + + +class IAtomProvider { +public: + virtual ~IAtomProvider() = default; + + virtual double get_lat0() const = 0; + virtual double get_omega() const = 0; + virtual const ModuleBase::Matrix3& get_latvec() const = 0; + + + + virtual int get_natom() const = 0; + virtual int get_na(int i) const = 0; + virtual int get_ntype() const = 0; + virtual ModuleBase::Vector3 get_tauu(int i,int j) const = 0; +}; + +#endif // UNITCELL_INTERFACE_H \ No newline at end of file diff --git a/source/source_cell/unitcell.h b/source/source_cell/unitcell.h index 3292a40a065..959735ceca6 100644 --- a/source/source_cell/unitcell.h +++ b/source/source_cell/unitcell.h @@ -7,6 +7,7 @@ #include "source_estate/magnetism.h" #include "source_io/module_output/output.h" #include "module_symmetry/symmetry.h" +#include "source_cell/module_neighbor_search/unitcell_interface.h" #ifdef __LCAO #include "source_basis/module_ao/ORB_read.h" @@ -14,8 +15,36 @@ #endif // provide the basic information about unitcell. -class UnitCell { +class UnitCell : public IAtomProvider { public: + double get_lat0() const override { + return lat0; + } + + double get_omega() const override { + return omega; + } + + const ModuleBase::Matrix3& get_latvec() const override { + return latvec; + } + + int get_natom() const override { + return nat; + } + + int get_na(int i) const override { + return atoms[i].na; + } + + int get_ntype() const override { + return ntype; + } + + ModuleBase::Vector3 get_tauu(int i, int j) const override { + return atoms[i].tau[j]; + } + Atom* atoms = nullptr; Sep_Cell sep_cell; From b98c44f463ff52068e233cb96845c397e0fa5bf1 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Tue, 21 Apr 2026 17:11:34 +0800 Subject: [PATCH 5/8] module_neighbor_search --- .../module_neighbor_search/bin_manager.cpp | 6 ++- .../module_neighbor_search/unitcell_plus.h | 52 +++++++++++++++++++ source/source_esolver/esolver_lj.cpp | 26 +++++++++- source/source_esolver/esolver_lj.h | 3 ++ 4 files changed, 84 insertions(+), 3 deletions(-) create mode 100644 source/source_cell/module_neighbor_search/unitcell_plus.h diff --git a/source/source_cell/module_neighbor_search/bin_manager.cpp b/source/source_cell/module_neighbor_search/bin_manager.cpp index cb4ba57f555..b916c90a8ec 100644 --- a/source/source_cell/module_neighbor_search/bin_manager.cpp +++ b/source/source_cell/module_neighbor_search/bin_manager.cpp @@ -163,9 +163,11 @@ void BinManager::build_atom_neighbors( } } } - } - + } int n = neigh_tmp.size(); + + //std::cout< + +class UnitCellPlus : public IAtomProvider +{ +public: + UnitCellPlus()=default; + ~UnitCellPlus()=default; + + double get_lat0() const override { + return lat0; + } + + double get_omega() const override { + return omega; + } + + const ModuleBase::Matrix3& get_latvec() const override { + return latvec; + } + + int get_natom() const override { + return nat; + } + + int get_na(int i) const override { + return na[i]; + } + + int get_ntype() const override { + return ntype; + } + + ModuleBase::Vector3 get_tauu(int i, int j) const override { + int count=0; + for(int k=0;k na; + int ntype; + ModuleBase::Matrix3 latvec; + std::vector> tau; + +}; \ No newline at end of file diff --git a/source/source_esolver/esolver_lj.cpp b/source/source_esolver/esolver_lj.cpp index d38cf3d9211..bc885f19b8d 100644 --- a/source/source_esolver/esolver_lj.cpp +++ b/source/source_esolver/esolver_lj.cpp @@ -7,9 +7,32 @@ #include "source_cell/module_neighbor_search/neighbor_search.h" + namespace ModuleESolver { + UnitCellPlus ESolver_LJ::change_from_ucell_to_ucell_plus(UnitCell& ucell) + { + UnitCellPlus ucell_plus; + ucell_plus.lat0 = ucell.lat0; + ucell_plus.omega = ucell.omega; + ucell_plus.nat = ucell.nat; + for(int i=0;i Date: Sun, 26 Apr 2026 23:24:16 +0800 Subject: [PATCH 6/8] module_neighbor_search --- .../module_neighbor_search/CMakeLists.txt | 5 + .../module_neighbor_search/neighbor_atom.h | 3 +- .../module_neighbor_search/neighbor_list.h | 7 +- .../neighbor_search.cpp | 89 ++++++-- .../module_neighbor_search/neighbor_search.h | 1 + .../test/CMakeLists.txt | 16 ++ .../test/test_neighbor_search.cpp | 196 ++++++++++++++++++ 7 files changed, 294 insertions(+), 23 deletions(-) create mode 100644 source/source_cell/module_neighbor_search/test/CMakeLists.txt create mode 100644 source/source_cell/module_neighbor_search/test/test_neighbor_search.cpp diff --git a/source/source_cell/module_neighbor_search/CMakeLists.txt b/source/source_cell/module_neighbor_search/CMakeLists.txt index 173d3463cd6..6f54f71513f 100644 --- a/source/source_cell/module_neighbor_search/CMakeLists.txt +++ b/source/source_cell/module_neighbor_search/CMakeLists.txt @@ -9,3 +9,8 @@ if(ENABLE_COVERAGE) add_coverage(neighbor_search) endif() +if(BUILD_TESTING) + if(ENABLE_MPI) + add_subdirectory(test) + endif() +endif() \ No newline at end of file diff --git a/source/source_cell/module_neighbor_search/neighbor_atom.h b/source/source_cell/module_neighbor_search/neighbor_atom.h index efc6a4af523..ed8d99c161a 100644 --- a/source/source_cell/module_neighbor_search/neighbor_atom.h +++ b/source/source_cell/module_neighbor_search/neighbor_atom.h @@ -12,7 +12,8 @@ class NeighborAtom int atom_type; int atom_index; int atom_id; - bool isghost; + //bool isghost; + bool is_inside; NeighborAtom(double x, double y, double z, int type, int index, int id) : position_x(x), position_y(y), position_z(z), diff --git a/source/source_cell/module_neighbor_search/neighbor_list.h b/source/source_cell/module_neighbor_search/neighbor_list.h index c54e6ff3f58..0db8229ff3e 100644 --- a/source/source_cell/module_neighbor_search/neighbor_list.h +++ b/source/source_cell/module_neighbor_search/neighbor_list.h @@ -84,9 +84,10 @@ class NeighborList { void initialize(int n, int pgsize) { nlocal = n; - allocator = PageAllocator(pgsize); - numneigh.resize(n); - firstneigh.resize(n); + allocator = PageAllocator(pgsize); + // ensure neighbor containers are sized and initialized + numneigh.assign(n, 0); + firstneigh.assign(n, nullptr); } void reset() { diff --git a/source/source_cell/module_neighbor_search/neighbor_search.cpp b/source/source_cell/module_neighbor_search/neighbor_search.cpp index b6a4c3336d4..9775df965d8 100644 --- a/source/source_cell/module_neighbor_search/neighbor_search.cpp +++ b/source/source_cell/module_neighbor_search/neighbor_search.cpp @@ -60,26 +60,69 @@ void NeighborSearch::init(const IAtomProvider& ucell, double sr, int mpi_rank) wide_y = (atoms.y_high - atoms.y_low) / ny; wide_z = (atoms.z_high - atoms.z_low) / nz; + int in_x, in_y, in_z; + for (int i = 0; i < all_atoms.size(); i++) { - - int in_x = std::min( - static_cast(std::floor((all_atoms[i].position_x - atoms.x_low) / wide_x)), - nx - 1 - ); - int in_y = std::min( - static_cast(std::floor((all_atoms[i].position_y - atoms.y_low) / wide_y)), - ny - 1 - ); - int in_z = std::min( - static_cast(std::floor((all_atoms[i].position_z - atoms.z_low) / wide_z)), - nz - 1 - ); - - - if (in_x==x && in_y==y && in_z==z&&all_atoms[i].position_x<=atoms.x_high&&all_atoms[i].position_y<=atoms.y_high&&all_atoms[i].position_z<=atoms.z_high) + if(wide_x==0) + { + if(all_atoms[i].position_x==atoms.x_low) + { + in_x = x; + } + else + { + in_x = std::numeric_limits::max(); + } + } + else + { + in_x = std::min( + static_cast(std::floor((all_atoms[i].position_x - atoms.x_low) / wide_x)), + nx - 1 + ); + } + if(wide_y==0) + { + if(all_atoms[i].position_y==atoms.y_low) + { + in_y = y; + } + else + { + in_y = std::numeric_limits::max(); + } + } + else + { + in_y = std::min( + static_cast(std::floor((all_atoms[i].position_y - atoms.y_low) / wide_y)), + ny - 1 + ); + } + if(wide_z==0) + { + if(all_atoms[i].position_z==atoms.z_low) + { + in_z = z; + } + else + { + in_z = std::numeric_limits::max(); + } + } + else + { + in_z = std::min( + static_cast(std::floor((all_atoms[i].position_z - atoms.z_low) / wide_z)), + nz - 1 + ); + } + //std::cout< +#include "../neighbor_search.h" +#include "../unitcell_plus.h" +#include "../bin_manager.h" +#include "../neighbor_list.h" + +TEST(NeighborSearchTest, TwoAtomsNeighbor) +{ + UnitCellPlus ucell; + + // ===== 基本晶胞信息 ===== + ucell.lat0 = 1.0; + ucell.omega = 1.0; + + ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0; + ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0; + ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1; + + // ===== 原子信息 ===== + ucell.ntype = 1; + ucell.na = {2}; + ucell.nat = 2; + + ucell.tau = { + {0.0, 0.0, 0.0}, + {0.5, 0.0, 0.0} + }; + + NeighborSearch ns; + + double cutoff = 1.0; + + ns.init(ucell, cutoff, 0); + ns.build_neighbors(); + + auto &list = ns.get_neighbor_list(); + + ASSERT_EQ(list.numneigh.size(), 2); + + EXPECT_EQ(list.numneigh[0], 8); + EXPECT_EQ(list.numneigh[1], 8); +} + +TEST(NeighborSearchUnit, UCellToInputAtoms) +{ + UnitCellPlus ucell; + ucell.lat0 = 1.0; + ucell.omega = 1.0; + ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0; + ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0; + ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1; + + ucell.ntype = 1; + ucell.na = {2}; + ucell.nat = 2; + ucell.tau = { + {0.0, 0.0, 0.0}, + {0.5, 0.0, 0.0} + }; + + NeighborSearch ns; + auto inputs = ns.ucell_to_input_atoms(ucell); + + EXPECT_EQ(inputs.n_atoms, 2); + ASSERT_EQ(inputs.InputAtom.size(), 2); + EXPECT_DOUBLE_EQ(inputs.InputAtom[0].position_x, 0.0); + EXPECT_DOUBLE_EQ(inputs.InputAtom[1].position_x, 0.5); + EXPECT_DOUBLE_EQ(inputs.x_low, 0.0); + EXPECT_DOUBLE_EQ(inputs.x_high, 0.5); +} + +TEST(NeighborSearchUnit, CheckExpandAndSetMembers) +{ + UnitCellPlus ucell; + ucell.lat0 = 1.0; + ucell.omega = 1.0; + ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0; + ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0; + ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1; + + ucell.ntype = 1; + ucell.na = {2}; + ucell.nat = 2; + ucell.tau = { + {0.0, 0.0, 0.0}, + {0.5, 0.0, 0.0} + }; + + NeighborSearch ns; + ns.search_radius = 1.0; // use search radius = 1 for Check_Expand_Condition + ns.Check_Expand_Condition(ucell); + + // For identity lattice with search_radius=1 expected ceil produce values + EXPECT_EQ(ns.glayerX, 2); + EXPECT_EQ(ns.glayerY, 2); + EXPECT_EQ(ns.glayerZ, 2); + EXPECT_EQ(ns.glayerX_minus, 1); + + // Now populate all_atoms and check count + ns.setMemberVariables(ucell); + int images_x = ns.glayerX + ns.glayerX_minus; // iterations in x + int images_y = ns.glayerY + ns.glayerY_minus; + int images_z = ns.glayerZ + ns.glayerZ_minus; + int expected = images_x * images_y * images_z * 2; // 2 atoms per cell + EXPECT_EQ(static_cast(ns.all_atoms.size()), expected); +} + +TEST(NeighborSearchUnit, DistanceBox) +{ + NeighborSearch ns; + // set a single cell region at x=0..1,y=0..1,z=0..1 + ns.x = 0; ns.y = 0; ns.z = 0; + ns.wide_x = 1.0; ns.wide_y = 1.0; ns.wide_z = 1.0; + + double inside = ns.distance(0.2, 0.5, 0.5, 0.0, 0.0, 0.0); + EXPECT_DOUBLE_EQ(inside, 0.0); + + double outside = ns.distance(2.0, 0.5, 0.5, 0.0, 0.0, 0.0); + // squared distance should be (2-1)^2 = 1 + EXPECT_DOUBLE_EQ(outside, 1.0); +} + +TEST(NeighborSearchUnit, DecomposeCases) +{ + NeighborSearch ns; + int nx, ny, nz; + + ns.decompose(8, nx, ny, nz); + EXPECT_EQ(nx * ny * nz, 8); + // expect somewhat balanced cube factors for 8 + EXPECT_EQ(nx, 2); + EXPECT_EQ(ny, 2); + EXPECT_EQ(nz, 2); + + ns.decompose(7, nx, ny, nz); + EXPECT_EQ(nx * ny * nz, 7); + EXPECT_EQ(nx, 1); + EXPECT_EQ(ny, 1); + EXPECT_EQ(nz, 7); +} + +TEST(BinManagerUnit, InitAndBinning) +{ + std::vector inside; + std::vector ghost; + + // two atoms close to each other + inside.emplace_back(0.0, 0.0, 0.0, 0, 0, 0); + inside.emplace_back(0.5, 0.0, 0.0, 0, 1, 1); + + BinManager bm; + bm.init_bins(1.0, inside, ghost); + + EXPECT_GE(bm.nbinx, 1); + EXPECT_GE(bm.nbiny, 1); + EXPECT_GE(bm.nbinz, 1); + EXPECT_EQ(static_cast(bm.bins.size()), bm.nbinx * bm.nbiny * bm.nbinz); + + bm.do_binning(inside, ghost); + + // compute bin index for first atom + int ix = std::min(std::max(int((inside[0].position_x - bm.x_min) / bm.bin_sizex), 0), bm.nbinx - 1); + int iy = std::min(std::max(int((inside[0].position_y - bm.y_min) / bm.bin_sizey), 0), bm.nbiny - 1); + int iz = std::min(std::max(int((inside[0].position_z - bm.z_min) / bm.bin_sizez), 0), bm.nbinz - 1); + int idx = ix * bm.nbiny * bm.nbinz + iy * bm.nbinz + iz; + + EXPECT_GE(bm.bins[idx].atoms.size(), 1u); +} + +TEST(BinManagerUnit, BuildNeighborsAndClear) +{ + std::vector atoms; + atoms.emplace_back(0.0, 0.0, 0.0, 0, 0, 0); + atoms.emplace_back(0.5, 0.0, 0.0, 0, 1, 1); + atoms.emplace_back(5.0, 0.0, 0.0, 0, 2, 2); + + std::vector inside = atoms; + std::vector ghost; + + BinManager bm; + bm.init_bins(1.0, inside, ghost); + bm.do_binning(inside, ghost); + + NeighborList nl; + nl.initialize(static_cast(atoms.size()), 1024); + + bm.build_atom_neighbors(nl, atoms); + + // atom 0 and 1 are neighbors; atom 2 is far + EXPECT_GE(nl.numneigh[0], 1); + EXPECT_GE(nl.numneigh[1], 1); + EXPECT_EQ(nl.numneigh[2], 0); + + bm.clear(); + EXPECT_EQ(bm.bins.size(), 0u); +} \ No newline at end of file From f8d50a91856fc998b782246788061917625e64c4 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Wed, 29 Apr 2026 13:25:35 +0800 Subject: [PATCH 7/8] module_neighbor_search --- .../module_neighbor_search/neighbor_list.h | 10 +- .../test/test_neighbor_search.cpp | 205 ++++++++++++++++++ 2 files changed, 212 insertions(+), 3 deletions(-) diff --git a/source/source_cell/module_neighbor_search/neighbor_list.h b/source/source_cell/module_neighbor_search/neighbor_list.h index 0db8229ff3e..9ed0a4dc409 100644 --- a/source/source_cell/module_neighbor_search/neighbor_list.h +++ b/source/source_cell/module_neighbor_search/neighbor_list.h @@ -48,6 +48,11 @@ class PageAllocator { int* allocate(int n) { if (n <= 0) return nullptr; + // reject requests larger than a single page + if (n > pgsize) { + std::cerr << "PageAllocator::allocate error: request " << n << " larger than page size " << pgsize << std::endl; + return nullptr; + } if (pages.empty()) new_page(); Page& p = pages.back(); if (p.offset + n > p.capacity) { @@ -60,9 +65,8 @@ class PageAllocator { } void reset() { - for (auto& p : pages) { - p.offset = 0; - } + pages.resize(1); + pages[0].offset = 0; } }; diff --git a/source/source_cell/module_neighbor_search/test/test_neighbor_search.cpp b/source/source_cell/module_neighbor_search/test/test_neighbor_search.cpp index 78eb0d5772c..56bc7f5b93b 100644 --- a/source/source_cell/module_neighbor_search/test/test_neighbor_search.cpp +++ b/source/source_cell/module_neighbor_search/test/test_neighbor_search.cpp @@ -193,4 +193,209 @@ TEST(BinManagerUnit, BuildNeighborsAndClear) bm.clear(); EXPECT_EQ(bm.bins.size(), 0u); +} + +// Additional, more exhaustive tests + +TEST(PageAllocatorTest, MultiPageAllocationAndReset) +{ + PageAllocator pa(4); + int* a = pa.allocate(3); + ASSERT_NE(a, nullptr); + int* b = pa.allocate(3); + ASSERT_NE(b, nullptr); + // since page size = 4, a uses 3, b must be allocated on a new page + EXPECT_NE(b, a + 3); + + pa.reset(); + int* c = pa.allocate(2); + ASSERT_NE(c, nullptr); + // after reset, allocation starts at the beginning of the last page (current implementation) + EXPECT_EQ(c, pa.pages[0].data.data()); +} + +TEST(PageAllocatorTest, ZeroAndExactPageAllocations) +{ + PageAllocator pa(8); + // zero or negative allocation should return nullptr + EXPECT_EQ(pa.allocate(0), nullptr); + EXPECT_EQ(pa.allocate(-1), nullptr); + + // allocate exactly page size + int* p = pa.allocate(8); + ASSERT_NE(p, nullptr); + // next tiny allocation must go to a new page + int* q = pa.allocate(1); + ASSERT_NE(q, nullptr); + EXPECT_NE(q, p + 8); +} + +TEST(PageAllocatorTest, LargeAllocationSpansPages) +{ + PageAllocator pa(16); + // request larger than single page -> should still succeed (allocates new page and fits) + int* p = pa.allocate(32); + EXPECT_EQ(p, nullptr); +} + +TEST(NeighborListUnit, InitializeZeroAndResize) +{ + NeighborList nl; + nl.initialize(0, 8); + EXPECT_EQ(nl.numneigh.size(), 0); + EXPECT_EQ(nl.firstneigh.size(), 0); + + // initialize with larger size + nl.initialize(5, 8); + EXPECT_EQ(nl.numneigh.size(), 5); + EXPECT_EQ(nl.firstneigh.size(), 5); +} + +TEST(BinManagerUnit, EmptyAtomsBuildNeighbors) +{ + std::vector atoms; + std::vector ghost; + + BinManager bm; + bm.init_bins(1.0, atoms, ghost); + + NeighborList nl; + nl.initialize(0, 16); + + // should not crash or assert when atoms is empty + bm.build_atom_neighbors(nl, atoms); + EXPECT_EQ(nl.numneigh.size(), 0); +} + +TEST(BinManagerUnit, BoundaryAndExactRadius) +{ + // inside atom at origin; other atoms placed on bin boundaries and at exact radius + std::vector atoms; + atoms.emplace_back(0.0, 0.0, 0.0, 0, 0, 0); + // exactly at search radius 1.0 along x direction + atoms.emplace_back(1.0, 0.0, 0.0, 0, 1, 1); + // slightly inside + atoms.emplace_back(0.9, 0.0, 0.0, 0, 2, 2); + + std::vector inside = atoms; + std::vector ghost; + + BinManager bm; + bm.init_bins(1.0, inside, ghost); + bm.do_binning(inside, ghost); + + NeighborList nl; + nl.initialize(static_cast(inside.size()), 64); + + bm.build_atom_neighbors(nl, inside); + + // atom at distance exactly 1.0 should be considered neighbor (dist2 == sradius2) + EXPECT_GE(nl.numneigh[0], 1); + // the exact atom itself must not be counted as its own neighbor + for (int i = 0; i < static_cast(inside.size()); ++i) { + for (int j = 0; j < nl.numneigh[i]; ++j) { + int id = nl.firstneigh[i][j]; + EXPECT_NE(id, inside[i].atom_id); + } + } +} + +TEST(NeighborListUnit, AllocateManyNeighbors) +{ + // verify neighbor_list allocator can handle many small allocations across pages + NeighborList nl; + nl.initialize(10, 8); // small page size to force multiple pages + + // simulate allocations of varying sizes + for (int i = 0; i < 10; ++i) { + int count = (i % 4) + 1; + int* ptr = nl.allocator.allocate(count); + ASSERT_NE(ptr, nullptr); + // write and read back + for (int k = 0; k < count; ++k) ptr[k] = i * 10 + k; + } +} + +TEST(NeighborSearchUnit, NonOrthogonalLatticeExpand) +{ + UnitCellPlus ucell; + ucell.lat0 = 1.0; + ucell.omega = 1.0; + // skewed lattice + ucell.latvec.e11 = 1; ucell.latvec.e12 = 0.3; ucell.latvec.e13 = 0.0; + ucell.latvec.e21 = 0.1; ucell.latvec.e22 = 1.0; ucell.latvec.e23 = 0.0; + ucell.latvec.e31 = 0.0; ucell.latvec.e32 = 0.0; ucell.latvec.e33 = 1.0; + + ucell.ntype = 1; + ucell.na = {1}; + ucell.nat = 1; + ucell.tau = {{0.0, 0.0, 0.0}}; + + NeighborSearch ns; + ns.search_radius = 2.5; + ns.Check_Expand_Condition(ucell); + // for skewed lattice, expansion layers should be >= 1 + EXPECT_GE(ns.glayerX, 1); + EXPECT_GE(ns.glayerY, 1); + EXPECT_GE(ns.glayerZ, 1); +} + + + +TEST(NeighborSearchUnit, UCellToInputAtomsMultipleTypes) +{ + UnitCellPlus ucell; + ucell.lat0 = 1.0; + ucell.omega = 1.0; + ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0; + ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0; + ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1; + + ucell.ntype = 2; + ucell.na = {1, 2}; + ucell.nat = 3; + ucell.tau = { + {0.0, 0.0, 0.0}, + {0.5, 0.0, 0.0}, + {0.0, 0.5, 0.0} + }; + + NeighborSearch ns; + auto inputs = ns.ucell_to_input_atoms(ucell); + + EXPECT_EQ(inputs.n_atoms, 3); + ASSERT_EQ(inputs.InputAtom.size(), 3); + EXPECT_DOUBLE_EQ(inputs.InputAtom[2].position_y, 0.5); +} + +TEST(NeighborSearchUnit, DecomposePrimeNumber) +{ + NeighborSearch ns; + int nx, ny, nz; + ns.decompose(13, nx, ny, nz); + EXPECT_EQ(nx * ny * nz, 13); + EXPECT_EQ(nx, 1); + EXPECT_EQ(ny, 1); + EXPECT_EQ(nz, 13); +} + +TEST(NeighborSearchIntegration, GhostAtomsCountedAsNeighbors) +{ + // inside atom at origin; ghost atom within search radius should be found + std::vector inside; + std::vector ghost; + + inside.emplace_back(0.0, 0.0, 0.0, 0, 0, 0); + ghost.emplace_back(0.4, 0.0, 0.0, 0, 1, 1); + + BinManager bm; + bm.init_bins(1.0, inside, ghost); + bm.do_binning(inside, ghost); + + NeighborList nl; + nl.initialize(static_cast(inside.size()), 32); + + bm.build_atom_neighbors(nl, inside); + + EXPECT_GE(nl.numneigh[0], 1); } \ No newline at end of file From 6fb0adafba57c26fba66db292b3343e05600d882 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Tue, 5 May 2026 21:08:35 +0800 Subject: [PATCH 8/8] module_neighbor_search --- .../test/CMakeLists.txt | 23 +- .../test/bin_manager_test.cpp | 229 +++++++++++++ .../test/neighbor_list_test.cpp | 91 +++++ .../test/neighbor_search_test.cpp | 319 ++++++++++++++++++ 4 files changed, 659 insertions(+), 3 deletions(-) create mode 100644 source/source_cell/module_neighbor_search/test/bin_manager_test.cpp create mode 100644 source/source_cell/module_neighbor_search/test/neighbor_list_test.cpp create mode 100644 source/source_cell/module_neighbor_search/test/neighbor_search_test.cpp diff --git a/source/source_cell/module_neighbor_search/test/CMakeLists.txt b/source/source_cell/module_neighbor_search/test/CMakeLists.txt index ef45038fc8e..a8462c59f7f 100644 --- a/source/source_cell/module_neighbor_search/test/CMakeLists.txt +++ b/source/source_cell/module_neighbor_search/test/CMakeLists.txt @@ -4,13 +4,30 @@ remove_definitions(-D__EXX) +# Split tests: create one test target per test source to keep failures isolated. + AddTest( TARGET MODULE_CELL_NEIGHBOR_neighbor_search - LIBS parameter ${math_libs} base device + LIBS parameter ${math_libs} base device + SOURCES + neighbor_search_test.cpp + ../neighbor_search.cpp + ../bin_manager.cpp +) + +AddTest( + TARGET MODULE_CELL_NEIGHBOR_bin_manager + LIBS parameter ${math_libs} base device SOURCES - test_neighbor_search.cpp - ../neighbor_search.cpp + bin_manager_test.cpp ../bin_manager.cpp ) +AddTest( + TARGET MODULE_CELL_NEIGHBOR_allocator_and_list + LIBS parameter ${math_libs} base device + SOURCES + neighbor_list_test.cpp +) + diff --git a/source/source_cell/module_neighbor_search/test/bin_manager_test.cpp b/source/source_cell/module_neighbor_search/test/bin_manager_test.cpp new file mode 100644 index 00000000000..51e5e5e7400 --- /dev/null +++ b/source/source_cell/module_neighbor_search/test/bin_manager_test.cpp @@ -0,0 +1,229 @@ +#include +#include "../bin_manager.h" +#include "../neighbor_list.h" + +TEST(BinManagerUnit, InitAndBinning) +{ + std::vector inside; + std::vector ghost; + + // two atoms close to each other + inside.emplace_back(0.0, 0.0, 0.0, 0, 0, 0); + inside.emplace_back(0.5, 0.0, 0.0, 0, 1, 1); + + BinManager bm; + bm.init_bins(1.0, inside, ghost); + + EXPECT_EQ(bm.nbinx, 1); + EXPECT_EQ(bm.nbiny, 1); + EXPECT_EQ(bm.nbinz, 1); + EXPECT_EQ(static_cast(bm.bins.size()), bm.nbinx * bm.nbiny * bm.nbinz); + + bm.do_binning(inside, ghost); + + // compute bin index for first atom + int ix = std::min(std::max(int((inside[0].position_x - bm.x_min) / bm.bin_sizex), 0), bm.nbinx - 1); + int iy = std::min(std::max(int((inside[0].position_y - bm.y_min) / bm.bin_sizey), 0), bm.nbiny - 1); + int iz = std::min(std::max(int((inside[0].position_z - bm.z_min) / bm.bin_sizez), 0), bm.nbinz - 1); + int idx = ix * bm.nbiny * bm.nbinz + iy * bm.nbinz + iz; + + EXPECT_GE(bm.bins[idx].atoms.size(), 1u); +} + +TEST(BinManagerUnit, InitBins) +{ + std::vector atoms; + atoms.emplace_back(0.0, 0.0, 0.0, 0, 0, 0); + atoms.emplace_back(0.5, 0.0, 0.0, 0, 1, 1); + atoms.emplace_back(4.9, 0.0, 0.0, 0, 2, 2); + + std::vector inside = atoms; + std::vector ghost; + + BinManager bm; + bm.init_bins(1.0, inside, ghost); + EXPECT_EQ(bm.nbinx, 5); + EXPECT_EQ(bm.nbiny, 1); + EXPECT_EQ(bm.nbinz, 1); +} + +TEST(BinManagerUnit, BuildNeighborsAndClear) +{ + std::vector atoms; + atoms.emplace_back(0.0, 0.0, 0.0, 0, 0, 0); + atoms.emplace_back(0.5, 0.0, 0.0, 0, 1, 1); + atoms.emplace_back(5.0, 0.0, 0.0, 0, 2, 2); + + std::vector inside = atoms; + std::vector ghost; + + BinManager bm; + bm.init_bins(1.0, inside, ghost); + EXPECT_EQ(bm.nbinx, 5); + EXPECT_EQ(bm.nbiny, 1); + EXPECT_EQ(bm.nbinz, 1); + EXPECT_EQ(static_cast(bm.bins.size()), bm.nbinx * bm.nbiny * bm.nbinz); + + bm.do_binning(inside, ghost); + + NeighborList nl; + nl.initialize(static_cast(atoms.size()), 1024); + + bm.build_atom_neighbors(nl, atoms); + + // atom 0 and 1 are neighbors; atom 2 is far + EXPECT_EQ(nl.numneigh[0], 1); + EXPECT_EQ(nl.numneigh[1], 1); + EXPECT_EQ(nl.numneigh[2], 0); + + bm.clear(); + EXPECT_EQ(bm.bins.size(), 0u); +} + +TEST(BinManagerUnit, EmptyAtomsBuildNeighbors) +{ + std::vector atoms; + std::vector ghost; + + BinManager bm; + bm.init_bins(1.0, atoms, ghost); + + NeighborList nl; + nl.initialize(0, 16); + + // should not crash or assert when atoms is empty + bm.build_atom_neighbors(nl, atoms); + EXPECT_EQ(nl.numneigh.size(), 0); +} + +TEST(BinManagerUnit, BoundaryAndExactRadius) +{ + // inside atom at origin; other atoms placed on bin boundaries and at exact radius + std::vector atoms; + atoms.emplace_back(0.0, 0.0, 0.0, 0, 0, 0); + // exactly at search radius 1.0 along x direction + atoms.emplace_back(1.0, 0.0, 0.0, 0, 1, 1); + // slightly inside + atoms.emplace_back(0.9, 0.0, 0.0, 0, 2, 2); + + std::vector inside = atoms; + std::vector ghost; + + BinManager bm; + bm.init_bins(1.0, inside, ghost); + bm.do_binning(inside, ghost); + + NeighborList nl; + nl.initialize(static_cast(inside.size()), 64); + + bm.build_atom_neighbors(nl, inside); + + // atom at distance exactly 1.0 should be considered neighbor (dist2 == sradius2) + EXPECT_EQ(nl.numneigh[0], 2); + // the exact atom itself must not be counted as its own neighbor + for (int i = 0; i < static_cast(inside.size()); ++i) { + for (int j = 0; j < nl.numneigh[i]; ++j) { + int id = nl.firstneigh[i][j]; + EXPECT_NE(id, inside[i].atom_id); + } + } +} + +TEST(BinManagerUnit, InitWithGhostOnly) +{ + // inside empty, ghosts present -> init_bins should still compute bounds from ghosts + std::vector inside; + std::vector ghost; + + ghost.emplace_back(-1.0, -1.0, -1.0, 0, 0, 0); + ghost.emplace_back(2.0, 0.0, 0.0, 0, 1, 1); + + BinManager bm; + bm.init_bins(1.0, inside, ghost); + + EXPECT_EQ(bm.nbinx, 3); + EXPECT_EQ(bm.nbiny, 1); + EXPECT_EQ(bm.nbinz, 1); +} + +TEST(BinManagerUnit, BuildNeighborsNoNeighborsFirstneighNull) +{ + // atoms all far apart => no neighbors; firstneigh entries should be nullptr + std::vector atoms; + atoms.emplace_back(0.0, 0.0, 0.0, 0, 0, 0); + atoms.emplace_back(100.0, 100.0, 100.0, 0, 1, 1); + + std::vector inside = atoms; + std::vector ghost; + + BinManager bm; + bm.init_bins(1.0, inside, ghost); + bm.do_binning(inside, ghost); + + NeighborList nl; + nl.initialize(static_cast(inside.size()), 8); + + bm.build_atom_neighbors(nl, inside); + + EXPECT_EQ(nl.numneigh[0], 0); + EXPECT_EQ(nl.numneigh[1], 0); + EXPECT_EQ(nl.firstneigh[0], nullptr); + EXPECT_EQ(nl.firstneigh[1], nullptr); +} + +TEST(BinManagerUnit, GhostAtomsAreCounted) +{ + // inside atom at origin; ghost atom within search radius should be found + std::vector inside; + std::vector ghost; + + inside.emplace_back(0.0, 0.0, 0.0, 0, 0, 0); + ghost.emplace_back(0.4, 0.0, 0.0, 0, 1, 3); + + BinManager bm; + bm.init_bins(1.0, inside, ghost); + bm.do_binning(inside, ghost); + + NeighborList nl; + nl.initialize(static_cast(inside.size()), 32); + + bm.build_atom_neighbors(nl, inside); + + EXPECT_EQ(nl.numneigh.size(), 1); + EXPECT_EQ(nl.numneigh[0], 1); + // ensure neighbor id matches ghost atom id + bool found = false; + if (nl.numneigh[0] > 0 && nl.firstneigh[0] != nullptr) { + for (int k = 0; k < nl.numneigh[0]; ++k) { + if (nl.firstneigh[0][k] == 3) found = true; + } + } + EXPECT_TRUE(found); +} + +TEST(BinManagerUnit, MultipleBinsNeighborSearch) +{ + // Create a 3x3x3 grid of atoms spaced by 1.0 so multiple bins exist + std::vector atoms; + int id = 0; + for (int x = 0; x < 3; ++x) + for (int y = 0; y < 3; ++y) + for (int z = 0; z < 3; ++z) + atoms.emplace_back(x * 1.0, y * 1.0, z * 1.0, 0, 0, id++); + + std::vector inside = atoms; + std::vector ghost; + + BinManager bm; + bm.init_bins(1.0, inside, ghost); + bm.do_binning(inside, ghost); + + NeighborList nl; + nl.initialize(static_cast(inside.size()), 16); + + bm.build_atom_neighbors(nl, inside); + + // center atom (1,1,1) should have multiple neighbors + int center_index = 13; // 1*9 + 1*3 + 1 + EXPECT_EQ(nl.numneigh[center_index], 6); +} diff --git a/source/source_cell/module_neighbor_search/test/neighbor_list_test.cpp b/source/source_cell/module_neighbor_search/test/neighbor_list_test.cpp new file mode 100644 index 00000000000..8a9b7160543 --- /dev/null +++ b/source/source_cell/module_neighbor_search/test/neighbor_list_test.cpp @@ -0,0 +1,91 @@ +#include +#include "../neighbor_list.h" + +// Full coverage tests for PageAllocator and NeighborList + +TEST(PageAllocator_ConstructorsAndNewPage, DefaultAndCustom) +{ + PageAllocator pa_def; // default page size + ASSERT_EQ(pa_def.pages.size(), 1u); + EXPECT_EQ(pa_def.pages[0].capacity, pa_def.pgsize); + + PageAllocator pa(4); + EXPECT_EQ(pa.pgsize, 4); + size_t before = pa.pages.size(); + pa.new_page(); + EXPECT_EQ(pa.pages.size(), before + 1); +} + +TEST(PageAllocator_AllocateEdgeCases, ZeroNegativeAndTooLarge) +{ + PageAllocator pa(8); + EXPECT_EQ(pa.allocate(0), nullptr); + EXPECT_EQ(pa.allocate(-5), nullptr); + // larger than single page -> rejected + EXPECT_EQ(pa.allocate(9), nullptr); +} + +TEST(PageAllocator_AllocationBehavior, ExactPageAndOverflow) +{ + PageAllocator pa(4); + int* p1 = pa.allocate(4); + ASSERT_NE(p1, nullptr); + int* p2 = pa.allocate(1); + ASSERT_NE(p2, nullptr); + EXPECT_NE(p2, p1 + 4); + + PageAllocator pa2(3); + int* a = pa2.allocate(2); + ASSERT_NE(a, nullptr); + int* b = pa2.allocate(2); + ASSERT_NE(b, nullptr); + EXPECT_NE(b, a + 2); +} + +TEST(PageAllocator_ResetAndPages, ClearAndReset) +{ + PageAllocator pa(4); + pa.allocate(3); + pa.allocate(3); + ASSERT_EQ(pa.pages.size(), 2u); + + pa.reset(); + EXPECT_EQ(pa.pages.size(), 1u); + EXPECT_EQ(pa.pages[0].offset, 0); + + pa.pages.clear(); + int* p = pa.allocate(1); + ASSERT_NE(p, nullptr); + EXPECT_EQ(pa.pages.back().offset, 1); +} + +TEST(NeighborList_InitializeAndReset, Behavior) +{ + NeighborList nl; + nl.initialize(0, 16); + EXPECT_EQ(nl.nlocal, 0); + EXPECT_EQ(nl.numneigh.size(), 0u); + EXPECT_EQ(nl.firstneigh.size(), 0u); + + nl.initialize(5, 8); + EXPECT_EQ(nl.nlocal, 5); + EXPECT_EQ(nl.numneigh.size(), 5u); + EXPECT_EQ(nl.firstneigh.size(), 5u); + for (int i = 0; i < 5; ++i) { + EXPECT_EQ(nl.numneigh[i], 0); + EXPECT_EQ(nl.firstneigh[i], nullptr); + } + EXPECT_EQ(nl.allocator.pgsize, 8); + + int* storage = nl.allocator.allocate(3); + ASSERT_NE(storage, nullptr); + nl.numneigh[2] = 3; + nl.firstneigh[2] = storage; + + nl.reset(); + EXPECT_EQ(nl.allocator.pages.size(), 1u); + EXPECT_EQ(nl.allocator.pages[0].offset, 0); + EXPECT_EQ(nl.numneigh.size(), 5u); +} + + diff --git a/source/source_cell/module_neighbor_search/test/neighbor_search_test.cpp b/source/source_cell/module_neighbor_search/test/neighbor_search_test.cpp new file mode 100644 index 00000000000..5358f5e87b7 --- /dev/null +++ b/source/source_cell/module_neighbor_search/test/neighbor_search_test.cpp @@ -0,0 +1,319 @@ +#include +#include "../neighbor_search.h" +#include "../unitcell_plus.h" + +TEST(NeighborSearchTest, TwoAtomsNeighbor) +{ + UnitCellPlus ucell; + + // ===== 基本晶胞信息 ===== + ucell.lat0 = 1.0; + ucell.omega = 1.0; + + ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0; + ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0; + ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1; + + // ===== 原子信息 ===== + ucell.ntype = 1; + ucell.na = {2}; + ucell.nat = 2; + + ucell.tau = { + {0.0, 0.0, 0.0}, + {0.5, 0.0, 0.0} + }; + + NeighborSearch ns; + + double cutoff = 1.0; + + ns.init(ucell, cutoff, 0); + ns.build_neighbors(); + + auto &list = ns.get_neighbor_list(); + + ASSERT_EQ(list.numneigh.size(), 2); + + EXPECT_EQ(list.numneigh[0], 8); + EXPECT_EQ(list.numneigh[1], 8); +} + +TEST(NeighborSearchTest, NoNeighbor) +{ + UnitCellPlus ucell; + + ucell.lat0 = 1.0; + ucell.omega = 1.0; + + ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0; + ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0; + ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1; + + ucell.ntype = 1; + ucell.na = {2}; + ucell.nat = 2; + + ucell.tau = { + {0.0, 0.0, 0.0}, + {5.0, 0.0, 0.0} + }; + + NeighborSearch ns; + + // use a smaller search radius to avoid counting periodic-image neighbors + ns.init(ucell, 0.1, 0); + ns.build_neighbors(); + + auto &list = ns.get_neighbor_list(); + + EXPECT_EQ(list.numneigh[0], 0); + EXPECT_EQ(list.numneigh[1], 0); +} + +TEST(NeighborSearchUnit, UCellToInputAtoms) +{ + UnitCellPlus ucell; + ucell.lat0 = 1.0; + ucell.omega = 1.0; + ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0; + ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0; + ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1; + + ucell.ntype = 1; + ucell.na = {2}; + ucell.nat = 2; + ucell.tau = { + {0.0, 0.0, 0.0}, + {0.5, 0.0, 0.0} + }; + + NeighborSearch ns; + auto inputs = ns.ucell_to_input_atoms(ucell); + + EXPECT_EQ(inputs.n_atoms, 2); + ASSERT_EQ(inputs.InputAtom.size(), 2); + EXPECT_DOUBLE_EQ(inputs.InputAtom[0].position_x, 0.0); + EXPECT_DOUBLE_EQ(inputs.InputAtom[1].position_x, 0.5); + EXPECT_DOUBLE_EQ(inputs.x_low, 0.0); + EXPECT_DOUBLE_EQ(inputs.x_high, 0.5); +} + +TEST(NeighborSearchUnit, CheckExpandAndSetMembers) +{ + UnitCellPlus ucell; + ucell.lat0 = 1.0; + ucell.omega = 1.0; + ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0; + ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0; + ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1; + + ucell.ntype = 1; + ucell.na = {2}; + ucell.nat = 2; + ucell.tau = { + {0.0, 0.0, 0.0}, + {0.5, 0.0, 0.0} + }; + + NeighborSearch ns; + ns.search_radius = 1.0; // use search radius = 1 for Check_Expand_Condition + ns.Check_Expand_Condition(ucell); + + // For identity lattice with search_radius=1 expected ceil produce values + EXPECT_EQ(ns.glayerX, 2); + EXPECT_EQ(ns.glayerY, 2); + EXPECT_EQ(ns.glayerZ, 2); + EXPECT_EQ(ns.glayerX_minus, 1); + + // Now populate all_atoms and check count + ns.setMemberVariables(ucell); + int images_x = ns.glayerX + ns.glayerX_minus; // iterations in x + int images_y = ns.glayerY + ns.glayerY_minus; + int images_z = ns.glayerZ + ns.glayerZ_minus; + int expected = images_x * images_y * images_z * 2; // 2 atoms per cell + EXPECT_EQ(static_cast(ns.all_atoms.size()), expected); +} + +TEST(NeighborSearchUnit, DistanceBox) +{ + NeighborSearch ns; + // set a single cell region at x=0..1,y=0..1,z=0..1 + ns.x = 0; ns.y = 0; ns.z = 0; + ns.wide_x = 1.0; ns.wide_y = 1.0; ns.wide_z = 1.0; + + double inside = ns.distance(0.2, 0.5, 0.5, 0.0, 0.0, 0.0); + EXPECT_DOUBLE_EQ(inside, 0.0); + + double outside = ns.distance(2.0, 0.5, 0.5, 0.0, 0.0, 0.0); + // squared distance should be (2-1)^2 = 1 + EXPECT_DOUBLE_EQ(outside, 1.0); +} + +TEST(NeighborSearchUnit, DecomposeCases) +{ + NeighborSearch ns; + int nx, ny, nz; + + ns.decompose(8, nx, ny, nz); + EXPECT_EQ(nx * ny * nz, 8); + // expect somewhat balanced cube factors for 8 + EXPECT_EQ(nx, 2); + EXPECT_EQ(ny, 2); + EXPECT_EQ(nz, 2); + + ns.decompose(7, nx, ny, nz); + EXPECT_EQ(nx * ny * nz, 7); + EXPECT_EQ(nx, 1); + EXPECT_EQ(ny, 1); + EXPECT_EQ(nz, 7); +} + +TEST(NeighborSearchUnit, UCellToInputAtomsMultipleTypes) +{ + UnitCellPlus ucell; + ucell.lat0 = 1.0; + ucell.omega = 1.0; + ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0; + ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0; + ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1; + + ucell.ntype = 2; + ucell.na = {1, 2}; + ucell.nat = 3; + ucell.tau = { + {0.0, 0.0, 0.0}, + {0.5, 0.0, 0.0}, + {0.0, 0.5, 0.0} + }; + + NeighborSearch ns; + auto inputs = ns.ucell_to_input_atoms(ucell); + + EXPECT_EQ(inputs.n_atoms, 3); + ASSERT_EQ(inputs.InputAtom.size(), 3); + EXPECT_DOUBLE_EQ(inputs.InputAtom[2].position_y, 0.5); +} + +TEST(NeighborSearchUnit, DecomposePrimeNumber) +{ + NeighborSearch ns; + int nx, ny, nz; + ns.decompose(13, nx, ny, nz); + EXPECT_EQ(nx * ny * nz, 13); + EXPECT_EQ(nx, 1); + EXPECT_EQ(ny, 1); + EXPECT_EQ(nz, 13); +} + +TEST(NeighborSearchUnit, NonOrthogonalLatticeExpand) +{ + UnitCellPlus ucell; + ucell.lat0 = 1.0; + ucell.omega = 1.0; + // skewed lattice + ucell.latvec.e11 = 1; ucell.latvec.e12 = 0.3; ucell.latvec.e13 = 0.0; + ucell.latvec.e21 = 0.1; ucell.latvec.e22 = 1.0; ucell.latvec.e23 = 0.0; + ucell.latvec.e31 = 0.0; ucell.latvec.e32 = 0.0; ucell.latvec.e33 = 1.0; + + ucell.ntype = 1; + ucell.na = {1}; + ucell.nat = 1; + ucell.tau = {{0.0, 0.0, 0.0}}; + + NeighborSearch ns; + ns.search_radius = 2.5; + ns.Check_Expand_Condition(ucell); + // for skewed lattice, expansion layers should be >= 1 + EXPECT_GE(ns.glayerX, 1); + EXPECT_GE(ns.glayerY, 1); + EXPECT_GE(ns.glayerZ, 1); +} + +// --- additional tests to cover remaining branches in neighbor_search.cpp --- + +TEST(NeighborSearchInit_WideZero_CentralInside, SingleAtomCell) +{ + UnitCellPlus ucell; + ucell.lat0 = 1.0; + ucell.omega = 1.0; + ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0; + ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0; + ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1; + + ucell.ntype = 1; + ucell.na = {1}; + ucell.nat = 1; + ucell.tau = {{0.0, 0.0, 0.0}}; + + NeighborSearch ns; + // choose sr small enough; with mpi_size fixed to 1 in init, wide_* become 0 + ns.init(ucell, 0.1, 0); + // central cell atom should be counted as inside + EXPECT_EQ(ns.inside_atoms.size(), 1); + EXPECT_EQ(ns.neighbor_list.nlocal, static_cast(ns.inside_atoms.size())); +} + +TEST(NeighborSearchInit_MpiRankIndexing, RankValues) +{ + UnitCellPlus ucell; + ucell.lat0 = 1.0; + ucell.omega = 1.0; + ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0; + ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0; + ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1; + + ucell.ntype = 1; + ucell.na = {1}; + ucell.nat = 1; + ucell.tau = {{0.0, 0.0, 0.0}}; + + NeighborSearch ns0; + ns0.init(ucell, 0.5, 0); + // with mpi_size fixed to 1 in init, nx=ny=nz=1; for rank 0 expect x=y=0,z=0 + EXPECT_EQ(ns0.x, 0); + EXPECT_EQ(ns0.y, 0); + EXPECT_EQ(ns0.z, 0); + + NeighborSearch ns1; + ns1.init(ucell, 0.5, 1); + // rank 1: z = mpi_rank /(nx*ny) = 1, x and y remain 0 when ny==1 + EXPECT_EQ(ns1.z, 1); + EXPECT_EQ(ns1.y, 0); + EXPECT_EQ(ns1.x, 0); +} + +TEST(NeighborSearchDistance_OutsideCases, VariousAxes) +{ + NeighborSearch ns; + ns.x = 0; ns.y = 0; ns.z = 0; + ns.wide_x = 2.0; ns.wide_y = 3.0; ns.wide_z = 4.0; + + // position inside box along x (no dx), but outside along y by above high bound + double d = ns.distance(0.5, 4.5, 1.0, 0.0, 0.0, 0.0); + // dy = position_y - (y_low + (y+1)*wide_y) = 4.5 - 3.0 = 1.5 -> squared 2.25 + // dx = 0, dz = 0 -> total 2.25 + EXPECT_DOUBLE_EQ(d, 2.25); + + // position left of low bound on x + double d2 = ns.distance(-1.0, 1.0, 1.0, 0.0, 0.0, 0.0); + // dx = x_low - position_x = 0 - (-1) = 1 -> squared 1 + EXPECT_DOUBLE_EQ(d2, 1.0); +} + +TEST(NeighborSearchDecompose_SmallSizes, TwoAndOne) +{ + NeighborSearch ns; + int nx, ny, nz; + ns.decompose(2, nx, ny, nz); + EXPECT_EQ(nx * ny * nz, 2); + // possible decomposition is nx=1, ny=1, nz=2 (or nx=1, ny=2, nz=1 depending on algorithm) + EXPECT_EQ(nx, 1); + + ns.decompose(1, nx, ny, nz); + EXPECT_EQ(nx, 1); + EXPECT_EQ(ny, 1); + EXPECT_EQ(nz, 1); +} + +// end of additional tests