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..6f54f71513f --- /dev/null +++ b/source/source_cell/module_neighbor_search/CMakeLists.txt @@ -0,0 +1,16 @@ +add_library( + neighbor_search + OBJECT + bin_manager.cpp + neighbor_search.cpp +) + +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/bin_manager.cpp b/source/source_cell/module_neighbor_search/bin_manager.cpp new file mode 100644 index 00000000000..b916c90a8ec --- /dev/null +++ b/source/source_cell/module_neighbor_search/bin_manager.cpp @@ -0,0 +1,191 @@ +#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); + + int 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(); + + //std::cout< +#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; + + 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..ed8d99c161a --- /dev/null +++ b/source/source_cell/module_neighbor_search/neighbor_atom.h @@ -0,0 +1,34 @@ +#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; + 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), + 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..9ed0a4dc409 --- /dev/null +++ b/source/source_cell/module_neighbor_search/neighbor_list.h @@ -0,0 +1,100 @@ +#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; + // 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) { + new_page(); + return allocate(n); + } + int* ptr = p.data.data() + p.offset; + p.offset += n; + return ptr; + } + + void reset() { + pages.resize(1); + pages[0].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); + // ensure neighbor containers are sized and initialized + numneigh.assign(n, 0); + firstneigh.assign(n, nullptr); + } + + 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..9775df965d8 --- /dev/null +++ b/source/source_cell/module_neighbor_search/neighbor_search.cpp @@ -0,0 +1,266 @@ +#include "source_cell/module_neighbor_search/neighbor_search.h" +#include +#include +#include + + +InputAtoms NeighborSearch::ucell_to_input_atoms(const IAtomProvider& ucell) +{ + InputAtoms input_atoms; + int atom_count = 0; + + 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.get_ntype(); i++) + { + for (int j = 0; j < ucell.get_na(i); j++) + { + NeighborAtom atom( + ucell.get_tauu(i,j).x, + ucell.get_tauu(i,j).y, + ucell.get_tauu(i,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 IAtomProvider& ucell, double sr, int mpi_rank) +{ + search_radius = sr / ucell.get_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; + + int in_x, in_y, in_z; + + for (int i = 0; i < all_atoms.size(); i++) + { + 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< 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; + + 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.get_ntype(); i++) + { + for (int j = 0; j < ucell.get_na(i); j++) + { + 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); + if(ix==0&&iy==0&&iz==0) + { + atom.is_inside = true; + } + else + { + atom.is_inside = false; + } + 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..29e0f771f6d --- /dev/null +++ b/source/source_cell/module_neighbor_search/neighbor_search.h @@ -0,0 +1,58 @@ +#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/module_neighbor_search/unitcell_interface.h" + +class NeighborSearch +{ +public: + NeighborSearch()=default; + ~NeighborSearch()=default; + + void init(const IAtomProvider& ucell, double sr, int mpi_rank); + + + void build_neighbors(); + + InputAtoms ucell_to_input_atoms(const IAtomProvider& ucell); + + void Check_Expand_Condition(const IAtomProvider& ucell); + + void setMemberVariables(const IAtomProvider& ucell); + NeighborList& get_neighbor_list() { return neighbor_list; } + + 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_cell/module_neighbor_search/test/CMakeLists.txt b/source/source_cell/module_neighbor_search/test/CMakeLists.txt new file mode 100644 index 00000000000..a8462c59f7f --- /dev/null +++ b/source/source_cell/module_neighbor_search/test/CMakeLists.txt @@ -0,0 +1,33 @@ +remove_definitions(-D__CUDA) +remove_definitions(-D__ROCM) +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 + 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 + 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 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 new file mode 100644 index 00000000000..56bc7f5b93b --- /dev/null +++ b/source/source_cell/module_neighbor_search/test/test_neighbor_search.cpp @@ -0,0 +1,401 @@ +#include +#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); +} + +// 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 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/module_neighbor_search/unitcell_plus.h b/source/source_cell/module_neighbor_search/unitcell_plus.h new file mode 100644 index 00000000000..0655b1fd623 --- /dev/null +++ b/source/source_cell/module_neighbor_search/unitcell_plus.h @@ -0,0 +1,52 @@ +#include "source_cell/module_neighbor_search/unitcell_interface.h" +#include + +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_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; diff --git a/source/source_esolver/esolver_lj.cpp b/source/source_esolver/esolver_lj.cpp index f71698586ea..bc885f19b8d 100644 --- a/source/source_esolver/esolver_lj.cpp +++ b/source/source_esolver/esolver_lj.cpp @@ -4,11 +4,35 @@ #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 { + 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 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 +141,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" diff --git a/source/source_esolver/esolver_lj.h b/source/source_esolver/esolver_lj.h index ca23a8bed80..57b39dce74c 100644 --- a/source/source_esolver/esolver_lj.h +++ b/source/source_esolver/esolver_lj.h @@ -2,6 +2,7 @@ #define ESOLVER_LJ_H #include "esolver.h" +#include "source_cell/module_neighbor_search/unitcell_plus.h" namespace ModuleESolver { @@ -14,6 +15,8 @@ namespace ModuleESolver classname = "ESolver_LJ"; } + UnitCellPlus change_from_ucell_to_ucell_plus(UnitCell& ucell); + void before_all_runners(UnitCell& ucell, const Input_para& inp) override; void runner(UnitCell& cell, const int istep) override; 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