diff --git a/source/source_cell/module_neighbor/sltk_grid.cpp b/source/source_cell/module_neighbor/sltk_grid.cpp index 2fe7a77b8ea..25355ea1d8c 100644 --- a/source/source_cell/module_neighbor/sltk_grid.cpp +++ b/source/source_cell/module_neighbor/sltk_grid.cpp @@ -4,6 +4,126 @@ #include "source_base/global_variable.h" #include "source_base/timer.h" +#include +#include +#include +#include + +namespace +{ +constexpr double minimum_cartesian_bin_edge = 0.1; +constexpr std::size_t max_cartesian_bin_count = 1000000; + +int count_bins_for_range(const double coordinate_range, const double box_edge_length) +{ + if (coordinate_range <= 0.0 || box_edge_length <= 0.0) + { + return 1; + } + + const double bin_count = std::floor(coordinate_range / box_edge_length) + 1.0; + if (!std::isfinite(bin_count) + || bin_count > static_cast(std::numeric_limits::max())) + { + return std::numeric_limits::max(); + } + return std::max(1, static_cast(bin_count)); +} + +std::size_t saturated_bin_count_product(const int nx, const int ny, const int nz) +{ + if (nx <= 0 || ny <= 0 || nz <= 0) + { + return 0; + } + + const std::size_t sx = static_cast(nx); + const std::size_t sy = static_cast(ny); + const std::size_t sz = static_cast(nz); + const std::size_t max_size = std::numeric_limits::max(); + if (sx > max_size / sy) + { + return max_size; + } + const std::size_t xy = sx * sy; + if (xy > max_size / sz) + { + return max_size; + } + return xy * sz; +} + +void set_box_counts_for_edge(const double range_x, + const double range_y, + const double range_z, + const double box_edge_length, + int& box_nx, + int& box_ny, + int& box_nz) +{ + box_nx = count_bins_for_range(range_x, box_edge_length); + box_ny = count_bins_for_range(range_y, box_edge_length); + box_nz = count_bins_for_range(range_z, box_edge_length); +} + +double estimate_bin_count_for_edge(const double range_x, + const double range_y, + const double range_z, + const double box_edge_length) +{ + const double nx = std::floor(range_x / box_edge_length) + 1.0; + const double ny = std::floor(range_y / box_edge_length) + 1.0; + const double nz = std::floor(range_z / box_edge_length) + 1.0; + return nx * ny * nz; +} + +double coarsen_edge_for_bin_cap(const double initial_edge, + const double range_x, + const double range_y, + const double range_z) +{ + double edge = initial_edge; + const double max_range = std::max(range_x, std::max(range_y, range_z)); + + for (int iteration = 0; iteration < 16; ++iteration) + { + const double estimated_count = estimate_bin_count_for_edge(range_x, range_y, range_z, edge); + if (std::isfinite(estimated_count) + && estimated_count <= static_cast(max_cartesian_bin_count)) + { + return edge; + } + + if (!std::isfinite(estimated_count)) + { + return std::max(edge, max_range); + } + + const double scale = std::max(2.0, + std::cbrt(estimated_count + / static_cast(max_cartesian_bin_count))); + edge *= scale; + } + + return std::max(edge, max_range); +} + +double squared_distance(const FAtom& lhs, const FAtom& rhs) +{ + const double delta_x = lhs.x - rhs.x; + const double delta_y = lhs.y - rhs.y; + const double delta_z = lhs.z - rhs.z; + return delta_x * delta_x + delta_y * delta_y + delta_z * delta_z; +} + +bool is_neighbor_by_current_rule(const double distance2, const double radius2) +{ + const bool is_not_zero_distance = distance2 != 0.0; + const bool is_within_search_radius = distance2 <= radius2; + return is_not_zero_distance && is_within_search_radius; +} +} // namespace + Grid::Grid(const int& test_grid_in) : test_grid(test_grid_in) { } @@ -37,6 +157,13 @@ void Grid::Check_Expand_Condition(const UnitCell& ucell) { // ModuleBase::TITLE(GlobalV::ofs_running, "Atom_input", "Check_Expand_Condition"); + glayerX = 0; + glayerX_minus = 0; + glayerY = 0; + glayerY_minus = 0; + glayerZ = 0; + glayerZ_minus = 0; + if (!pbc) { return; @@ -199,27 +326,32 @@ void Grid::setMemberVariables(std::ofstream& ofs_in, // output data to ofs ModuleBase::GlobalFunc::OUT(ofs_in, "Min coordinates of atoms", x_min, y_min, z_min); ModuleBase::GlobalFunc::OUT(ofs_in, "Max coordinates of atoms", x_max, y_max, z_max); - this->box_edge_length = sradius + 0.1; // To avoid edge cases, the size of the box is slightly increased. + // Keep each Cartesian bin at least as wide as the search radius so a 27-bin scan is sufficient. + // Coarsen very sparse AABB grids to avoid allocating millions of empty bins. + this->box_edge_length = std::max(sradius, minimum_cartesian_bin_edge); -/* warning box algorithm - this->box_nx = std::ceil((this->x_max - this->x_min) / box_edge_length) + 1; - this->box_ny = std::ceil((this->y_max - this->y_min) / box_edge_length) + 1; - this->box_nz = std::ceil((this->z_max - this->z_min) / box_edge_length) + 1; - ModuleBase::GlobalFunc::OUT(ofs_in, "BoxNumber", box_nx, box_ny, box_nz); - - atoms_in_box.resize(this->box_nx); - for (int i = 0; i < this->box_nx; i++) + const bool has_expanded_images = (glayerX + glayerX_minus) > 0 && (glayerY + glayerY_minus) > 0 + && (glayerZ + glayerZ_minus) > 0; + if (has_expanded_images) { - atoms_in_box[i].resize(this->box_ny); - for (int j = 0; j < this->box_ny; j++) + const double range_x = std::max(0.0, this->x_max - this->x_min); + const double range_y = std::max(0.0, this->y_max - this->y_min); + const double range_z = std::max(0.0, this->z_max - this->z_min); + + this->box_edge_length = coarsen_edge_for_bin_cap(this->box_edge_length, range_x, range_y, range_z); + set_box_counts_for_edge(range_x, range_y, range_z, this->box_edge_length, box_nx, box_ny, box_nz); + while (saturated_bin_count_product(box_nx, box_ny, box_nz) > max_cartesian_bin_count) { - atoms_in_box[i][j].resize(this->box_nz); + this->box_edge_length *= 2.0; + set_box_counts_for_edge(range_x, range_y, range_z, this->box_edge_length, box_nx, box_ny, box_nz); } } - */ - this->box_nx = glayerX + glayerX_minus; - this->box_ny = glayerY + glayerY_minus; - this->box_nz = glayerZ + glayerZ_minus; + else + { + this->box_nx = 0; + this->box_ny = 0; + this->box_nz = 0; + } ModuleBase::GlobalFunc::OUT(ofs_in, "Number of needed cells", box_nx, box_ny, box_nz); atoms_in_box.resize(this->box_nx); @@ -245,18 +377,17 @@ void Grid::setMemberVariables(std::ofstream& ofs_in, // output data to ofs 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; FAtom atom(x, y, z, i, j, ix, iy, iz); - int box_i_x, box_i_y, box_i_z; - //this->getBox(box_i_x, box_i_y, box_i_z, x, y, z); - box_i_x = ix + glayerX_minus; - box_i_y = iy + glayerY_minus; - box_i_z = iz + glayerZ_minus; + int box_i_x = 0; + int box_i_y = 0; + int box_i_z = 0; + this->getBox(box_i_x, box_i_y, box_i_z, x, y, z); this->atoms_in_box[box_i_x][box_i_y][box_i_z].push_back(atom); } } } } } - + this->all_adj_info.resize(ucell.ntype); for (int i = 0; i < ucell.ntype; i++) { @@ -267,66 +398,73 @@ void Grid::setMemberVariables(std::ofstream& ofs_in, // output data to ofs void Grid::Construct_Adjacent(const UnitCell& ucell) { ModuleBase::timer::start("Grid", "constru_adj"); - - for (int i_type = 0; i_type < ucell.ntype; i_type++) + for (int i_type = 0; i_type < ucell.ntype; i_type++) { for (int j_atom = 0; j_atom < ucell.atoms[i_type].na; j_atom++) { - FAtom atom(ucell.atoms[i_type].tau[j_atom].x, - ucell.atoms[i_type].tau[j_atom].y, - ucell.atoms[i_type].tau[j_atom].z, - i_type, - j_atom, - 0, 0 ,0); - - this->Construct_Adjacent_near_box(atom); + ucell.atoms[i_type].tau[j_atom].y, + ucell.atoms[i_type].tau[j_atom].z, + i_type, + j_atom, + 0, + 0, + 0); + + std::vector& adjacent_atoms = this->all_adj_info[i_type][j_atom]; + this->Construct_Adjacent_near_box(atom, adjacent_atoms); } } ModuleBase::timer::end("Grid", "constru_adj"); } -void Grid::Construct_Adjacent_near_box(const FAtom& fatom) +void Grid::Construct_Adjacent_near_box(const FAtom& fatom, std::vector& adjacent_atoms) { - ModuleBase::timer::start("Grid", "adj_near_box"); + if (box_nx <= 0 || box_ny <= 0 || box_nz <= 0) + { + return; + } + int box_i_x=0; int box_i_y=0; int box_i_z=0; this->getBox(box_i_x, box_i_y, box_i_z, fatom.x, fatom.y, fatom.z); - for (int box_i_x_adj = 0; box_i_x_adj < glayerX + glayerX_minus; box_i_x_adj++) + const int box_x_begin = std::max(0, box_i_x - 1); + const int box_x_end = std::min(box_nx - 1, box_i_x + 1); + const int box_y_begin = std::max(0, box_i_y - 1); + const int box_y_end = std::min(box_ny - 1, box_i_y + 1); + const int box_z_begin = std::max(0, box_i_z - 1); + const int box_z_end = std::min(box_nz - 1, box_i_z + 1); + + for (int box_i_x_adj = box_x_begin; box_i_x_adj <= box_x_end; box_i_x_adj++) { - for (int box_i_y_adj = 0; box_i_y_adj < glayerY + glayerY_minus; box_i_y_adj++) + for (int box_i_y_adj = box_y_begin; box_i_y_adj <= box_y_end; box_i_y_adj++) { - for (int box_i_z_adj = 0; box_i_z_adj < glayerZ + glayerZ_minus; box_i_z_adj++) + for (int box_i_z_adj = box_z_begin; box_i_z_adj <= box_z_end; box_i_z_adj++) { - for (auto &fatom2 : this->atoms_in_box[box_i_x_adj][box_i_y_adj][box_i_z_adj]) + for (auto& fatom2 : this->atoms_in_box[box_i_x_adj][box_i_y_adj][box_i_z_adj]) { - this->Construct_Adjacent_final(fatom, &fatom2); + this->Construct_Adjacent_final(fatom, &fatom2, adjacent_atoms); } } } } - ModuleBase::timer::end("Grid", "adj_near_box"); } void Grid::Construct_Adjacent_final(const FAtom& fatom1, - FAtom* fatom2) + FAtom* fatom2, + std::vector& adjacent_atoms) const { - double delta_x = fatom1.x - fatom2->x; - double delta_y = fatom1.y - fatom2->y; - double delta_z = fatom1.z - fatom2->z; - - double dr = delta_x * delta_x + delta_y * delta_y + delta_z * delta_z; - + const double distance2 = squared_distance(fatom1, *fatom2); // 20241204 zhanghaochong // dr == 0 means the same atom // the atom itself is neighbour atom, but the order itself must on last in the list. // so we will add itself on find atom function, and skip here. // I dont know why, but if we add self here, test 701_LJ_MD_Anderson will assert - if (dr != 0.0 && dr <= this->sradius2) + if (is_neighbor_by_current_rule(distance2, this->sradius2)) { - all_adj_info[fatom1.type][fatom1.natom].push_back(fatom2); + adjacent_atoms.push_back(fatom2); } } diff --git a/source/source_cell/module_neighbor/sltk_grid.h b/source/source_cell/module_neighbor/sltk_grid.h index 4b2da582800..69d2ab40753 100644 --- a/source/source_cell/module_neighbor/sltk_grid.h +++ b/source/source_cell/module_neighbor/sltk_grid.h @@ -4,10 +4,12 @@ #include "source_cell/unitcell.h" #include "sltk_atom.h" +#include #include #include #include #include +#include typedef std::vector AtomMap; @@ -22,6 +24,19 @@ class Grid Grid& operator=(Grid&&) = default; + /** + * @brief Build the expanded-cell neighbor table for one UnitCell. + * + * @param [out] ofs Output stream used by ABACUS running logs. + * @param [in] ucell Unit cell whose atom positions are in lat0-scaled + * Cartesian coordinates. + * @param [in] radius_in Search radius in lat0 units. + * @param [in] boundary Whether periodic boundary images are included. + * + * @details The current implementation first estimates periodic image + * layers, then builds atom images in a Cartesian cell list, and finally + * constructs all_adj_info for every center atom. + */ void init(std::ofstream& ofs, const UnitCell& ucell, const double radius_in, const bool boundary = true); // Data @@ -37,18 +52,53 @@ class Grid double y_max=0.0; double z_max=0.0; - // The algorithm for searching neighboring atoms uses a "box" partitioning method. - // Each box has an edge length of sradius, and the number of boxes in each direction is recorded here. + // Cartesian cell-list partitioning. Each bin edge is at least sradius so a + // 27-bin stencil is sufficient for Euclidean cutoff search. double box_edge_length=0.0; int box_nx=0; int box_ny=0; int box_nz=0; - void getBox(int& bx, int& by, int& bz, const double& x, const double& y, const double& z) + void getBox(int& bx, int& by, int& bz, const double& x, const double& y, const double& z) const { - bx = std::floor((x - x_min) / box_edge_length); - by = std::floor((y - y_min) / box_edge_length); - bz = std::floor((z - z_min) / box_edge_length); + if (box_edge_length <= 0.0 || box_nx <= 0 || box_ny <= 0 || box_nz <= 0) + { + bx = 0; + by = 0; + bz = 0; + return; + } + + bx = static_cast(std::floor((x - x_min) / box_edge_length)); + by = static_cast(std::floor((y - y_min) / box_edge_length)); + bz = static_cast(std::floor((z - z_min) / box_edge_length)); + + if (bx < 0) + { + bx = 0; + } + else if (bx >= box_nx) + { + bx = box_nx - 1; + } + + if (by < 0) + { + by = 0; + } + else if (by >= box_ny) + { + by = box_ny - 1; + } + + if (bz < 0) + { + bz = 0; + } + else if (bz >= box_nz) + { + bz = box_nz - 1; + } } // Stores the atoms after box partitioning. std::vector>> atoms_in_box; @@ -96,12 +146,32 @@ class Grid private: int test_grid; + /** + * @brief Build expanded atom images and resize neighbor storage. + */ void setMemberVariables(std::ofstream& ofs_in, const UnitCell& ucell); + /** + * @brief Construct neighbor pointer lists for all atoms in the original cell. + */ void Construct_Adjacent(const UnitCell& ucell); - void Construct_Adjacent_near_box(const FAtom& fatom); - void Construct_Adjacent_final(const FAtom& fatom1, FAtom* fatom2); + /** + * @brief Scan the Cartesian cell-list stencil for one center atom. + */ + void Construct_Adjacent_near_box(const FAtom& fatom, std::vector& adjacent_atoms); + + /** + * @brief Append fatom2 to adjacent_atoms when the current distance rule + * accepts the pair. + */ + void Construct_Adjacent_final(const FAtom& fatom1, + FAtom* fatom2, + std::vector& adjacent_atoms) const; + + /** + * @brief Estimate positive and negative periodic image layers needed for sradius. + */ void Check_Expand_Condition(const UnitCell& ucell); int glayerX=0; int glayerX_minus=0; diff --git a/source/source_cell/module_neighbor/test/CMakeLists.txt b/source/source_cell/module_neighbor/test/CMakeLists.txt index 3d891aa20c3..27a5ca5494b 100644 --- a/source/source_cell/module_neighbor/test/CMakeLists.txt +++ b/source/source_cell/module_neighbor/test/CMakeLists.txt @@ -22,4 +22,12 @@ AddTest( SOURCES sltk_atom_arrange_test.cpp ../sltk_atom_arrange.cpp ../sltk_grid_driver.cpp ../sltk_grid.cpp ../sltk_atom.cpp ../../../source_io/module_output/output.cpp -) \ No newline at end of file +) + +AddTest( + TARGET MODULE_CELL_NEIGHBOR_sltk_bruteforce + LIBS parameter ${math_libs} base device cell_info + SOURCES sltk_bruteforce_test.cpp ../sltk_grid_driver.cpp ../sltk_grid.cpp + ../sltk_atom.cpp + ../../../source_io/module_output/output.cpp +) diff --git a/source/source_cell/module_neighbor/test/sltk_bruteforce_test.cpp b/source/source_cell/module_neighbor/test/sltk_bruteforce_test.cpp new file mode 100644 index 00000000000..f8dd5e1eec1 --- /dev/null +++ b/source/source_cell/module_neighbor/test/sltk_bruteforce_test.cpp @@ -0,0 +1,645 @@ +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +#include "../sltk_grid_driver.h" +#include "prepare_unitcell.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef __LCAO +InfoNonlocal::InfoNonlocal() +{ +} +InfoNonlocal::~InfoNonlocal() +{ +} +LCAO_Orbitals::LCAO_Orbitals() +{ +} +LCAO_Orbitals::~LCAO_Orbitals() +{ +} +#endif + +Magnetism::Magnetism() +{ + this->tot_mag = 0.0; + this->abs_mag = 0.0; + this->start_mag = nullptr; +} +Magnetism::~Magnetism() +{ + delete[] this->start_mag; +} + +namespace +{ +struct NeighborKey +{ + int type = 0; + int atom = 0; + int box_x = 0; + int box_y = 0; + int box_z = 0; + + bool operator<(const NeighborKey& other) const + { + return std::tie(type, atom, box_x, box_y, box_z) + < std::tie(other.type, other.atom, other.box_x, other.box_y, other.box_z); + } + + bool operator==(const NeighborKey& other) const + { + return std::tie(type, atom, box_x, box_y, box_z) + == std::tie(other.type, other.atom, other.box_x, other.box_y, other.box_z); + } +}; + +struct ExpandLayers +{ + int plus_x = 0; + int minus_x = 0; + int plus_y = 0; + int minus_y = 0; + int plus_z = 0; + int minus_z = 0; +}; + +ModuleBase::Matrix3 MatrixFromRows(const std::vector& rows) +{ + ModuleBase::Matrix3 latvec; + latvec.e11 = rows[0]; + latvec.e12 = rows[1]; + latvec.e13 = rows[2]; + latvec.e21 = rows[3]; + latvec.e22 = rows[4]; + latvec.e23 = rows[5]; + latvec.e31 = rows[6]; + latvec.e32 = rows[7]; + latvec.e33 = rows[8]; + return latvec; +} + +std::valarray ValarrayFromVector(const std::vector& values) +{ + return std::valarray(values.data(), values.size()); +} + +std::valarray ValarrayFromVector(const std::vector& values) +{ + return std::valarray(values.data(), values.size()); +} + +std::unique_ptr MakeCell(const double lat0, + const ModuleBase::Matrix3& latvec, + const std::vector& labels, + const std::vector& atoms_per_type, + const std::vector& coords) +{ + std::vector pp_files; + std::vector pp_types; + std::vector orb_files; + std::vector masses; + for (const std::string& label: labels) + { + pp_files.push_back(label + ".upf"); + pp_types.push_back("upf201"); + orb_files.push_back(label + ".orb"); + masses.push_back(1.0); + } + + const std::vector latvec_rows = {latvec.e11, + latvec.e12, + latvec.e13, + latvec.e21, + latvec.e22, + latvec.e23, + latvec.e31, + latvec.e32, + latvec.e33}; + UcellTestPrepare prep("custom", + 2, + false, + false, + false, + "volume", + lat0, + ValarrayFromVector(latvec_rows), + labels, + pp_files, + pp_types, + orb_files, + ValarrayFromVector(atoms_per_type), + masses, + "Cartesian", + ValarrayFromVector(coords)); + return std::unique_ptr(prep.SetUcellInfo()); +} + +double CrossNorm(const ModuleBase::Vector3& lhs, const ModuleBase::Vector3& rhs) +{ + const double x = lhs.y * rhs.z - lhs.z * rhs.y; + const double y = lhs.z * rhs.x - lhs.x * rhs.z; + const double z = lhs.x * rhs.y - lhs.y * rhs.x; + return std::sqrt(x * x + y * y + z * z); +} + +ExpandLayers BruteForceLayers(const UnitCell& ucell, const double radius_lat0) +{ + const ModuleBase::Vector3 a1(ucell.latvec.e11, ucell.latvec.e12, ucell.latvec.e13); + const ModuleBase::Vector3 a2(ucell.latvec.e21, ucell.latvec.e22, ucell.latvec.e23); + const ModuleBase::Vector3 a3(ucell.latvec.e31, ucell.latvec.e32, ucell.latvec.e33); + const double lat0_volume = ucell.lat0 * ucell.lat0 * ucell.lat0; + + const int dx = std::ceil(CrossNorm(a2, a3) * radius_lat0 / ucell.omega * lat0_volume); + const int dy = std::ceil(CrossNorm(a3, a1) * radius_lat0 / ucell.omega * lat0_volume); + const int dz = std::ceil(CrossNorm(a1, a2) * radius_lat0 / ucell.omega * lat0_volume); + + return {dx + 1, dx, dy + 1, dy, dz + 1, dz}; +} + +std::vector BruteForceNeighbors(const UnitCell& ucell, + const double radius_lat0, + const int center_type, + const int center_atom) +{ + const ExpandLayers layers = BruteForceLayers(ucell, radius_lat0); + const ModuleBase::Vector3 a1(ucell.latvec.e11, ucell.latvec.e12, ucell.latvec.e13); + const ModuleBase::Vector3 a2(ucell.latvec.e21, ucell.latvec.e22, ucell.latvec.e23); + const ModuleBase::Vector3 a3(ucell.latvec.e31, ucell.latvec.e32, ucell.latvec.e33); + const ModuleBase::Vector3& center = ucell.atoms[center_type].tau[center_atom]; + const double radius2 = radius_lat0 * radius_lat0; + + std::vector result; + for (int ix = -layers.minus_x; ix < layers.plus_x; ++ix) + { + for (int iy = -layers.minus_y; iy < layers.plus_y; ++iy) + { + for (int iz = -layers.minus_z; iz < layers.plus_z; ++iz) + { + const ModuleBase::Vector3 shift(ix * a1.x + iy * a2.x + iz * a3.x, + ix * a1.y + iy * a2.y + iz * a3.y, + ix * a1.z + iy * a2.z + iz * a3.z); + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + const ModuleBase::Vector3 image(ucell.atoms[it].tau[ia].x + shift.x, + ucell.atoms[it].tau[ia].y + shift.y, + ucell.atoms[it].tau[ia].z + shift.z); + const double dx = center.x - image.x; + const double dy = center.y - image.y; + const double dz = center.z - image.z; + const double dr2 = dx * dx + dy * dy + dz * dz; + const bool central_self + = (it == center_type && ia == center_atom && ix == 0 && iy == 0 && iz == 0); + if (!central_self && dr2 <= radius2) + { + result.push_back({it, ia, ix, iy, iz}); + } + } + } + } + } + } + std::sort(result.begin(), result.end()); + return result; +} + +std::vector GridNeighbors(const Grid_Driver& grid, + const UnitCell& ucell, + const int center_type, + const int center_atom) +{ + AdjacentAtomInfo adjs; + grid.Find_atom(ucell, center_type, center_atom, &adjs); + EXPECT_EQ(adjs.ntype.size(), static_cast(adjs.adj_num + 1)); + EXPECT_EQ(adjs.natom.size(), static_cast(adjs.adj_num + 1)); + EXPECT_EQ(adjs.box.size(), static_cast(adjs.adj_num + 1)); + EXPECT_EQ(adjs.adjacent_tau.size(), static_cast(adjs.adj_num + 1)); + EXPECT_EQ(adjs.ntype[adjs.adj_num], center_type); + EXPECT_EQ(adjs.natom[adjs.adj_num], center_atom); + EXPECT_EQ(adjs.box[adjs.adj_num].x, 0); + EXPECT_EQ(adjs.box[adjs.adj_num].y, 0); + EXPECT_EQ(adjs.box[adjs.adj_num].z, 0); + + std::vector result; + for (int i = 0; i < adjs.adj_num; ++i) + { + result.push_back({adjs.ntype[i], adjs.natom[i], adjs.box[i].x, adjs.box[i].y, adjs.box[i].z}); + } + std::sort(result.begin(), result.end()); + return result; +} + +void ExpectGridMatchesBruteForce(const UnitCell& ucell, + const double radius_lat0, + const std::string& output_name) +{ + Grid_Driver grid(false, false); + std::ofstream ofs(output_name); + grid.init(ofs, ucell, radius_lat0, true); + ofs.close(); + + const ExpandLayers layers = BruteForceLayers(ucell, radius_lat0); + EXPECT_EQ(grid.getGlayerX(), layers.plus_x); + EXPECT_EQ(grid.getGlayerX_minus(), layers.minus_x); + EXPECT_EQ(grid.getGlayerY(), layers.plus_y); + EXPECT_EQ(grid.getGlayerY_minus(), layers.minus_y); + EXPECT_EQ(grid.getGlayerZ(), layers.plus_z); + EXPECT_EQ(grid.getGlayerZ_minus(), layers.minus_z); + + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + EXPECT_EQ(GridNeighbors(grid, ucell, it, ia), BruteForceNeighbors(ucell, radius_lat0, it, ia)) + << "center type=" << it << " atom=" << ia; + } + } + std::remove(output_name.c_str()); +} + +bool ContainsNeighbor(const std::vector& neighbors, const NeighborKey& key) +{ + return std::find(neighbors.begin(), neighbors.end(), key) != neighbors.end(); +} + +void ExpectReciprocalGridNeighbors(const UnitCell& ucell, + const double radius_lat0, + const std::string& output_name) +{ + Grid_Driver grid(false, false); + std::ofstream ofs(output_name); + grid.init(ofs, ucell, radius_lat0, true); + ofs.close(); + + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + const std::vector neighbors = GridNeighbors(grid, ucell, it, ia); + for (const NeighborKey& neighbor: neighbors) + { + const std::vector reverse_neighbors = GridNeighbors(grid, ucell, neighbor.type, neighbor.atom); + const NeighborKey reverse_key{it, ia, -neighbor.box_x, -neighbor.box_y, -neighbor.box_z}; + EXPECT_TRUE(ContainsNeighbor(reverse_neighbors, reverse_key)) + << "missing reverse neighbor for center type=" << it << " atom=" << ia; + } + } + } + std::remove(output_name.c_str()); +} +} // namespace + +TEST(SltkBruteForceTest, OrthogonalMultiTypeCellMatchesIndependentBruteForce) +{ + const auto ucell = MakeCell(20.0, + MatrixFromRows({1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}), + {"A", "B"}, + {2, 2}, + {0.05, 0.05, 0.05, + 0.55, 0.50, 0.50, + 0.95, 0.05, 0.05, + 0.12, 0.06, 0.05}); + + ExpectGridMatchesBruteForce(*ucell, 0.18, "bruteforce_orthogonal.out"); +} + +TEST(SltkBruteForceTest, PbcFaceEdgeAndCornerImagesMatchIndependentBruteForce) +{ + const auto ucell = MakeCell(20.0, + MatrixFromRows({1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}), + {"A"}, + {4}, + {0.02, 0.02, 0.02, + 0.98, 0.02, 0.02, + 0.98, 0.98, 0.02, + 0.98, 0.98, 0.98}); + + ExpectGridMatchesBruteForce(*ucell, 0.08, "bruteforce_pbc_edge.out"); +} + +TEST(SltkBruteForceTest, TriclinicCellMatchesIndependentBruteForce) +{ + const auto ucell = MakeCell(25.0, + MatrixFromRows({1.0, 0.22, 0.07, 0.15, 1.0, 0.18, 0.08, 0.25, 1.0}), + {"A", "B"}, + {2, 1}, + {0.03, 0.96, 0.04, + 0.52, 0.48, 0.51, + 0.94, 0.08, 0.06}); + + ExpectGridMatchesBruteForce(*ucell, 0.24, "bruteforce_triclinic.out"); +} + +TEST(SltkBruteForceTest, CutoffBoundaryIncludesEqualDistanceAndRejectsOutside) +{ + const auto ucell = MakeCell(20.0, + MatrixFromRows({1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}), + {"A", "B"}, + {2, 2}, + {0.50, 0.50, 0.50, + 0.75, 0.50, 0.50, + 0.25, 0.50, 0.50, + 0.751, 0.50, 0.50}); + + ExpectGridMatchesBruteForce(*ucell, 0.25, "bruteforce_cutoff_boundary.out"); +} + +TEST(SltkBruteForceTest, DenseClusterAcrossPeriodicBoundaryKeepsImageIdentities) +{ + const auto ucell = MakeCell(30.0, + MatrixFromRows({1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}), + {"A", "B", "C"}, + {3, 2, 2}, + {0.015, 0.020, 0.025, + 0.985, 0.020, 0.025, + 0.020, 0.985, 0.025, + 0.985, 0.985, 0.025, + 0.030, 0.025, 0.985, + 0.975, 0.975, 0.975, + 0.500, 0.500, 0.500}); + + ExpectGridMatchesBruteForce(*ucell, 0.09, "bruteforce_dense_boundary_cluster.out"); +} + +TEST(SltkBruteForceTest, StronglySkewedTriclinicCellMatchesIndependentBruteForce) +{ + const auto ucell = MakeCell(18.0, + MatrixFromRows({1.0, 0.46, 0.19, 0.38, 1.0, 0.31, 0.12, 0.27, 1.0}), + {"A", "B"}, + {3, 3}, + {0.06, 0.08, 0.91, + 0.18, 0.88, 0.12, + 0.43, 0.47, 0.51, + 0.92, 0.09, 0.08, + 0.78, 0.82, 0.87, + 0.51, 0.49, 0.48}); + + ExpectGridMatchesBruteForce(*ucell, 0.22, "bruteforce_strongly_skewed.out"); +} + +TEST(SltkBruteForceTest, DeterministicIrregularSmallCellMatchesIndependentBruteForce) +{ + const auto ucell = MakeCell(22.0, + MatrixFromRows({1.0, 0.17, 0.03, 0.09, 1.0, 0.21, 0.04, 0.13, 1.0}), + {"A", "B", "C"}, + {2, 3, 2}, + {0.137, 0.941, 0.219, + 0.773, 0.113, 0.617, + 0.951, 0.887, 0.071, + 0.064, 0.481, 0.902, + 0.382, 0.284, 0.336, + 0.621, 0.742, 0.158, + 0.229, 0.059, 0.809}); + + ExpectGridMatchesBruteForce(*ucell, 0.19, "bruteforce_irregular_small.out"); +} + +TEST(SltkBruteForceTest, SameAtomPeriodicImagesAcrossMultipleShellsAreKept) +{ + const auto ucell = MakeCell(12.0, + MatrixFromRows({1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}), + {"A"}, + {1}, + {0.50, 0.50, 0.50}); + + ExpectGridMatchesBruteForce(*ucell, 1.05, "bruteforce_same_atom_images.out"); +} + +TEST(SltkBruteForceTest, AnisotropicCellUsesDifferentExpansionLayersPerDirection) +{ + const auto ucell = MakeCell(16.0, + MatrixFromRows({0.55, 0.00, 0.00, 0.00, 1.70, 0.00, 0.00, 0.00, 1.25}), + {"A", "B"}, + {2, 2}, + {0.08, 0.10, 0.12, + 0.91, 0.11, 0.13, + 0.49, 1.48, 0.60, + 0.52, 0.24, 1.10}); + + ExpectGridMatchesBruteForce(*ucell, 0.34, "bruteforce_anisotropic.out"); +} + +TEST(SltkBruteForceTest, NearCutoffAcrossPeriodicFaceKeepsOnlyInsideImage) +{ + const auto ucell = MakeCell(10.0, + MatrixFromRows({1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}), + {"A"}, + {4}, + {0.0100, 0.5000, 0.5000, + 0.9899, 0.5000, 0.5000, + 0.6500, 0.5000, 0.5000, + 0.6502, 0.5000, 0.5000}); + + ExpectGridMatchesBruteForce(*ucell, 0.35, "bruteforce_periodic_cutoff_face.out"); +} + +TEST(SltkBruteForceTest, ManyTypesWithUnevenAtomCountsPreserveTypeAndAtomIdentity) +{ + const auto ucell = MakeCell(24.0, + MatrixFromRows({1.0, 0.11, 0.02, 0.05, 1.0, 0.16, 0.03, 0.08, 1.0}), + {"A", "B", "C", "D"}, + {1, 3, 1, 2}, + {0.04, 0.05, 0.06, + 0.96, 0.05, 0.06, + 0.50, 0.50, 0.50, + 0.21, 0.78, 0.32, + 0.08, 0.94, 0.07, + 0.92, 0.91, 0.93, + 0.43, 0.18, 0.88}); + + ExpectGridMatchesBruteForce(*ucell, 0.17, "bruteforce_many_types_uneven.out"); +} + +TEST(SltkBruteForceTest, VacuumLikeSparseCellDoesNotInventNeighbors) +{ + const auto ucell = MakeCell(80.0, + MatrixFromRows({1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}), + {"A", "B"}, + {2, 1}, + {0.10, 0.10, 0.10, + 0.80, 0.80, 0.80, + 0.45, 0.45, 0.45}); + + ExpectGridMatchesBruteForce(*ucell, 0.05, "bruteforce_vacuum_sparse.out"); +} + +TEST(SltkBruteForceTest, LeadingEmptyAtomTypeDoesNotBreakIndexing) +{ + const auto ucell = MakeCell(20.0, + MatrixFromRows({1.0, 0.08, 0.02, 0.03, 1.0, 0.06, 0.01, 0.04, 1.0}), + {"Empty", "A", "B"}, + {0, 2, 1}, + {0.05, 0.05, 0.05, + 0.95, 0.05, 0.05, + 0.52, 0.50, 0.50}); + + ExpectGridMatchesBruteForce(*ucell, 0.16, "bruteforce_leading_empty_type.out"); +} + +TEST(SltkBruteForceTest, VerySmallNonzeroSeparationIsNotMistakenForSelf) +{ + const auto ucell = MakeCell(20.0, + MatrixFromRows({1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}), + {"A", "B"}, + {1, 1}, + {0.40, 0.40, 0.40, + 0.400001, 0.40, 0.40}); + + ExpectGridMatchesBruteForce(*ucell, 0.00001, "bruteforce_tiny_separation.out"); +} + +TEST(SltkBruteForceTest, SubAngstromScaleSeparationsRemainDistinctDownToDoublePrecision) +{ + const std::vector separations = {1.0e-8, 1.0e-10, 1.0e-12, 1.0e-14, 1.0e-15}; + for (const double delta: separations) + { + const auto ucell = MakeCell(20.0, + MatrixFromRows({1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}), + {"A", "B"}, + {1, 1}, + {0.40, 0.40, 0.40, + 0.40 + delta, 0.40, 0.40}); + + ExpectGridMatchesBruteForce(*ucell, delta * 2.0, "bruteforce_double_precision_separation.out"); + } +} + +TEST(SltkBruteForceTest, DirectionalTinySeparationsInSkewedCellRemainDistinct) +{ + const std::vector separations = {1.0e-8, 1.0e-10, 1.0e-12, 1.0e-14}; + for (const double delta: separations) + { + const auto ucell = MakeCell(20.0, + MatrixFromRows({1.0, 0.31, 0.07, 0.18, 1.0, 0.23, 0.05, 0.29, 1.0}), + {"A"}, + {8}, + {0.40, 0.40, 0.40, + 0.40 + delta, 0.40, 0.40, + 0.40, 0.40 + delta, 0.40, + 0.40, 0.40, 0.40 + delta, + 0.40 + delta, 0.40 + delta, 0.40, + 0.40 + delta, 0.40, 0.40 + delta, + 0.40, 0.40 + delta, 0.40 + delta, + 0.40 + delta, 0.40 + delta, 0.40 + delta}); + + ExpectGridMatchesBruteForce(*ucell, delta * 3.0, "bruteforce_directional_tiny_separation.out"); + } +} + +TEST(SltkBruteForceTest, DeterministicFuzzSmallCellsMatchIndependentBruteForce) +{ + const std::vector lattices = { + MatrixFromRows({1.0, 0.00, 0.00, 0.00, 1.0, 0.00, 0.00, 0.00, 1.0}), + MatrixFromRows({1.0, 0.18, 0.04, 0.07, 1.0, 0.11, 0.03, 0.09, 1.0}), + MatrixFromRows({0.82, 0.28, 0.05, 0.16, 1.13, 0.18, 0.04, 0.22, 0.91}), + }; + const std::vector radii = {0.13, 0.19, 0.27}; + + for (size_t ilat = 0; ilat < lattices.size(); ++ilat) + { + for (size_t ir = 0; ir < radii.size(); ++ir) + { + std::vector coords; + coords.reserve(18); + for (int i = 0; i < 6; ++i) + { + const double seed = static_cast((ilat + 1) * 37 + (ir + 1) * 19 + i * 23); + coords.push_back(std::fmod(0.137 * seed + 0.071 * i, 0.92) + 0.04); + coords.push_back(std::fmod(0.193 * seed + 0.053 * i, 0.92) + 0.04); + coords.push_back(std::fmod(0.271 * seed + 0.029 * i, 0.92) + 0.04); + } + + const auto ucell = MakeCell(24.0, + lattices[ilat], + {"A", "B", "C"}, + {2, 3, 1}, + coords); + ExpectGridMatchesBruteForce(*ucell, radii[ir], "bruteforce_deterministic_fuzz.out"); + } + } +} + +TEST(SltkBruteForceTest, NearZeroAndOneFractionalBoundariesMatchIndependentBruteForce) +{ + const double eps = 1.0e-12; + const auto ucell = MakeCell(20.0, + MatrixFromRows({1.0, 0.16, 0.03, 0.05, 1.0, 0.12, 0.02, 0.07, 1.0}), + {"A", "B"}, + {4, 4}, + {eps, 0.50, 0.50, + 1.0 - eps, 0.50, 0.50, + 0.50, eps, 0.50, + 0.50, 1.0 - eps, 0.50, + 0.50, 0.50, eps, + 0.50, 0.50, 1.0 - eps, + eps, eps, eps, + 1.0 - eps, 1.0 - eps, 1.0 - eps}); + + ExpectGridMatchesBruteForce(*ucell, 0.09, "bruteforce_near_boundary_normalization.out"); +} + +TEST(SltkBruteForceTest, ThinStronglySkewedCellMatchesIndependentBruteForce) +{ + const auto ucell = MakeCell(18.0, + MatrixFromRows({1.0, 0.65, 0.00, 0.15, 1.0, 0.00, 0.03, 0.04, 0.22}), + {"A", "B"}, + {3, 3}, + {0.05, 0.06, 0.08, + 0.95, 0.07, 0.09, + 0.48, 0.52, 0.11, + 0.08, 0.91, 0.17, + 0.89, 0.86, 0.19, + 0.51, 0.49, 0.03}); + + ExpectGridMatchesBruteForce(*ucell, 0.11, "bruteforce_thin_strongly_skewed.out"); +} + +TEST(SltkBruteForceTest, NeighborRelationIsReciprocalWithOppositePeriodicBox) +{ + const auto ucell = MakeCell(18.0, + MatrixFromRows({1.0, 0.25, 0.05, 0.12, 1.0, 0.17, 0.04, 0.20, 1.0}), + {"A", "B", "C"}, + {2, 2, 1}, + {0.04, 0.05, 0.06, + 0.96, 0.07, 0.06, + 0.10, 0.91, 0.08, + 0.88, 0.86, 0.92, + 0.52, 0.48, 0.50}); + + ExpectReciprocalGridNeighbors(*ucell, 0.19, "bruteforce_reciprocal.out"); +} + +TEST(SltkBruteForceTest, DISABLED_PeriodicBoundaryEquivalentPositionsShouldRemainDistinctNeighbors) +{ + const auto ucell = MakeCell(20.0, + MatrixFromRows({1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}), + {"A", "B"}, + {1, 1}, + {0.0, 0.50, 0.50, + 1.0, 0.50, 0.50}); + + ExpectGridMatchesBruteForce(*ucell, 0.10, "bruteforce_periodic_equivalent_positions.out"); +} + +TEST(SltkBruteForceTest, DISABLED_DistinctAtomsAtSamePositionShouldBeReportedAsNeighbors) +{ + const auto ucell = MakeCell(20.0, + MatrixFromRows({1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}), + {"A", "B"}, + {1, 1}, + {0.40, 0.40, 0.40, + 0.40, 0.40, 0.40}); + + ExpectGridMatchesBruteForce(*ucell, 0.10, "bruteforce_overlapping_atoms.out"); +} diff --git a/source/source_cell/module_neighbor/test/sltk_grid_test.cpp b/source/source_cell/module_neighbor/test/sltk_grid_test.cpp index 044feafc2de..374643de02a 100644 --- a/source/source_cell/module_neighbor/test/sltk_grid_test.cpp +++ b/source/source_cell/module_neighbor/test/sltk_grid_test.cpp @@ -6,6 +6,12 @@ #include "prepare_unitcell.h" #include "source_io/module_parameter/parameter.h" #undef private +#include +#include +#include +#include +#include +#include #include "source_cell/read_stru.h" #ifdef __LCAO InfoNonlocal::InfoNonlocal() @@ -48,6 +54,7 @@ Magnetism::~Magnetism() void SetGlobalV() { PARAM.input.test_grid = 0; + PARAM.input.test_deconstructor = 0; } class SltkGridTest : public testing::Test @@ -74,6 +81,104 @@ class SltkGridTest : public testing::Test using SltkGridDeathTest = SltkGridTest; +namespace +{ +void SetLattice(UnitCell* ucell, + const double lat0, + const ModuleBase::Matrix3& latvec) +{ + ucell->lat0 = lat0; + ucell->lat0_angstrom = lat0 * ModuleBase::BOHR_TO_A; + ucell->tpiba = ModuleBase::TWO_PI / lat0; + ucell->tpiba2 = ucell->tpiba * ucell->tpiba; + ucell->latvec = latvec; + ucell->a1.set(latvec.e11, latvec.e12, latvec.e13); + ucell->a2.set(latvec.e21, latvec.e22, latvec.e23); + ucell->a3.set(latvec.e31, latvec.e32, latvec.e33); + ucell->GT = latvec.Inverse(); + ucell->G = ucell->GT.Transpose(); + ucell->GGT = ucell->G * ucell->GT; + ucell->invGGT = ucell->GGT.Inverse(); + ucell->omega = std::abs(latvec.Det()) * lat0 * lat0 * lat0; +} + +ModuleBase::Matrix3 IdentityLatvec() +{ + ModuleBase::Matrix3 latvec; + latvec.e11 = 1.0; latvec.e12 = 0.0; latvec.e13 = 0.0; + latvec.e21 = 0.0; latvec.e22 = 1.0; latvec.e23 = 0.0; + latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = 1.0; + return latvec; +} + +ModuleBase::Matrix3 SkewedLatvec() +{ + ModuleBase::Matrix3 latvec; + latvec.e11 = 1.0; latvec.e12 = 0.3; latvec.e13 = 0.0; + latvec.e21 = 0.1; latvec.e22 = 1.0; latvec.e23 = 0.0; + latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = 1.0; + return latvec; +} + +void SetTwoAtomPositions(UnitCell* ucell, + const ModuleBase::Vector3& tau0, + const ModuleBase::Vector3& tau1) +{ + ASSERT_EQ(ucell->ntype, 1); + ASSERT_GE(ucell->atoms[0].na, 2); + ucell->atoms[0].tau[0] = tau0; + ucell->atoms[0].tau[1] = tau1; +} + + +std::size_t CountExpandedAtomsInBoxes(const Grid& grid) +{ + std::size_t count = 0; + for (const auto& boxes_yz: grid.atoms_in_box) + { + for (const auto& boxes_z: boxes_yz) + { + for (const auto& atoms: boxes_z) + { + count += atoms.size(); + } + } + } + return count; +} + +std::size_t ExpectedExpandedAtomCount(const Grid& grid, const UnitCell& ucell) +{ + const std::size_t layer_x = static_cast(grid.getGlayerX() + grid.getGlayerX_minus()); + const std::size_t layer_y = static_cast(grid.getGlayerY() + grid.getGlayerY_minus()); + const std::size_t layer_z = static_cast(grid.getGlayerZ() + grid.getGlayerZ_minus()); + return static_cast(ucell.nat) * layer_x * layer_y * layer_z; +} + +std::size_t SaturatedAllocatedBoxCount(const Grid& grid) +{ + if (grid.box_nx <= 0 || grid.box_ny <= 0 || grid.box_nz <= 0) + { + return 0; + } + const std::size_t sx = static_cast(grid.box_nx); + const std::size_t sy = static_cast(grid.box_ny); + const std::size_t sz = static_cast(grid.box_nz); + const std::size_t max_size = std::numeric_limits::max(); + if (sx > max_size / sy) + { + return max_size; + } + const std::size_t xy = sx * sy; + if (xy > max_size / sz) + { + return max_size; + } + return xy * sz; +} + +} // namespace + TEST_F(SltkGridTest, Init) { ofs.open("test.out"); @@ -88,6 +193,14 @@ TEST_F(SltkGridTest, Init) EXPECT_EQ(LatGrid.getGlayerX_minus(), 5); EXPECT_EQ(LatGrid.getGlayerY_minus(), 5); EXPECT_EQ(LatGrid.getGlayerZ_minus(), 5); + EXPECT_GE(LatGrid.box_edge_length, LatGrid.sradius); + EXPECT_GT(LatGrid.box_nx, 0); + EXPECT_GT(LatGrid.box_ny, 0); + EXPECT_GT(LatGrid.box_nz, 0); + EXPECT_EQ(LatGrid.atoms_in_box.size(), static_cast(LatGrid.box_nx)); + EXPECT_EQ(CountExpandedAtomsInBoxes(LatGrid), ExpectedExpandedAtomCount(LatGrid, *ucell)); + EXPECT_EQ(LatGrid.all_adj_info.size(), static_cast(ucell->ntype)); + ASSERT_EQ(LatGrid.all_adj_info[0].size(), static_cast(ucell->atoms[0].na)); ofs.close(); remove("test.out"); } @@ -108,6 +221,12 @@ TEST_F(SltkGridTest, InitSmall) EXPECT_DOUBLE_EQ(LatGrid.sradius2, 0.5 * 0.5); EXPECT_DOUBLE_EQ(LatGrid.sradius, radius); EXPECT_DOUBLE_EQ(LatGrid.sradius, 0.5); + EXPECT_GT(LatGrid.box_nx, 0); + EXPECT_GT(LatGrid.box_ny, 0); + EXPECT_GT(LatGrid.box_nz, 0); + EXPECT_EQ(LatGrid.atoms_in_box.size(), static_cast(LatGrid.box_nx)); + ASSERT_FALSE(LatGrid.atoms_in_box.empty()); + EXPECT_EQ(LatGrid.atoms_in_box[0].size(), static_cast(LatGrid.box_ny)); /* // minimal value of x, y, z EXPECT_DOUBLE_EQ(LatGrid.true_cell_x, 1); @@ -122,6 +241,147 @@ TEST_F(SltkGridTest, InitSmall) remove("test.out"); } +TEST_F(SltkGridTest, InitOrthogonalCellSetsExpectedLayersAndBoxCounts) +{ + SetLattice(ucell, 1.0, IdentityLatvec()); + SetTwoAtomPositions(ucell, + ModuleBase::Vector3(0.0, 0.0, 0.0), + ModuleBase::Vector3(0.5, 0.0, 0.0)); + + ofs.open("orthogonal_grid.out"); + Grid lat_grid(PARAM.input.test_grid); + lat_grid.init(ofs, *ucell, 1.0, true); + + EXPECT_EQ(lat_grid.getGlayerX(), 2); + EXPECT_EQ(lat_grid.getGlayerY(), 2); + EXPECT_EQ(lat_grid.getGlayerZ(), 2); + EXPECT_EQ(lat_grid.getGlayerX_minus(), 1); + EXPECT_EQ(lat_grid.getGlayerY_minus(), 1); + EXPECT_EQ(lat_grid.getGlayerZ_minus(), 1); + EXPECT_EQ(lat_grid.box_nx, 3); + EXPECT_EQ(lat_grid.box_ny, 3); + EXPECT_EQ(lat_grid.box_nz, 3); + EXPECT_DOUBLE_EQ(lat_grid.sradius, 1.0); + EXPECT_DOUBLE_EQ(lat_grid.sradius2, 1.0); + EXPECT_EQ(lat_grid.atoms_in_box.size(), 3u); + + ofs.close(); + remove("orthogonal_grid.out"); +} + +TEST_F(SltkGridTest, InitNonOrthogonalCellBuildsPositiveExpansion) +{ + SetLattice(ucell, 1.0, SkewedLatvec()); + SetTwoAtomPositions(ucell, + ModuleBase::Vector3(0.0, 0.0, 0.0), + ModuleBase::Vector3(0.45, 0.20, 0.0)); + + ofs.open("skewed_grid.out"); + Grid lat_grid(PARAM.input.test_grid); + lat_grid.init(ofs, *ucell, 0.75, true); + + EXPECT_GE(lat_grid.getGlayerX(), 1); + EXPECT_GE(lat_grid.getGlayerY(), 1); + EXPECT_GE(lat_grid.getGlayerZ(), 1); + EXPECT_GE(lat_grid.box_edge_length, lat_grid.sradius); + EXPECT_GT(lat_grid.box_nx, 0); + EXPECT_GT(lat_grid.box_ny, 0); + EXPECT_GT(lat_grid.box_nz, 0); + EXPECT_GT(lat_grid.atoms_in_box.size(), 0u); + EXPECT_EQ(CountExpandedAtomsInBoxes(lat_grid), ExpectedExpandedAtomCount(lat_grid, *ucell)); + EXPECT_EQ(lat_grid.all_adj_info.size(), static_cast(ucell->ntype)); + + ofs.close(); + remove("skewed_grid.out"); +} + +TEST_F(SltkGridTest, InitWithPbcDisabledKeepsNoPeriodicExpansion) +{ + SetLattice(ucell, 1.0, IdentityLatvec()); + SetTwoAtomPositions(ucell, + ModuleBase::Vector3(0.0, 0.0, 0.0), + ModuleBase::Vector3(0.5, 0.0, 0.0)); + + ofs.open("no_pbc_grid.out"); + Grid lat_grid(PARAM.input.test_grid); + lat_grid.init(ofs, *ucell, 1.0, false); + + EXPECT_FALSE(lat_grid.pbc); + EXPECT_EQ(lat_grid.getGlayerX(), 0); + EXPECT_EQ(lat_grid.getGlayerY(), 0); + EXPECT_EQ(lat_grid.getGlayerZ(), 0); + EXPECT_EQ(lat_grid.getGlayerX_minus(), 0); + EXPECT_EQ(lat_grid.getGlayerY_minus(), 0); + EXPECT_EQ(lat_grid.getGlayerZ_minus(), 0); + EXPECT_EQ(lat_grid.box_nx, 0); + EXPECT_EQ(lat_grid.box_ny, 0); + EXPECT_EQ(lat_grid.box_nz, 0); + EXPECT_EQ(lat_grid.all_adj_info.size(), static_cast(ucell->ntype)); + + ofs.close(); + remove("no_pbc_grid.out"); +} + +TEST_F(SltkGridTest, CellListCoarsensBinsForSparseHugeAabb) +{ + ModuleBase::Matrix3 large_latvec; + large_latvec.e11 = 100000.0; large_latvec.e12 = 0.0; large_latvec.e13 = 0.0; + large_latvec.e21 = 0.0; large_latvec.e22 = 100000.0; large_latvec.e23 = 0.0; + large_latvec.e31 = 0.0; large_latvec.e32 = 0.0; large_latvec.e33 = 100000.0; + SetLattice(ucell, 1.0, large_latvec); + SetTwoAtomPositions(ucell, + ModuleBase::Vector3(0.0, 0.0, 0.0), + ModuleBase::Vector3(0.5, 0.5, 0.5)); + + ofs.open("sparse_huge_aabb_grid.out"); + Grid lat_grid(PARAM.input.test_grid); + lat_grid.init(ofs, *ucell, 0.01, true); + + EXPECT_GT(lat_grid.box_edge_length, 0.1); + EXPECT_GT(lat_grid.box_nx, 0); + EXPECT_GT(lat_grid.box_ny, 0); + EXPECT_GT(lat_grid.box_nz, 0); + EXPECT_LE(SaturatedAllocatedBoxCount(lat_grid), 1000000u); + EXPECT_EQ(CountExpandedAtomsInBoxes(lat_grid), ExpectedExpandedAtomCount(lat_grid, *ucell)); + + ofs.close(); + remove("sparse_huge_aabb_grid.out"); +} + +TEST_F(SltkGridTest, ReinitWithPbcDisabledResetsExpansionState) +{ + SetLattice(ucell, 1.0, IdentityLatvec()); + SetTwoAtomPositions(ucell, + ModuleBase::Vector3(0.0, 0.0, 0.0), + ModuleBase::Vector3(0.5, 0.0, 0.0)); + + Grid lat_grid(PARAM.input.test_grid); + ofs.open("reinit_pbc_grid.out"); + lat_grid.init(ofs, *ucell, 1.0, true); + ofs.close(); + EXPECT_GT(lat_grid.getGlayerX(), 0); + EXPECT_GT(lat_grid.box_nx, 0); + + ofs.open("reinit_no_pbc_grid.out"); + lat_grid.init(ofs, *ucell, 1.0, false); + ofs.close(); + + EXPECT_FALSE(lat_grid.pbc); + EXPECT_EQ(lat_grid.getGlayerX(), 0); + EXPECT_EQ(lat_grid.getGlayerY(), 0); + EXPECT_EQ(lat_grid.getGlayerZ(), 0); + EXPECT_EQ(lat_grid.getGlayerX_minus(), 0); + EXPECT_EQ(lat_grid.getGlayerY_minus(), 0); + EXPECT_EQ(lat_grid.getGlayerZ_minus(), 0); + EXPECT_EQ(lat_grid.box_nx, 0); + EXPECT_EQ(lat_grid.box_ny, 0); + EXPECT_EQ(lat_grid.box_nz, 0); + EXPECT_EQ(CountExpandedAtomsInBoxes(lat_grid), 0u); + + remove("reinit_pbc_grid.out"); + remove("reinit_no_pbc_grid.out"); +} + /* // This test cannot pass because setAtomLinkArray() is unsuccessful // if expand_flag is false