From 1145b388def32dd6e8a7102cb9ca730a9de23f4a Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Mon, 3 Aug 2026 17:45:28 +0800 Subject: [PATCH 01/10] feat: add distributed MdCell domain decomposition --- source/Makefile.Objects | 1 - source/source_cell/CMakeLists.txt | 2 + .../source_cell/distributed_mdcell_reader.cpp | 324 ++++++++++ .../source_cell/distributed_mdcell_reader.h | 16 + source/source_cell/md_cell.cpp | 579 ++++++++++++++++++ source/source_cell/md_cell.h | 142 +++++ .../module_neighlist/CMakeLists.txt | 1 - .../module_neighlist/atom_provider.h | 74 --- .../module_neighlist/domain_decomposition.cpp | 141 ++++- .../module_neighlist/domain_decomposition.h | 14 +- .../source_cell/module_neighlist/local_atom.h | 20 +- .../module_neighlist/neighbor_atom.h | 8 +- .../module_neighlist/neighbor_search.cpp | 96 +-- .../module_neighlist/neighbor_search.h | 26 +- .../module_neighlist/neighbor_types.h | 1 - .../module_neighlist/test/CMakeLists.txt | 49 +- .../test/distributed_mdcell_reader_test.cpp | 134 ++++ .../test/md_cell_migrate_mpi_test.cpp | 101 +++ .../test/neighbor_search_mpi_benchmark.cpp | 413 ------------- .../test/neighbor_search_test.cpp | 279 +++++---- .../module_neighlist/unitcell_lite.cpp | 95 --- .../module_neighlist/unitcell_lite.h | 169 ----- source/source_cell/unitcell.h | 31 +- source/source_esolver/esolver_lj.cpp | 30 +- source/source_esolver/esolver_lj.h | 3 - source/source_md/test/CMakeLists.txt | 2 +- 26 files changed, 1730 insertions(+), 1021 deletions(-) create mode 100644 source/source_cell/distributed_mdcell_reader.cpp create mode 100644 source/source_cell/distributed_mdcell_reader.h create mode 100644 source/source_cell/md_cell.cpp create mode 100644 source/source_cell/md_cell.h delete mode 100644 source/source_cell/module_neighlist/atom_provider.h create mode 100644 source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp create mode 100644 source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp delete mode 100644 source/source_cell/module_neighlist/test/neighbor_search_mpi_benchmark.cpp delete mode 100644 source/source_cell/module_neighlist/unitcell_lite.cpp delete mode 100644 source/source_cell/module_neighlist/unitcell_lite.h diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 771396233dc..c5c503b33bb 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -429,7 +429,6 @@ OBJS_NEIGHBOR_SEARCH=neighbor_search.o\ bin_manager.o\ domain_decomposition.o\ page_allocator.o\ - unitcell_lite.o\ OBJS_ORBITAL=ORB_atomic.o\ diff --git a/source/source_cell/CMakeLists.txt b/source/source_cell/CMakeLists.txt index ca64ccf60df..3aa11e61dfa 100644 --- a/source/source_cell/CMakeLists.txt +++ b/source/source_cell/CMakeLists.txt @@ -15,6 +15,8 @@ add_library( read_pp_upf201.cpp read_pp_blps.cpp read_pp_vwr.cpp + distributed_mdcell_reader.cpp + md_cell.cpp unitcell.cpp read_atoms.cpp read_atoms_helper.cpp diff --git a/source/source_cell/distributed_mdcell_reader.cpp b/source/source_cell/distributed_mdcell_reader.cpp new file mode 100644 index 00000000000..314ac05eea7 --- /dev/null +++ b/source/source_cell/distributed_mdcell_reader.cpp @@ -0,0 +1,324 @@ +#include "source_cell/distributed_mdcell_reader.h" + +#include "source_base/constants.h" +#include "source_base/vector3.h" +#include "source_cell/md_cell.h" + +#ifdef __MPI +#include "source_cell/module_neighlist/domain_decomposition.h" +#endif + +#include +#include +#include +#include +#include +#include + +namespace +{ +struct StruMetadata +{ + double lat0; + double omega; + ModuleBase::Matrix3 latvec; + ModuleBase::Matrix3 gt; + std::vector labels; + std::vector masses; + MdStruMetadata stru_metadata; +}; + +std::string trim_copy(const std::string& value) +{ + std::size_t begin = 0; + while (begin < value.size() && std::isspace(static_cast(value[begin]))) + { + ++begin; + } + std::size_t end = value.size(); + while (end > begin && std::isspace(static_cast(value[end - 1]))) + { + --end; + } + return value.substr(begin, end - begin); +} + +std::string strip_comment(const std::string& line) +{ + const std::size_t pos = line.find('#'); + return trim_copy(pos == std::string::npos ? line : line.substr(0, pos)); +} + +std::string next_data_line(std::ifstream& ifs, const char* context) +{ + std::string line; + while (std::getline(ifs, line)) + { + line = strip_comment(line); + if (!line.empty()) + { + return line; + } + } + throw std::runtime_error(std::string("Unexpected EOF while reading ") + context + "."); +} + +void expect_keyword(std::ifstream& ifs, const char* keyword) +{ + const std::string line = next_data_line(ifs, keyword); + if (line != keyword) + { + throw std::runtime_error(std::string("Expected keyword '") + keyword + "', got '" + line + "'."); + } +} + +double parse_double(const std::string& token, const char* context) +{ + char* end = NULL; + const double value = std::strtod(token.c_str(), &end); + if (end == token.c_str() || *end != '\0') + { + throw std::runtime_error(std::string("Failed to parse double for ") + context + ": " + token); + } + return value; +} + +int parse_int(const std::string& token, const char* context) +{ + char* end = NULL; + const long value = std::strtol(token.c_str(), &end, 10); + if (end == token.c_str() || *end != '\0') + { + throw std::runtime_error(std::string("Failed to parse int for ") + context + ": " + token); + } + return static_cast(value); +} + +ModuleBase::Vector3 wrap_fractional(const ModuleBase::Vector3& frac) +{ + ModuleBase::Vector3 wrapped = frac; + wrapped.x -= std::floor(wrapped.x); + wrapped.y -= std::floor(wrapped.y); + wrapped.z -= std::floor(wrapped.z); + if (wrapped.x >= 1.0 - 1.0e-12 || wrapped.x < 1.0e-12) wrapped.x = 0.0; + if (wrapped.y >= 1.0 - 1.0e-12 || wrapped.y < 1.0e-12) wrapped.y = 0.0; + if (wrapped.z >= 1.0 - 1.0e-12 || wrapped.z < 1.0e-12) wrapped.z = 0.0; + return wrapped; +} + +StruMetadata parse_stru_metadata(std::ifstream& ifs) +{ + StruMetadata metadata; + metadata.lat0 = 1.0; + metadata.omega = 0.0; + + expect_keyword(ifs, "ATOMIC_SPECIES"); + while (true) + { + const std::streampos mark = ifs.tellg(); + const std::string line = next_data_line(ifs, "ATOMIC_SPECIES body"); + if (line == "LATTICE_CONSTANT") + { + ifs.seekg(mark); + break; + } + if (line == "NUMERICAL_ORBITAL") + { + for (std::size_t it = 0; it < metadata.stru_metadata.species.size(); ++it) + { + metadata.stru_metadata.species[it].orbital_file = next_data_line(ifs, "NUMERICAL_ORBITAL body"); + } + continue; + } + if (line == "NUMERICAL_DESCRIPTOR") + { + metadata.stru_metadata.descriptor_file = next_data_line(ifs, "NUMERICAL_DESCRIPTOR body"); + continue; + } + + std::istringstream iss(line); + std::string label; + std::string mass_token; + iss >> label >> mass_token; + if (label.empty() || mass_token.empty()) + { + throw std::runtime_error("Invalid ATOMIC_SPECIES line: " + line); + } + + metadata.labels.push_back(label); + metadata.masses.push_back(parse_double(mass_token, "atomic mass")); + MdStruSpecies species; + species.label = label; + species.mass = metadata.masses.back(); + iss >> species.pseudo_file >> species.pseudo_type; + metadata.stru_metadata.species.push_back(species); + } + + expect_keyword(ifs, "LATTICE_CONSTANT"); + metadata.lat0 = parse_double(next_data_line(ifs, "LATTICE_CONSTANT value"), "lattice constant"); + + expect_keyword(ifs, "LATTICE_VECTORS"); + for (int row = 0; row < 3; ++row) + { + std::istringstream iss(next_data_line(ifs, "LATTICE_VECTORS row")); + double x = 0.0; + double y = 0.0; + double z = 0.0; + iss >> x >> y >> z; + if (!iss) + { + throw std::runtime_error("Invalid LATTICE_VECTORS row."); + } + if (row == 0) { metadata.latvec.e11 = x; metadata.latvec.e12 = y; metadata.latvec.e13 = z; } + if (row == 1) { metadata.latvec.e21 = x; metadata.latvec.e22 = y; metadata.latvec.e23 = z; } + if (row == 2) { metadata.latvec.e31 = x; metadata.latvec.e32 = y; metadata.latvec.e33 = z; } + } + metadata.gt = metadata.latvec.Inverse(); + metadata.omega = std::abs(metadata.latvec.Det()) * metadata.lat0 * metadata.lat0 * metadata.lat0; + return metadata; +} + +std::vector read_owned_atoms(std::ifstream& ifs, + StruMetadata& metadata, + double cutoff_bohr, + double skin_bohr, + int& nat) +{ + int rank = 0; +#ifdef __MPI + DomainDecomposition decomposition; + decomposition.init(MPI_COMM_WORLD, metadata.latvec, metadata.lat0, cutoff_bohr, skin_bohr); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); +#endif + + expect_keyword(ifs, "ATOMIC_POSITIONS"); + const std::string coord_type = next_data_line(ifs, "ATOMIC_POSITIONS type"); + const bool is_cartesian = coord_type == "Cartesian"; + const bool is_direct = coord_type == "Direct"; + if (!is_cartesian && !is_direct) + { + throw std::runtime_error("Only Direct and Cartesian ATOMIC_POSITIONS are supported for LJ MD."); + } + + std::vector owned_atoms; + nat = 0; + for (std::size_t it = 0; it < metadata.labels.size(); ++it) + { + const std::string label = next_data_line(ifs, "atom label"); + if (label != metadata.labels[it]) + { + throw std::runtime_error("ATOMIC_POSITIONS label order does not match ATOMIC_SPECIES."); + } + std::istringstream magnetism(next_data_line(ifs, "magnetism")); + magnetism >> metadata.stru_metadata.species[it].start_mag; + const int nat_type = parse_int(next_data_line(ifs, "atom count"), "atom count"); + metadata.stru_metadata.species[it].atom_count = nat_type; + + for (int ia = 0; ia < nat_type; ++ia) + { + std::istringstream iss(next_data_line(ifs, "atom line")); + double c1 = 0.0; + double c2 = 0.0; + double c3 = 0.0; + iss >> c1 >> c2 >> c3; + if (!iss) + { + throw std::runtime_error("Invalid atomic coordinate line."); + } + + ModuleBase::Vector3 frac; + ModuleBase::Vector3 cart; + if (is_cartesian) + { + cart.set(c1, c2, c3); + frac = wrap_fractional(cart * metadata.gt); + cart = frac * metadata.latvec; + } + else + { + frac = wrap_fractional(ModuleBase::Vector3(c1, c2, c3)); + cart = frac * metadata.latvec; + } + + ModuleBase::Vector3 mbl(1, 1, 1); + ModuleBase::Vector3 vel(0.0, 0.0, 0.0); + std::string token; + while (iss >> token) + { + if (token == "m") + { + std::string mx; + std::string my; + std::string mz; + iss >> mx >> my >> mz; + if (!iss) throw std::runtime_error("Invalid move flag record in STRU."); + mbl.set(parse_int(mx, "move flag x"), parse_int(my, "move flag y"), parse_int(mz, "move flag z")); + } + else if (token == "v" || token == "vel" || token == "velocity") + { + std::string vx; + std::string vy; + std::string vz; + iss >> vx >> vy >> vz; + if (!iss) throw std::runtime_error("Invalid velocity record in STRU."); + vel.set(parse_double(vx, "velocity x"), + parse_double(vy, "velocity y"), + parse_double(vz, "velocity z")); + } + } + + int owner = 0; +#ifdef __MPI + owner = decomposition.owner_rank_from_frac(frac); +#endif + if (owner == rank) + { + owned_atoms.push_back(LocalAtom(cart, + frac, + vel, + ModuleBase::Vector3(0.0, 0.0, 0.0), + mbl, + metadata.masses[it] / ModuleBase::AU_to_MASS, + static_cast(it), + ia, + owner, + false)); + } + ++nat; + } + } + return owned_atoms; +} +} // namespace + +MdCell DistributedMdCellReader::read_lj_stru(const std::string& stru_file, + double cutoff_bohr, + double skin_bohr) +{ + if (cutoff_bohr <= 0.0) + { + throw std::runtime_error("MdCell requires a positive LJ cutoff from Parameter."); + } + + std::ifstream ifs(stru_file.c_str(), std::ios::in); + if (!ifs) + { + throw std::runtime_error("Failed to open STRU file: " + stru_file); + } + + StruMetadata metadata = parse_stru_metadata(ifs); + int nat = 0; + const std::vector owned_atoms = read_owned_atoms(ifs, metadata, cutoff_bohr, skin_bohr, nat); + MdCell mdcell(metadata.latvec, + metadata.gt, + metadata.lat0, + metadata.omega, + nat, + owned_atoms, + metadata.labels, + metadata.masses, + cutoff_bohr, + skin_bohr); + mdcell.set_stru_metadata(metadata.stru_metadata); + return mdcell; +} diff --git a/source/source_cell/distributed_mdcell_reader.h b/source/source_cell/distributed_mdcell_reader.h new file mode 100644 index 00000000000..1af6fa66708 --- /dev/null +++ b/source/source_cell/distributed_mdcell_reader.h @@ -0,0 +1,16 @@ +#ifndef DISTRIBUTED_MDCELL_READER_H +#define DISTRIBUTED_MDCELL_READER_H + +#include + +class MdCell; + +class DistributedMdCellReader +{ +public: + static MdCell read_lj_stru(const std::string& stru_file, + double cutoff_bohr, + double skin_bohr); +}; + +#endif diff --git a/source/source_cell/md_cell.cpp b/source/source_cell/md_cell.cpp new file mode 100644 index 00000000000..6ae4fd138c8 --- /dev/null +++ b/source/source_cell/md_cell.cpp @@ -0,0 +1,579 @@ +#include "source_cell/md_cell.h" + +#include "source_cell/unitcell.h" +#include "source_io/module_parameter/parameter.h" + +#include +#include +#include + +double MdCell::wrap_fractional_(double value) +{ + value -= std::floor(value); + if (value >= 1.0 - 1.0e-12 || value < 1.0e-12) + { + return 0.0; + } + return value; +} + +double MdCell::infer_cutoff_from_parameter_(const Parameter& param) +{ + double cutoff = 0.0; + const std::vector& lj_rcut = param.inp.mdp.lj_rcut; + for (std::size_t i = 0; i < lj_rcut.size(); ++i) + { + cutoff = std::max(cutoff, lj_rcut[i] * ModuleBase::ANGSTROM_AU); + } + return cutoff; +} + +void MdCell::clear_forces_(std::vector& atoms) +{ + for (std::size_t i = 0; i < atoms.size(); ++i) + { + atoms[i].force.set(0.0, 0.0, 0.0); + } +} + +void MdCell::sync_backing_unitcell_geometry_() +{ + if (backing_unitcell_ == nullptr) + { + return; + } + + backing_unitcell_->latvec = latvec_; + backing_unitcell_->omega = omega_; + backing_unitcell_->GT = gt_; + backing_unitcell_->G = gt_.Transpose(); + backing_unitcell_->GGT = backing_unitcell_->G * backing_unitcell_->GT; + backing_unitcell_->invGGT = backing_unitcell_->GGT.Inverse(); + backing_unitcell_->lat0_angstrom = lat0_ * ModuleBase::BOHR_TO_A; + backing_unitcell_->tpiba = ModuleBase::TWO_PI / lat0_; + backing_unitcell_->tpiba2 = backing_unitcell_->tpiba * backing_unitcell_->tpiba; + backing_unitcell_->a1.set(latvec_.e11, latvec_.e12, latvec_.e13); + backing_unitcell_->a2.set(latvec_.e21, latvec_.e22, latvec_.e23); + backing_unitcell_->a3.set(latvec_.e31, latvec_.e32, latvec_.e33); + backing_unitcell_->cell_parameter_updated = true; +} + +void MdCell::sync_backing_unitcell_owned_atoms_() +{ + if (backing_unitcell_ == nullptr) + { + return; + } + + for (std::size_t i = 0; i < owned_atoms_.size(); ++i) + { + const LocalAtom& atom = owned_atoms_[i]; + backing_unitcell_->atoms[atom.type].tau[atom.type_index] = atom.cart; + backing_unitcell_->atoms[atom.type].taud[atom.type_index] = atom.frac; + backing_unitcell_->atoms[atom.type].vel[atom.type_index] = atom.vel; + backing_unitcell_->atoms[atom.type].mbl[atom.type_index] = atom.mbl; + } +} + +void MdCell::initialize_from_ucell_serial_(UnitCell& ucell, double cutoff, double skin) +{ + backing_unitcell_ = &ucell; + nat_ = ucell.nat; + lat0_ = ucell.lat0; + omega_ = ucell.omega; + latvec_ = ucell.latvec; + gt_ = ucell.GT; + type_labels_.resize(static_cast(ucell.ntype)); + type_masses_.resize(static_cast(ucell.ntype)); + stru_metadata_.species.resize(static_cast(ucell.ntype)); + for (int it = 0; it < ucell.ntype; ++it) + { + MdStruSpecies& species = stru_metadata_.species[static_cast(it)]; + species.label = ucell.atoms[it].label; + species.mass = ucell.atoms[it].mass; + type_labels_[static_cast(it)] = species.label; + type_masses_[static_cast(it)] = species.mass; + if (static_cast(it) < ucell.pseudo_fn.size()) species.pseudo_file = ucell.pseudo_fn[it]; + if (static_cast(it) < ucell.pseudo_type.size()) species.pseudo_type = ucell.pseudo_type[it]; + if (static_cast(it) < ucell.orbital_fn.size()) species.orbital_file = ucell.orbital_fn[it]; + if (static_cast(it) < ucell.magnet.start_mag.size()) species.start_mag = ucell.magnet.start_mag[it]; + species.atom_count = ucell.atoms[it].na; + } + stru_metadata_.descriptor_file = ucell.descriptor_file; + init_vel_ = ucell.init_vel; + cutoff_ = cutoff; + skin_ = skin; + owned_atoms_.clear(); + ghost_atoms_.clear(); + + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + owned_atoms_.push_back(LocalAtom(ucell.atoms[it].tau[ia], + ucell.atoms[it].taud[ia], + ucell.atoms[it].vel[ia], + ModuleBase::Vector3(0.0, 0.0, 0.0), + ucell.atoms[it].mbl[ia], + ucell.atoms[it].mass / ModuleBase::AU_to_MASS, + it, + ia, + 0, + false)); + } + } + exchange_ghost_atoms(); +} + +#ifdef __MPI +void MdCell::initialize_from_ucell_(UnitCell& ucell, MPI_Comm comm, double cutoff, double skin) +{ + backing_unitcell_ = &ucell; + nat_ = ucell.nat; + lat0_ = ucell.lat0; + omega_ = ucell.omega; + latvec_ = ucell.latvec; + gt_ = ucell.GT; + type_labels_.resize(static_cast(ucell.ntype)); + type_masses_.resize(static_cast(ucell.ntype)); + stru_metadata_.species.resize(static_cast(ucell.ntype)); + for (int it = 0; it < ucell.ntype; ++it) + { + MdStruSpecies& species = stru_metadata_.species[static_cast(it)]; + species.label = ucell.atoms[it].label; + species.mass = ucell.atoms[it].mass; + type_labels_[static_cast(it)] = species.label; + type_masses_[static_cast(it)] = species.mass; + if (static_cast(it) < ucell.pseudo_fn.size()) species.pseudo_file = ucell.pseudo_fn[it]; + if (static_cast(it) < ucell.pseudo_type.size()) species.pseudo_type = ucell.pseudo_type[it]; + if (static_cast(it) < ucell.orbital_fn.size()) species.orbital_file = ucell.orbital_fn[it]; + if (static_cast(it) < ucell.magnet.start_mag.size()) species.start_mag = ucell.magnet.start_mag[it]; + species.atom_count = ucell.atoms[it].na; + } + stru_metadata_.descriptor_file = ucell.descriptor_file; + init_vel_ = ucell.init_vel; + comm_ = comm; + cutoff_ = cutoff; + skin_ = skin; + + owned_atoms_.clear(); + ghost_atoms_.clear(); + MPI_Comm_rank(comm_, &rank_); + MPI_Comm_size(comm_, &size_); + + decomp_.init(comm_, latvec_, lat0_, cutoff_, skin_); + decomp_.split_owned_atoms_from_ucell(ucell, owned_atoms_); + clear_forces_(owned_atoms_); + exchange_ghost_atoms(); +} + +void MdCell::initialize_from_owned_atoms_(MPI_Comm comm, double cutoff, double skin) +{ + comm_ = comm; + cutoff_ = cutoff; + skin_ = skin; + MPI_Comm_rank(comm_, &rank_); + MPI_Comm_size(comm_, &size_); + decomp_.init(comm_, latvec_, lat0_, cutoff_, skin_); + clear_forces_(owned_atoms_); + exchange_ghost_atoms(); +} + +#endif + +#ifndef __MPI +void MdCell::initialize_from_owned_atoms_(double cutoff, double skin) +{ + cutoff_ = cutoff; + skin_ = skin; + clear_forces_(owned_atoms_); + exchange_ghost_atoms(); +} +#endif + +MdCell::MdCell(UnitCell& ucell, const Parameter& param) +{ + const double cutoff = infer_cutoff_from_parameter_(param); +#ifdef __MPI + initialize_from_ucell_(ucell, MPI_COMM_WORLD, cutoff, 0.0); +#else + initialize_from_ucell_serial_(ucell, cutoff, 0.0); +#endif +} + +MdCell::MdCell(const ModuleBase::Matrix3& latvec, + const ModuleBase::Matrix3& gt, + double lat0, + double omega, + int nat, + const std::vector& owned_atoms, + const std::vector& type_labels, + const std::vector& type_masses, + double cutoff, + double skin) +{ + latvec_ = latvec; + gt_ = gt; + lat0_ = lat0; + omega_ = omega; + nat_ = nat; + owned_atoms_ = owned_atoms; + type_labels_ = type_labels; + type_masses_ = type_masses; + init_vel_ = true; +#ifdef __MPI + initialize_from_owned_atoms_(MPI_COMM_WORLD, cutoff, skin); +#else + initialize_from_owned_atoms_(cutoff, skin); +#endif +} + +#ifdef __MPI +MdCell::MdCell(const ModuleBase::Matrix3& latvec, + const ModuleBase::Matrix3& gt, + double lat0, + double omega, + int nat, + const std::vector& owned_atoms, + const std::vector& type_labels, + const std::vector& type_masses, + MPI_Comm comm, + double cutoff, + double skin) +{ + latvec_ = latvec; + gt_ = gt; + lat0_ = lat0; + omega_ = omega; + nat_ = nat; + owned_atoms_ = owned_atoms; + type_labels_ = type_labels; + type_masses_ = type_masses; + init_vel_ = true; + initialize_from_owned_atoms_(comm, cutoff, skin); +} + +MdCell::MdCell(UnitCell& ucell, MPI_Comm comm, double cutoff, double skin) +{ + initialize_from_ucell_(ucell, comm, cutoff, skin); +} + +int MdCell::mpi_rank() const +{ + return rank_; +} + +int MdCell::mpi_size() const +{ + return size_; +} + +MPI_Comm MdCell::communicator() const +{ + return comm_; +} + +const DomainDecomposition& MdCell::decomposition() const +{ + return decomp_; +} +#endif + +void MdCell::exchange_ghost_atoms() +{ +#ifdef __MPI + decomp_.exchange_ghost_atoms(owned_atoms_, ghost_atoms_); + clear_forces_(ghost_atoms_); + return; +#endif + + ghost_atoms_.clear(); + + if (cutoff_ <= 0.0) + { + return; + } + + const ModuleBase::Vector3 a1(latvec_.e11, latvec_.e12, latvec_.e13); + const ModuleBase::Vector3 a2(latvec_.e21, latvec_.e22, latvec_.e23); + const ModuleBase::Vector3 a3(latvec_.e31, latvec_.e32, latvec_.e33); + const ModuleBase::Vector3 a2xa3(a2.y * a3.z - a2.z * a3.y, + a2.z * a3.x - a2.x * a3.z, + a2.x * a3.y - a2.y * a3.x); + const ModuleBase::Vector3 a3xa1(a3.y * a1.z - a3.z * a1.y, + a3.z * a1.x - a3.x * a1.z, + a3.x * a1.y - a3.y * a1.x); + const ModuleBase::Vector3 a1xa2(a1.y * a2.z - a1.z * a2.y, + a1.z * a2.x - a1.x * a2.z, + a1.x * a2.y - a1.y * a2.x); + const double volume = std::abs(a1.x * a2xa3.x + a1.y * a2xa3.y + a1.z * a2xa3.z); + if (volume <= 0.0) + { + throw std::runtime_error("MdCell requires a nonzero cell volume for periodic ghosts."); + } + + const double search_radius = (cutoff_ + skin_) / lat0_; + const int layers[3] = { + static_cast(std::ceil(a2xa3.norm() * search_radius / volume)), + static_cast(std::ceil(a3xa1.norm() * search_radius / volume)), + static_cast(std::ceil(a1xa2.norm() * search_radius / volume)) + }; + for (int ix = -layers[0]; ix <= layers[0]; ++ix) + { + for (int iy = -layers[1]; iy <= layers[1]; ++iy) + { + for (int iz = -layers[2]; iz <= layers[2]; ++iz) + { + if (ix == 0 && iy == 0 && iz == 0) + { + continue; + } + for (std::size_t iat = 0; iat < owned_atoms_.size(); ++iat) + { + LocalAtom image = owned_atoms_[iat]; + const ModuleBase::Vector3 shifted_frac(image.frac.x + ix, + image.frac.y + iy, + image.frac.z + iz); + image.cart = shifted_frac * latvec_; + image.force.set(0.0, 0.0, 0.0); + image.is_ghost = true; + ghost_atoms_.push_back(image); + } + } + } + } +} + +void MdCell::migrate_owned_atoms() +{ +#ifdef __MPI + decomp_.migrate_owned_atoms(owned_atoms_); + exchange_ghost_atoms(); + return; +#endif + for (std::size_t i = 0; i < owned_atoms_.size(); ++i) + { + LocalAtom& atom = owned_atoms_[i]; + atom.frac = atom.cart * gt_; + atom.frac.x = wrap_fractional_(atom.frac.x); + atom.frac.y = wrap_fractional_(atom.frac.y); + atom.frac.z = wrap_fractional_(atom.frac.z); + atom.cart = atom.frac * latvec_; + } + sync_backing_unitcell_owned_atoms_(); +} + +void MdCell::set_lattice_vectors(const ModuleBase::Matrix3& latvec) +{ + latvec_ = latvec; + gt_ = latvec_.Inverse(); + omega_ = std::abs(latvec_.Det()) * lat0_ * lat0_ * lat0_; +#ifdef __MPI + if (comm_ != MPI_COMM_NULL) + { + decomp_.init(comm_, latvec_, lat0_, cutoff_, skin_); + } +#endif + sync_backing_unitcell_geometry_(); +} + +void MdCell::refresh_cart_from_frac() +{ + for (std::size_t i = 0; i < owned_atoms_.size(); ++i) + { + owned_atoms_[i].frac.x = wrap_fractional_(owned_atoms_[i].frac.x); + owned_atoms_[i].frac.y = wrap_fractional_(owned_atoms_[i].frac.y); + owned_atoms_[i].frac.z = wrap_fractional_(owned_atoms_[i].frac.z); + owned_atoms_[i].cart = owned_atoms_[i].frac * latvec_; + } + sync_backing_unitcell_owned_atoms_(); + exchange_ghost_atoms(); +} + +const std::vector& MdCell::owned_atoms() const +{ + return owned_atoms_; +} + +const std::vector& MdCell::ghost_atoms() const +{ + return ghost_atoms_; +} + +const std::vector& MdCell::type_labels() const +{ + return type_labels_; +} + +const std::vector& MdCell::type_masses() const +{ + return type_masses_; +} + +const MdStruMetadata& MdCell::stru_metadata() const +{ + return stru_metadata_; +} + +void MdCell::set_stru_metadata(const MdStruMetadata& metadata) +{ + stru_metadata_ = metadata; +} + +std::vector& MdCell::mutable_owned_atoms() +{ + return owned_atoms_; +} + +std::vector& MdCell::mutable_ghost_atoms() +{ + return ghost_atoms_; +} + +int MdCell::nlocal() const +{ + return static_cast(owned_atoms_.size()); +} + +int MdCell::nghost() const +{ + return static_cast(ghost_atoms_.size()); +} + +bool MdCell::init_vel() const +{ + return init_vel_; +} + +void MdCell::set_init_vel(bool init_vel) +{ + init_vel_ = init_vel; +} + +double MdCell::cutoff() const +{ + return cutoff_; +} + +double MdCell::skin() const +{ + return skin_; +} + +bool MdCell::has_backing_unitcell() const +{ + return backing_unitcell_ != nullptr; +} + +UnitCell& MdCell::backing_unitcell() +{ + assert(backing_unitcell_ != nullptr); + return *backing_unitcell_; +} + +const UnitCell& MdCell::backing_unitcell() const +{ + assert(backing_unitcell_ != nullptr); + return *backing_unitcell_; +} + +void MdCell::sync_backing_unitcell() +{ + if (backing_unitcell_ == nullptr) + { + return; + } + + sync_backing_unitcell_geometry_(); + +#ifdef __MPI + if (size_ > 1) + { + std::vector type_offset(backing_unitcell_->ntype + 1, 0); + for (int it = 0; it < backing_unitcell_->ntype; ++it) + { + type_offset[it + 1] = type_offset[it] + backing_unitcell_->atoms[it].na; + } + + std::vector cart(3 * nat_, 0.0); + std::vector frac(3 * nat_, 0.0); + std::vector vel(3 * nat_, 0.0); + std::vector mbl(3 * nat_, 0); + std::vector owner(nat_, 0); + + for (std::size_t i = 0; i < owned_atoms_.size(); ++i) + { + const LocalAtom& atom = owned_atoms_[i]; + const int iat = type_offset[atom.type] + atom.type_index; + cart[3 * iat] = atom.cart.x; + cart[3 * iat + 1] = atom.cart.y; + cart[3 * iat + 2] = atom.cart.z; + frac[3 * iat] = atom.frac.x; + frac[3 * iat + 1] = atom.frac.y; + frac[3 * iat + 2] = atom.frac.z; + vel[3 * iat] = atom.vel.x; + vel[3 * iat + 1] = atom.vel.y; + vel[3 * iat + 2] = atom.vel.z; + mbl[3 * iat] = atom.mbl.x; + mbl[3 * iat + 1] = atom.mbl.y; + mbl[3 * iat + 2] = atom.mbl.z; + owner[iat] = 1; + } + + MPI_Allreduce(MPI_IN_PLACE, cart.data(), 3 * nat_, MPI_DOUBLE, MPI_SUM, comm_); + MPI_Allreduce(MPI_IN_PLACE, frac.data(), 3 * nat_, MPI_DOUBLE, MPI_SUM, comm_); + MPI_Allreduce(MPI_IN_PLACE, vel.data(), 3 * nat_, MPI_DOUBLE, MPI_SUM, comm_); + MPI_Allreduce(MPI_IN_PLACE, mbl.data(), 3 * nat_, MPI_INT, MPI_SUM, comm_); + MPI_Allreduce(MPI_IN_PLACE, owner.data(), nat_, MPI_INT, MPI_SUM, comm_); + + for (int it = 0; it < backing_unitcell_->ntype; ++it) + { + for (int ia = 0; ia < backing_unitcell_->atoms[it].na; ++ia) + { + const int iat = type_offset[it] + ia; + if (owner[iat] != 1) + { + throw std::runtime_error("MdCell backing UnitCell atom ownership is invalid."); + } + backing_unitcell_->atoms[it].tau[ia].set(cart[3 * iat], cart[3 * iat + 1], cart[3 * iat + 2]); + backing_unitcell_->atoms[it].taud[ia].set(frac[3 * iat], frac[3 * iat + 1], frac[3 * iat + 2]); + backing_unitcell_->atoms[it].vel[ia].set(vel[3 * iat], vel[3 * iat + 1], vel[3 * iat + 2]); + backing_unitcell_->atoms[it].mbl[ia].set(mbl[3 * iat], mbl[3 * iat + 1], mbl[3 * iat + 2]); + } + } + return; + } +#endif + + sync_backing_unitcell_owned_atoms_(); +} + +BaseCell::Kind MdCell::get_kind() const +{ + return Kind::md_cell; +} + +int MdCell::get_nat() const +{ + return nat_; +} + +double MdCell::get_lat0() const +{ + return lat0_; +} + +double MdCell::get_omega() const +{ + return omega_; +} + +const ModuleBase::Matrix3& MdCell::get_latvec() const +{ + return latvec_; +} + +const ModuleBase::Matrix3& MdCell::get_GT() const +{ + return gt_; +} diff --git a/source/source_cell/md_cell.h b/source/source_cell/md_cell.h new file mode 100644 index 00000000000..27fe453a22d --- /dev/null +++ b/source/source_cell/md_cell.h @@ -0,0 +1,142 @@ +#ifndef MD_CELL_H +#define MD_CELL_H + +#include "source_cell/base_cell.h" +#include "source_cell/module_neighlist/local_atom.h" + +#ifdef __MPI +#include "source_cell/module_neighlist/domain_decomposition.h" +#endif + +#include +#include + +class Parameter; +class UnitCell; + +struct MdStruSpecies +{ + std::string label; + double mass = 0.0; + std::string pseudo_file; + std::string pseudo_type; + std::string orbital_file; + double start_mag = 0.0; + int atom_count = 0; +}; + +struct MdStruMetadata +{ + std::vector species; + std::string descriptor_file; +}; + +class MdCell : public BaseCell +{ +public: + MdCell(UnitCell& ucell, const Parameter& param); + MdCell(const ModuleBase::Matrix3& latvec, + const ModuleBase::Matrix3& gt, + double lat0, + double omega, + int nat, + const std::vector& owned_atoms, + const std::vector& type_labels, + const std::vector& type_masses, + double cutoff, + double skin); + +#ifdef __MPI + MdCell(const ModuleBase::Matrix3& latvec, + const ModuleBase::Matrix3& gt, + double lat0, + double omega, + int nat, + const std::vector& owned_atoms, + const std::vector& type_labels, + const std::vector& type_masses, + MPI_Comm comm, + double cutoff, + double skin); + + MdCell(UnitCell& ucell, MPI_Comm comm, double cutoff, double skin); + + int mpi_rank() const; + int mpi_size() const; + MPI_Comm communicator() const; + + const DomainDecomposition& decomposition() const; +#endif + + void exchange_ghost_atoms(); + void migrate_owned_atoms(); + void set_lattice_vectors(const ModuleBase::Matrix3& latvec); + void refresh_cart_from_frac(); + + const std::vector& owned_atoms() const; + const std::vector& ghost_atoms() const; + const std::vector& type_labels() const; + const std::vector& type_masses() const; + const MdStruMetadata& stru_metadata() const; + void set_stru_metadata(const MdStruMetadata& metadata); + std::vector& mutable_owned_atoms(); + std::vector& mutable_ghost_atoms(); + + int nlocal() const; + int nghost() const; + bool init_vel() const; + void set_init_vel(bool init_vel); + double cutoff() const; + double skin() const; + bool has_backing_unitcell() const; + UnitCell& backing_unitcell(); + const UnitCell& backing_unitcell() const; + void sync_backing_unitcell(); + +private: + Kind get_kind() const override; + int get_nat() const override; + double get_lat0() const override; + double get_omega() const override; + const ModuleBase::Matrix3& get_latvec() const override; + const ModuleBase::Matrix3& get_GT() const override; + + static double infer_cutoff_from_parameter_(const Parameter& param); +#ifdef __MPI + void initialize_from_ucell_(UnitCell& ucell, MPI_Comm comm, double cutoff, double skin); +#endif + void initialize_from_ucell_serial_(UnitCell& ucell, double cutoff, double skin); + void sync_backing_unitcell_geometry_(); + void sync_backing_unitcell_owned_atoms_(); + void clear_forces_(std::vector& atoms); + static double wrap_fractional_(double value); +#ifdef __MPI + void initialize_from_owned_atoms_(MPI_Comm comm, double cutoff, double skin); +#else + void initialize_from_owned_atoms_(double cutoff, double skin); +#endif + + int nat_ = 0; + double lat0_ = 0.0; + double omega_ = 0.0; + ModuleBase::Matrix3 latvec_; + ModuleBase::Matrix3 gt_; + std::vector owned_atoms_; + std::vector ghost_atoms_; + std::vector type_labels_; + std::vector type_masses_; + MdStruMetadata stru_metadata_; + bool init_vel_ = false; + double cutoff_ = 0.0; + double skin_ = 0.0; + UnitCell* backing_unitcell_ = nullptr; + +#ifdef __MPI + MPI_Comm comm_ = MPI_COMM_NULL; + int rank_ = 0; + int size_ = 1; + DomainDecomposition decomp_; +#endif +}; + +#endif diff --git a/source/source_cell/module_neighlist/CMakeLists.txt b/source/source_cell/module_neighlist/CMakeLists.txt index dc3e1e7c500..3e6c282f024 100644 --- a/source/source_cell/module_neighlist/CMakeLists.txt +++ b/source/source_cell/module_neighlist/CMakeLists.txt @@ -3,7 +3,6 @@ set(neighbor_search_sources domain_decomposition.cpp neighbor_search.cpp page_allocator.cpp - unitcell_lite.cpp ) add_library( diff --git a/source/source_cell/module_neighlist/atom_provider.h b/source/source_cell/module_neighlist/atom_provider.h deleted file mode 100644 index 3087148dac5..00000000000 --- a/source/source_cell/module_neighlist/atom_provider.h +++ /dev/null @@ -1,74 +0,0 @@ -#ifndef ATOM_PROVIDER_H -#define ATOM_PROVIDER_H - -#include "source_base/vector3.h" -#include "source_base/matrix3.h" - -/** - * @brief Interface for providing atom and lattice information. - * - * This abstract interface defines the minimum set of methods needed by - * the neighbor search module to access atom positions and lattice parameters. - * Any class implementing this interface can be used with NeighborSearch. - * - * @see UnitCell - * @see UnitCellLite - */ -class AtomProvider -{ -public: - /** - * @brief Default destructor. - */ - virtual ~AtomProvider() = default; - - /** - * @brief Get the lattice constant. - * @return Lattice constant in Bohr. - */ - virtual double get_lat0() const = 0; - - /** - * @brief Get the volume of the unit cell. - * @return Unit cell volume in Bohr^3. - */ - virtual double get_omega() const = 0; - - /** - * @brief Get the lattice vectors. - * @return Const reference to the 3x3 lattice vector matrix. - */ - virtual const ModuleBase::Matrix3& get_latvec() const = 0; - - /** - * @brief Get the total number of atoms. - * @return Total atom count. - */ - virtual int get_natom() const = 0; - - /** - * @brief Get the number of atoms of a specific type. - * @param i Type index. - * @return Number of atoms of type i. - */ - virtual int get_na(int i) const = 0; - - /** - * @brief Get the number of atom types. - * @return Number of atom types. - */ - virtual int get_ntype() const = 0; - - /** - * @brief Get the Cartesian coordinates of a specific atom. - * - * Returns the position of the j-th atom of type i. - * - * @param i Type index. - * @param j Atom index within type i. - * @return Cartesian position vector. - */ - virtual ModuleBase::Vector3 get_tau(int i, int j) const = 0; -}; - -#endif // ATOM_PROVIDER_H \ No newline at end of file diff --git a/source/source_cell/module_neighlist/domain_decomposition.cpp b/source/source_cell/module_neighlist/domain_decomposition.cpp index abbf529a682..443a7d6313b 100644 --- a/source/source_cell/module_neighlist/domain_decomposition.cpp +++ b/source/source_cell/module_neighlist/domain_decomposition.cpp @@ -2,6 +2,8 @@ #ifdef __MPI +#include "source_cell/unitcell.h" + #include #include #include @@ -204,26 +206,33 @@ int DomainDecomposition::owner_rank_from_frac(const ModuleBase::Vector3& return rank_from_coords(owner_coords); } -void DomainDecomposition::split_owned_atoms_from_ucell(const AtomProvider& ucell, +void DomainDecomposition::split_owned_atoms_from_ucell(const UnitCell& ucell, std::vector& owned_atoms) const { owned_atoms.clear(); - owned_atoms.reserve(static_cast(ucell.get_natom() / std::max(1, size_) + 1)); + owned_atoms.reserve(static_cast(ucell.nat / std::max(1, size_) + 1)); - ModuleNeighList::GlobalAtomId global_id = 0; - for (int it = 0; it < ucell.get_ntype(); ++it) + for (int it = 0; it < ucell.ntype; ++it) { - for (int ia = 0; ia < ucell.get_na(it); ++ia) + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { - const ModuleBase::Vector3 original_cart = ucell.get_tau(it, ia); - const ModuleBase::Vector3 frac = wrapped_frac_from_cart(original_cart); - const int owner = owner_rank_from_frac(frac); - if (owner == rank_) - { - const ModuleBase::Vector3 wrapped_cart = frac * latvec_; - owned_atoms.push_back(LocalAtom(wrapped_cart, frac, it, ia, global_id, owner, false)); - } - ++global_id; + const ModuleBase::Vector3 original_cart = ucell.atoms[it].tau[ia]; + const ModuleBase::Vector3 frac = wrapped_frac_from_cart(original_cart); + const int owner = owner_rank_from_frac(frac); + if (owner == rank_) + { + const ModuleBase::Vector3 wrapped_cart = frac * latvec_; + owned_atoms.push_back(LocalAtom(wrapped_cart, + frac, + ucell.atoms[it].vel[ia], + ModuleBase::Vector3(0.0, 0.0, 0.0), + ucell.atoms[it].mbl[ia], + ucell.atoms[it].mass / ModuleBase::AU_to_MASS, + it, + ia, + owner, + false)); + } } } } @@ -314,12 +323,21 @@ DomainDecomposition::PackedAtom DomainDecomposition::pack_atom( packed.frac[0] = atom.frac.x; packed.frac[1] = atom.frac.y; packed.frac[2] = atom.frac.z; + packed.vel[0] = atom.vel.x; + packed.vel[1] = atom.vel.y; + packed.vel[2] = atom.vel.z; + packed.force[0] = atom.force.x; + packed.force[1] = atom.force.y; + packed.force[2] = atom.force.z; + packed.mbl[0] = atom.mbl.x; + packed.mbl[1] = atom.mbl.y; + packed.mbl[2] = atom.mbl.z; + packed.mass = atom.mass; packed.image_shift[0] = image_shift[0]; packed.image_shift[1] = image_shift[1]; packed.image_shift[2] = image_shift[2]; packed.type = atom.type; packed.type_index = atom.type_index; - packed.global_id = atom.global_id; packed.owner_rank = atom.owner_rank; return packed; } @@ -331,15 +349,40 @@ LocalAtom DomainDecomposition::unpack_ghost_atom(const PackedAtom& packed) const packed.frac[1] + packed.image_shift[1], packed.frac[2] + packed.image_shift[2]); const ModuleBase::Vector3 cart = image_frac * latvec_; + const ModuleBase::Vector3 vel(packed.vel[0], packed.vel[1], packed.vel[2]); + const ModuleBase::Vector3 force(packed.force[0], packed.force[1], packed.force[2]); + const ModuleBase::Vector3 mbl(packed.mbl[0], packed.mbl[1], packed.mbl[2]); return LocalAtom(cart, frac, + vel, + force, + mbl, + packed.mass, packed.type, packed.type_index, - packed.global_id, packed.owner_rank, true); } +LocalAtom DomainDecomposition::unpack_owned_atom(const PackedAtom& packed) const +{ + const ModuleBase::Vector3 frac(packed.frac[0], packed.frac[1], packed.frac[2]); + const ModuleBase::Vector3 cart = frac * latvec_; + const ModuleBase::Vector3 vel(packed.vel[0], packed.vel[1], packed.vel[2]); + const ModuleBase::Vector3 force(packed.force[0], packed.force[1], packed.force[2]); + const ModuleBase::Vector3 mbl(packed.mbl[0], packed.mbl[1], packed.mbl[2]); + return LocalAtom(cart, + frac, + vel, + force, + mbl, + packed.mass, + packed.type, + packed.type_index, + packed.owner_rank, + false); +} + void DomainDecomposition::exchange_ghost_atoms(const std::vector& owned_atoms, std::vector& ghost_atoms) const { @@ -352,6 +395,7 @@ void DomainDecomposition::exchange_ghost_atoms(const std::vector& own const int span_y = 2 * nlayer[1] + 1; const int span_z = 2 * nlayer[2] + 1; const int lookup_size = (2 * nlayer[0] + 1) * span_y * span_z; + //assert(lookup_size==slots.size()); std::vector slot_lookup(static_cast(lookup_size), -1); for (std::size_t islot = 0; islot < slots.size(); ++islot) { @@ -492,4 +536,69 @@ void DomainDecomposition::exchange_ghost_atoms(const std::vector& own } } +void DomainDecomposition::migrate_owned_atoms(std::vector& owned_atoms) const +{ + std::vector > send_atoms(static_cast(size_)); + for (std::size_t i = 0; i < owned_atoms.size(); ++i) + { + LocalAtom atom = owned_atoms[i]; + atom.frac = wrapped_frac_from_cart(atom.cart); + atom.cart = atom.frac * latvec_; + atom.owner_rank = owner_rank_from_frac(atom.frac); + atom.is_ghost = false; + const std::array no_shift = {{0, 0, 0}}; + send_atoms[static_cast(atom.owner_rank)].push_back(pack_atom(atom, no_shift)); + } + + std::vector send_counts(static_cast(size_), 0); + std::vector recv_counts(static_cast(size_), 0); + for (int irank = 0; irank < size_; ++irank) + { + send_counts[static_cast(irank)] + = static_cast(send_atoms[static_cast(irank)].size() * sizeof(PackedAtom)); + } + MPI_Alltoall(&send_counts[0], 1, MPI_INT, &recv_counts[0], 1, MPI_INT, comm_); + + std::vector send_displs(static_cast(size_), 0); + std::vector recv_displs(static_cast(size_), 0); + int total_send_bytes = 0; + int total_recv_bytes = 0; + for (int irank = 0; irank < size_; ++irank) + { + send_displs[static_cast(irank)] = total_send_bytes; + recv_displs[static_cast(irank)] = total_recv_bytes; + total_send_bytes += send_counts[static_cast(irank)]; + total_recv_bytes += recv_counts[static_cast(irank)]; + } + + std::vector send_buffer(static_cast(total_send_bytes / static_cast(sizeof(PackedAtom)))); + int send_index = 0; + for (int irank = 0; irank < size_; ++irank) + { + const std::vector& atoms = send_atoms[static_cast(irank)]; + for (std::size_t i = 0; i < atoms.size(); ++i) + { + send_buffer[static_cast(send_index++)] = atoms[i]; + } + } + + std::vector recv_buffer(static_cast(total_recv_bytes / static_cast(sizeof(PackedAtom)))); + MPI_Alltoallv(total_send_bytes > 0 ? reinterpret_cast(&send_buffer[0]) : 0, + &send_counts[0], + &send_displs[0], + MPI_BYTE, + total_recv_bytes > 0 ? reinterpret_cast(&recv_buffer[0]) : 0, + &recv_counts[0], + &recv_displs[0], + MPI_BYTE, + comm_); + + owned_atoms.clear(); + owned_atoms.reserve(recv_buffer.size()); + for (std::size_t i = 0; i < recv_buffer.size(); ++i) + { + owned_atoms.push_back(unpack_owned_atom(recv_buffer[i])); + } +} + #endif // __MPI diff --git a/source/source_cell/module_neighlist/domain_decomposition.h b/source/source_cell/module_neighlist/domain_decomposition.h index 9b74729abd9..3b00c3d56fe 100644 --- a/source/source_cell/module_neighlist/domain_decomposition.h +++ b/source/source_cell/module_neighlist/domain_decomposition.h @@ -3,7 +3,8 @@ #ifdef __MPI -#include "source_cell/module_neighlist/atom_provider.h" +#include "source_base/matrix3.h" +#include "source_base/vector3.h" #include "source_cell/module_neighlist/local_atom.h" #include @@ -11,6 +12,8 @@ #include +class UnitCell; + /** * @brief MPI domain decomposition for distributed neighbor-search input. * @@ -32,11 +35,12 @@ class DomainDecomposition int owner_rank_from_frac(const ModuleBase::Vector3& frac) const; - void split_owned_atoms_from_ucell(const AtomProvider& ucell, + void split_owned_atoms_from_ucell(const UnitCell& ucell, std::vector& owned_atoms) const; void exchange_ghost_atoms(const std::vector& owned_atoms, std::vector& ghost_atoms) const; + void migrate_owned_atoms(std::vector& owned_atoms) const; const std::array& dims() const; const std::array& coords() const; @@ -47,10 +51,13 @@ class DomainDecomposition struct PackedAtom { double frac[3]; + double vel[3]; + double force[3]; + int mbl[3]; + double mass; int image_shift[3]; int type; int type_index; - ModuleNeighList::GlobalAtomId global_id; int owner_rank; }; @@ -99,6 +106,7 @@ class DomainDecomposition void build_ghost_exchange_slots(std::vector& slots) const; PackedAtom pack_atom(const LocalAtom& atom, const std::array& image_shift) const; LocalAtom unpack_ghost_atom(const PackedAtom& packed) const; + LocalAtom unpack_owned_atom(const PackedAtom& packed) const; }; #endif // __MPI diff --git a/source/source_cell/module_neighlist/local_atom.h b/source/source_cell/module_neighlist/local_atom.h index f48a8da8f75..5969a8508e1 100644 --- a/source/source_cell/module_neighlist/local_atom.h +++ b/source/source_cell/module_neighlist/local_atom.h @@ -16,18 +16,24 @@ struct LocalAtom { ModuleBase::Vector3 cart; ModuleBase::Vector3 frac; + ModuleBase::Vector3 vel; + ModuleBase::Vector3 force; + ModuleBase::Vector3 mbl; + double mass; int type; int type_index; - ModuleNeighList::GlobalAtomId global_id; int owner_rank; bool is_ghost; LocalAtom() : cart(0.0, 0.0, 0.0), frac(0.0, 0.0, 0.0), + vel(0.0, 0.0, 0.0), + force(0.0, 0.0, 0.0), + mbl(1, 1, 1), + mass(1.0), type(0), type_index(0), - global_id(-1), owner_rank(0), is_ghost(false) { @@ -35,16 +41,22 @@ struct LocalAtom LocalAtom(const ModuleBase::Vector3& cart_in, const ModuleBase::Vector3& frac_in, + const ModuleBase::Vector3& vel_in, + const ModuleBase::Vector3& force_in, + const ModuleBase::Vector3& mbl_in, + double mass_in, int type_in, int type_index_in, - ModuleNeighList::GlobalAtomId global_id_in, int owner_rank_in, bool is_ghost_in) : cart(cart_in), frac(frac_in), + vel(vel_in), + force(force_in), + mbl(mbl_in), + mass(mass_in), type(type_in), type_index(type_index_in), - global_id(global_id_in), owner_rank(owner_rank_in), is_ghost(is_ghost_in) { diff --git a/source/source_cell/module_neighlist/neighbor_atom.h b/source/source_cell/module_neighlist/neighbor_atom.h index 3f62d30571a..5e805d0e21f 100644 --- a/source/source_cell/module_neighlist/neighbor_atom.h +++ b/source/source_cell/module_neighlist/neighbor_atom.h @@ -33,9 +33,6 @@ class NeighborAtom /// Rank-local atom ID used by the neighbor list. ModuleNeighList::LocalAtomIndex atom_id; - /// Global atom ID in the primary cell. Rank-local images share this ID. - ModuleNeighList::GlobalAtomId global_id; - /// MPI rank that owns the primary atom. int owner_rank; @@ -56,8 +53,7 @@ class NeighborAtom int index, ModuleNeighList::LocalAtomIndex id) : position_x(x), position_y(y), position_z(z), - atom_type(type), atom_index(index), atom_id(id), - global_id(id), owner_rank(0) {} + atom_type(type), atom_index(index), atom_id(id), owner_rank(0) {} NeighborAtom(double x, double y, @@ -65,7 +61,6 @@ class NeighborAtom int type, int index, ModuleNeighList::LocalAtomIndex id, - ModuleNeighList::GlobalAtomId global_id_in, int owner_rank_in) : position_x(x), position_y(y), @@ -73,7 +68,6 @@ class NeighborAtom atom_type(type), atom_index(index), atom_id(id), - global_id(global_id_in), owner_rank(owner_rank_in) { } diff --git a/source/source_cell/module_neighlist/neighbor_search.cpp b/source/source_cell/module_neighlist/neighbor_search.cpp index 74e21cfac69..cd5c5125ab4 100644 --- a/source/source_cell/module_neighlist/neighbor_search.cpp +++ b/source/source_cell/module_neighlist/neighbor_search.cpp @@ -1,4 +1,7 @@ #include "source_cell/module_neighlist/neighbor_search.h" +#include "source_cell/md_cell.h" +#include "source_cell/unitcell.h" + #include #include #include @@ -45,7 +48,6 @@ void NeighborSearch::init_distributed(const std::vector& owned_atoms, ghost_atoms_.clear(); all_atoms_.clear(); bin_manager_.clear(); - search_radius_ = sr / lat0; const std::size_t total_atoms = ModuleNeighList::checked_size_sum(owned_atoms.size(), @@ -55,52 +57,40 @@ void NeighborSearch::init_distributed(const std::vector& owned_atoms, { throw std::overflow_error("NeighborSearch distributed atom count exceeds local atom index range."); } - all_atoms_.reserve(total_atoms); inside_atoms_.reserve(owned_atoms.size()); ghost_atoms_.reserve(ghost_atoms.size()); - - for (size_t iat = 0; iat < owned_atoms.size(); ++iat) + for (std::size_t iat = 0; iat < owned_atoms.size(); ++iat) { const LocalAtom& local = owned_atoms[iat]; - NeighborAtom atom(local.cart.x, - local.cart.y, - local.cart.z, - local.type, - local.type_index, - ModuleNeighList::checked_local_atom_index(all_atoms_.size(), - "NeighborSearch owned atom id"), - local.global_id, + NeighborAtom atom(local.cart.x, local.cart.y, local.cart.z, local.type, local.type_index, + ModuleNeighList::checked_local_atom_index(all_atoms_.size(), "NeighborSearch owned atom id"), local.owner_rank); all_atoms_.push_back(atom); inside_atoms_.push_back(atom); } - - for (size_t iat = 0; iat < ghost_atoms.size(); ++iat) + for (std::size_t iat = 0; iat < ghost_atoms.size(); ++iat) { const LocalAtom& local = ghost_atoms[iat]; - NeighborAtom atom(local.cart.x, - local.cart.y, - local.cart.z, - local.type, - local.type_index, - ModuleNeighList::checked_local_atom_index(all_atoms_.size(), - "NeighborSearch ghost atom id"), - local.global_id, + NeighborAtom atom(local.cart.x, local.cart.y, local.cart.z, local.type, local.type_index, + ModuleNeighList::checked_local_atom_index(all_atoms_.size(), "NeighborSearch ghost atom id"), local.owner_rank); all_atoms_.push_back(atom); ghost_atoms_.push_back(atom); } + neighbor_list_.initialize(inside_atoms_.size(), + ModuleNeighList::checked_size_product(all_atoms_.size(), neighbor_reserve_factor, + "NeighborSearch page size")); +} - const std::size_t page_size = ModuleNeighList::checked_size_product(all_atoms_.size(), - neighbor_reserve_factor, - "NeighborSearch page size"); - neighbor_list_.initialize(inside_atoms_.size(), page_size); +void NeighborSearch::init_from_mdcell_(const MdCell& cell, double sr) +{ + init_distributed(cell.owned_atoms(), cell.ghost_atoms(), sr, cell.lat0()); } -void NeighborSearch::init(const AtomProvider& ucell, double sr) +void NeighborSearch::init_from_unitcell_(const UnitCell& ucell, double sr) { - search_radius_ = sr / ucell.get_lat0(); + search_radius_ = sr / ucell.lat0; // clear possible residual data from previous runs inside_atoms_.clear(); @@ -108,17 +98,17 @@ void NeighborSearch::init(const AtomProvider& ucell, double sr) all_atoms_.clear(); bin_manager_.clear(); - for (int i = 0; i < ucell.get_ntype(); i++) + for (int i = 0; i < ucell.ntype; i++) { - for (int j = 0; j < ucell.get_na(i); j++) + for (int j = 0; j < ucell.atoms[i].na; j++) { const ModuleNeighList::LocalAtomIndex atom_count = ModuleNeighList::checked_local_atom_index(all_atoms_.size(), "NeighborSearch atom id"); NeighborAtom atom( - ucell.get_tau(i,j).x, - ucell.get_tau(i,j).y, - ucell.get_tau(i,j).z, + ucell.atoms[i].tau[j].x, + ucell.atoms[i].tau[j].y, + ucell.atoms[i].tau[j].z, i, j, atom_count @@ -144,6 +134,20 @@ void NeighborSearch::init(const AtomProvider& ucell, double sr) neighbor_list_.initialize(inside_atoms_.size(), page_size); } +void NeighborSearch::init(BaseCell& cell, double sr) +{ + if (cell.kind() == BaseCell::Kind::md_cell) + { + MdCell& md_cell = static_cast(cell); + init_from_mdcell_(md_cell, sr); + return; + } + + assert(cell.kind() == BaseCell::Kind::unit_cell); + UnitCell& ucell = static_cast(cell); + init_from_unitcell_(ucell, sr); +} + void NeighborSearch::build_neighbors() { bin_manager_.init_bins(search_radius_, all_atoms_); @@ -163,11 +167,11 @@ double NeighborSearch::cross_product_norm(double a1, double a2, double a3, return sqrt(c1 * c1 + c2 * c2 + c3 * c3); } -void NeighborSearch::check_expand_condition(const AtomProvider& ucell, int& glayerX_minus, int& glayerX, int& glayerY_minus, int& glayerY, int& glayerZ_minus, int& glayerZ) +void NeighborSearch::check_expand_condition(const UnitCell& ucell, int& glayerX_minus, int& glayerX, int& glayerY_minus, int& glayerY, int& glayerZ_minus, int& glayerZ) { - const auto& lat = ucell.get_latvec(); - const double omega = ucell.get_omega(); - const double lat0 = ucell.get_lat0(); + const auto& lat = ucell.latvec; + const double omega = ucell.omega; + const double lat0 = ucell.lat0; const double lat0_cubed = lat0 * lat0 * lat0; double a23_norm = cross_product_norm(lat.e21, lat.e22, lat.e23, lat.e31, lat.e32, lat.e33); @@ -187,11 +191,11 @@ void NeighborSearch::check_expand_condition(const AtomProvider& ucell, int& glay glayerZ_minus = extend_d33; } -void NeighborSearch::set_member_variables(const AtomProvider& ucell, int glayerX_minus, int glayerX, int glayerY_minus, int glayerY, int glayerZ_minus, int glayerZ) +void NeighborSearch::set_member_variables(const UnitCell& ucell, int glayerX_minus, int glayerX, int glayerY_minus, int glayerY, int glayerZ_minus, int glayerZ) { - ModuleBase::Vector3 vec1(ucell.get_latvec().e11, ucell.get_latvec().e12, ucell.get_latvec().e13); - ModuleBase::Vector3 vec2(ucell.get_latvec().e21, ucell.get_latvec().e22, ucell.get_latvec().e23); - ModuleBase::Vector3 vec3(ucell.get_latvec().e31, ucell.get_latvec().e32, ucell.get_latvec().e33); + ModuleBase::Vector3 vec1(ucell.latvec.e11, ucell.latvec.e12, ucell.latvec.e13); + ModuleBase::Vector3 vec2(ucell.latvec.e21, ucell.latvec.e22, ucell.latvec.e23); + ModuleBase::Vector3 vec3(ucell.latvec.e31, ucell.latvec.e32, ucell.latvec.e33); for (int ix = -glayerX_minus; ix < glayerX; ix++) { @@ -203,13 +207,13 @@ void NeighborSearch::set_member_variables(const AtomProvider& ucell, int glayerX { continue; } - for (int i = 0; i < ucell.get_ntype(); i++) + for (int i = 0; i < ucell.ntype; i++) { - for (int j = 0; j < ucell.get_na(i); j++) + for (int j = 0; j < ucell.atoms[i].na; j++) { - double atom_x = ucell.get_tau(i,j).x + vec1[0] * ix + vec2[0] * iy + vec3[0] * iz; - double atom_y = ucell.get_tau(i,j).y + vec1[1] * ix + vec2[1] * iy + vec3[1] * iz; - double atom_z = ucell.get_tau(i,j).z + vec1[2] * ix + vec2[2] * iy + vec3[2] * iz; + double atom_x = ucell.atoms[i].tau[j].x + vec1[0] * ix + vec2[0] * iy + vec3[0] * iz; + double atom_y = ucell.atoms[i].tau[j].y + vec1[1] * ix + vec2[1] * iy + vec3[1] * iz; + double atom_z = ucell.atoms[i].tau[j].z + vec1[2] * ix + vec2[2] * iy + vec3[2] * iz; const ModuleNeighList::LocalAtomIndex atom_count = ModuleNeighList::checked_local_atom_index(all_atoms_.size(), diff --git a/source/source_cell/module_neighlist/neighbor_search.h b/source/source_cell/module_neighlist/neighbor_search.h index b95ace5bc6d..fae73246d83 100644 --- a/source/source_cell/module_neighlist/neighbor_search.h +++ b/source/source_cell/module_neighlist/neighbor_search.h @@ -4,8 +4,11 @@ #include "source_cell/module_neighlist/neighbor_atom.h" #include "source_cell/module_neighlist/bin_manager.h" #include "source_cell/module_neighlist/neighbor_list.h" -#include "source_cell/module_neighlist/atom_provider.h" #include "source_cell/module_neighlist/local_atom.h" +#include "source_cell/base_cell.h" + +class MdCell; +class UnitCell; /** * @brief Neighbor search algorithm for building atom neighbor lists. @@ -43,19 +46,8 @@ class NeighborSearch * @param ucell Unit cell providing atom positions and lattice info. * @param sr Search radius (cutoff distance) in Bohr. */ - void init(const AtomProvider& ucell, double sr); + void init(BaseCell& cell, double sr); - /** - * @brief Initialize from rank-local owned atoms and exchanged ghost atoms. - * - * This distributed entry point does not inspect a global UnitCell. The - * caller is responsible for domain ownership and ghost exchange. - * - * @param owned_atoms Atoms owned by this rank and used as list centers. - * @param ghost_atoms Cutoff halo atoms received from neighboring ranks. - * @param sr Search radius (cutoff distance) in Bohr. - * @param lat0 Lattice constant in Bohr. - */ void init_distributed(const std::vector& owned_atoms, const std::vector& ghost_atoms, double sr, @@ -121,7 +113,11 @@ class NeighborSearch * * @param ucell Unit cell providing lattice vectors. */ - void check_expand_condition(const AtomProvider& ucell, int& glayerX_minus, int& glayerX, int& glayerY_minus, int& glayerY, int& glayerZ_minus, int& glayerZ); + void init_from_unitcell_(const UnitCell& ucell, double sr); + + void init_from_mdcell_(const MdCell& cell, double sr); + + void check_expand_condition(const UnitCell& ucell, int& glayerX_minus, int& glayerX, int& glayerY_minus, int& glayerY, int& glayerZ_minus, int& glayerZ); /** * @brief Set member variables by generating periodic images. @@ -131,7 +127,7 @@ class NeighborSearch * * @param ucell Unit cell providing atom positions. */ - void set_member_variables(const AtomProvider& ucell, int glayerX_minus, int glayerX, int glayerY_minus, int glayerY, int glayerZ_minus, int glayerZ); + void set_member_variables(const UnitCell& ucell, int glayerX_minus, int glayerX, int glayerY_minus, int glayerY, int glayerZ_minus, int glayerZ); // ========== Data members ========== diff --git a/source/source_cell/module_neighlist/neighbor_types.h b/source/source_cell/module_neighlist/neighbor_types.h index a3a95aeb317..739f72445e5 100644 --- a/source/source_cell/module_neighlist/neighbor_types.h +++ b/source/source_cell/module_neighlist/neighbor_types.h @@ -10,7 +10,6 @@ namespace ModuleNeighList { -using GlobalAtomId = std::int64_t; using LocalAtomIndex = std::int32_t; using NeighborCount = std::int32_t; diff --git a/source/source_cell/module_neighlist/test/CMakeLists.txt b/source/source_cell/module_neighlist/test/CMakeLists.txt index 281b4fc5fcb..ae0b535c9f4 100644 --- a/source/source_cell/module_neighlist/test/CMakeLists.txt +++ b/source/source_cell/module_neighlist/test/CMakeLists.txt @@ -12,9 +12,10 @@ AddTest( SOURCES neighbor_search_test.cpp ../neighbor_search.cpp + ../../md_cell.cpp + ../domain_decomposition.cpp ../bin_manager.cpp ../page_allocator.cpp - ../unitcell_lite.cpp ) AddTest( @@ -35,31 +36,43 @@ AddTest( ) if(ENABLE_MPI) - add_executable(MODULE_CELL_NEIGHBOR_neighbor_search_mpi_benchmark - neighbor_search_mpi_benchmark.cpp + add_executable(MODULE_CELL_NEIGHBOR_mdcell_migrate_mpi + md_cell_migrate_mpi_test.cpp + ../../md_cell.cpp ../domain_decomposition.cpp - ../neighbor_search.cpp - ../bin_manager.cpp - ../page_allocator.cpp - ../unitcell_lite.cpp + ) + target_link_libraries(MODULE_CELL_NEIGHBOR_mdcell_migrate_mpi + PRIVATE + parameter base device MPI::MPI_CXX GTest::gtest_main GTest::gmock_main abacus::linalg_libs + ) + install(TARGETS MODULE_CELL_NEIGHBOR_mdcell_migrate_mpi DESTINATION ${CMAKE_BINARY_DIR}/tests) + add_test(NAME MODULE_CELL_NEIGHBOR_mdcell_migrate_mpi + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 + $ + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + ) + + add_executable(MODULE_CELL_NEIGHBOR_distributed_mdcell_reader + distributed_mdcell_reader_test.cpp + ../../distributed_mdcell_reader.cpp + ../../md_cell.cpp + ../domain_decomposition.cpp + ../../../source_base/global_variable.cpp ../../../source_base/matrix.cpp ../../../source_base/matrix3.cpp ../../../source_base/tool_quit.cpp ) - target_include_directories(MODULE_CELL_NEIGHBOR_neighbor_search_mpi_benchmark PRIVATE ${ABACUS_SOURCE_DIR}) - target_compile_definitions(MODULE_CELL_NEIGHBOR_neighbor_search_mpi_benchmark PRIVATE __NORMAL) - target_link_libraries(MODULE_CELL_NEIGHBOR_neighbor_search_mpi_benchmark + target_include_directories(MODULE_CELL_NEIGHBOR_distributed_mdcell_reader PRIVATE ${ABACUS_SOURCE_DIR}) + target_compile_definitions(MODULE_CELL_NEIGHBOR_distributed_mdcell_reader PRIVATE __NORMAL) + target_link_libraries(MODULE_CELL_NEIGHBOR_distributed_mdcell_reader PRIVATE - Threads::Threads MPI::MPI_CXX + Threads::Threads MPI::MPI_CXX GTest::gtest GTest::gmock ) - if(ENABLE_OPENMP) - target_link_libraries(MODULE_CELL_NEIGHBOR_neighbor_search_mpi_benchmark PRIVATE OpenMP::OpenMP_CXX) - endif() - install(TARGETS MODULE_CELL_NEIGHBOR_neighbor_search_mpi_benchmark DESTINATION ${CMAKE_BINARY_DIR}/tests) - add_test(NAME MODULE_CELL_NEIGHBOR_neighbor_search_mpi_benchmark_np4 + install(TARGETS MODULE_CELL_NEIGHBOR_distributed_mdcell_reader DESTINATION ${CMAKE_BINARY_DIR}/tests) + add_test(NAME MODULE_CELL_NEIGHBOR_distributed_mdcell_reader_np4 COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 4 - $ - 12 12 12 2 1.75 1.0 0.2 1 + $ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} ) + endif() diff --git a/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp b/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp new file mode 100644 index 00000000000..7691e9d19ef --- /dev/null +++ b/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp @@ -0,0 +1,134 @@ +#include + +#include "source_cell/distributed_mdcell_reader.h" +#include "source_cell/md_cell.h" +#include "source_base/constants.h" +#include "source_cell/module_neighlist/domain_decomposition.h" + +#include +#include +#include +#include + +namespace +{ +void write_cartesian_stru_case(const std::string& stru_file) +{ + std::ofstream ofs(stru_file.c_str()); + ofs << "ATOMIC_SPECIES\n"; + ofs << "He 4.0026 auto auto\n\n"; + ofs << "LATTICE_CONSTANT\n"; + ofs << "1.0\n\n"; + ofs << "LATTICE_VECTORS\n"; + ofs << "4.0 0.0 0.0\n"; + ofs << "0.0 4.0 0.0\n"; + ofs << "0.0 0.0 4.0\n\n"; + ofs << "ATOMIC_POSITIONS\n"; + ofs << "Cartesian\n\n"; + ofs << "He\n"; + ofs << "0.0\n"; + ofs << "4\n"; + ofs << "0.40 0.40 0.40 m 1 1 1 v 0.01 0.00 0.00\n"; + ofs << "2.40 0.40 0.40 m 1 0 1 v 0.02 0.00 0.00\n"; + ofs << "0.40 2.40 0.40 m 0 1 1 v 0.03 0.00 0.00\n"; + ofs << "2.40 2.40 0.40 m 1 1 0 v 0.04 0.00 0.00\n"; +} + +ModuleBase::Matrix3 make_lattice() +{ + ModuleBase::Matrix3 latvec; + latvec.e11 = 4.0; + latvec.e12 = 0.0; + latvec.e13 = 0.0; + latvec.e21 = 0.0; + latvec.e22 = 4.0; + latvec.e23 = 0.0; + latvec.e31 = 0.0; + latvec.e32 = 0.0; + latvec.e33 = 4.0; + return latvec; +} +} // namespace + +TEST(DistributedMdCellReaderTest, ReadOwnedAtomsFromSTRUWithoutUnitCell) +{ + int rank = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + const std::string stru_file = "distributed_mdcell_reader_cartesian.STRU"; + if (rank == 0) + { + write_cartesian_stru_case(stru_file); + } + MPI_Barrier(MPI_COMM_WORLD); + + MdCell mdcell = DistributedMdCellReader::read_lj_stru(stru_file, + 1.0 * ModuleBase::ANGSTROM_AU, + 0.0); + + EXPECT_EQ(mdcell.type_labels().size(), 1U); + EXPECT_EQ(mdcell.type_labels()[0], "He"); + ASSERT_EQ(mdcell.type_masses().size(), 1U); + EXPECT_DOUBLE_EQ(mdcell.type_masses()[0], 4.0026); + EXPECT_EQ(mdcell.nat(), 4); + + DomainDecomposition decomp; + decomp.init(MPI_COMM_WORLD, make_lattice(), 1.0, 1.0 * ModuleBase::ANGSTROM_AU, 0.0); + + long long local_count = static_cast(mdcell.owned_atoms().size()); + long long global_count = 0; + MPI_Allreduce(&local_count, &global_count, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); + EXPECT_EQ(global_count, 4); + + std::set > local_ids; + for (std::size_t iat = 0; iat < mdcell.owned_atoms().size(); ++iat) + { + const LocalAtom& atom = mdcell.owned_atoms()[iat]; + EXPECT_EQ(decomp.owner_rank_from_frac(atom.frac), rank); + local_ids.insert(std::make_pair(atom.type, atom.type_index)); + EXPECT_GE(atom.type, 0); + EXPECT_DOUBLE_EQ(atom.force.x, 0.0); + EXPECT_DOUBLE_EQ(atom.force.y, 0.0); + EXPECT_DOUBLE_EQ(atom.force.z, 0.0); + } + EXPECT_EQ(local_ids.size(), mdcell.owned_atoms().size()); + + bool saw_v01 = false; + bool saw_v04 = false; + for (std::size_t iat = 0; iat < mdcell.owned_atoms().size(); ++iat) + { + const LocalAtom& atom = mdcell.owned_atoms()[iat]; + if (atom.type == 0 && atom.type_index == 0) + { + saw_v01 = true; + EXPECT_DOUBLE_EQ(atom.vel.x, 0.01); + EXPECT_EQ(atom.mbl.x, 1); + EXPECT_EQ(atom.mbl.y, 1); + EXPECT_EQ(atom.mbl.z, 1); + EXPECT_DOUBLE_EQ(atom.mass, 4.0026 / ModuleBase::AU_to_MASS); + } + if (atom.type == 0 && atom.type_index == 3) + { + saw_v04 = true; + EXPECT_DOUBLE_EQ(atom.vel.x, 0.04); + EXPECT_EQ(atom.mbl.x, 1); + EXPECT_EQ(atom.mbl.y, 1); + EXPECT_EQ(atom.mbl.z, 0); + } + } + + const int saw_flags[2] = {saw_v01 ? 1 : 0, saw_v04 ? 1 : 0}; + int reduced_flags[2] = {0, 0}; + MPI_Allreduce(saw_flags, reduced_flags, 2, MPI_INT, MPI_MAX, MPI_COMM_WORLD); + EXPECT_EQ(reduced_flags[0], 1); + EXPECT_EQ(reduced_flags[1], 1); +} + +int main(int argc, char** argv) +{ + MPI_Init(&argc, &argv); + ::testing::InitGoogleTest(&argc, argv); + const int result = RUN_ALL_TESTS(); + MPI_Finalize(); + return result; +} diff --git a/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp b/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp new file mode 100644 index 00000000000..406b88e02e8 --- /dev/null +++ b/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp @@ -0,0 +1,101 @@ +#include + +#include "source_cell/md_cell.h" + +#include + +#include +#include + +namespace +{ +void ensure_mpi_initialized() +{ + int initialized = 0; + MPI_Initialized(&initialized); + if (!initialized) + { + int provided = 0; + MPI_Init_thread(NULL, NULL, MPI_THREAD_SINGLE, &provided); + } +} + +ModuleBase::Matrix3 make_lattice() +{ + ModuleBase::Matrix3 latvec; + latvec.e11 = 1.0; + latvec.e22 = 1.0; + latvec.e33 = 1.0; + return latvec; +} +} + +TEST(MdCellMigrateMpiTest, AtomCrossingDomainMigratesToNewOwner) +{ + int rank = 0; + int size = 1; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); + ASSERT_GE(size, 2); + + const ModuleBase::Matrix3 latvec = make_lattice(); + std::vector owned_atoms; + if (rank < 2) + { + const ModuleBase::Vector3 frac(rank == 0 ? 0.2 : 0.7, 0.2, 0.2); + owned_atoms.push_back(LocalAtom(frac, + frac, + ModuleBase::Vector3(0.0, 0.0, 0.0), + ModuleBase::Vector3(0.0, 0.0, 0.0), + ModuleBase::Vector3(1, 1, 1), + 1.0, + 0, + rank, + rank, + false)); + } + MdCell mdcell(latvec, + latvec.Inverse(), + 1.0, + 1.0, + 2, + owned_atoms, + std::vector(1, "X"), + std::vector(1, 1.0), + MPI_COMM_WORLD, + 0.1, + 0.0); + + ASSERT_EQ(mdcell.mpi_size(), size); + if (size == 2) + { + if (rank == 0 && mdcell.nlocal() == 1) + { + mdcell.mutable_owned_atoms()[0].cart.x = 0.8; + } + if (rank == 1 && mdcell.nlocal() == 1) + { + mdcell.mutable_owned_atoms()[0].cart.x = 0.3; + } + mdcell.migrate_owned_atoms(); + + long long local_count = mdcell.nlocal(); + long long global_count = 0; + MPI_Allreduce(&local_count, &global_count, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); + EXPECT_EQ(global_count, 2); + + for (int i = 0; i < mdcell.nlocal(); ++i) + { + EXPECT_EQ(mdcell.owned_atoms()[static_cast(i)].owner_rank, rank); + } + } +} + +int main(int argc, char** argv) +{ + MPI_Init(&argc, &argv); + ::testing::InitGoogleTest(&argc, argv); + const int result = RUN_ALL_TESTS(); + MPI_Finalize(); + return result; +} diff --git a/source/source_cell/module_neighlist/test/neighbor_search_mpi_benchmark.cpp b/source/source_cell/module_neighlist/test/neighbor_search_mpi_benchmark.cpp deleted file mode 100644 index 0837d38ebec..00000000000 --- a/source/source_cell/module_neighlist/test/neighbor_search_mpi_benchmark.cpp +++ /dev/null @@ -1,413 +0,0 @@ -#include "source_cell/module_neighlist/neighbor_search.h" -#include "source_cell/module_neighlist/domain_decomposition.h" -#include "source_cell/module_neighlist/neighbor_types.h" -#include "source_cell/module_neighlist/unitcell_lite.h" - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace -{ -int read_int_arg(int argc, char** argv, int index, int fallback) -{ - return argc <= index ? fallback : std::atoi(argv[index]); -} - -double read_double_arg(int argc, char** argv, int index, double fallback) -{ - return argc <= index ? fallback : std::atof(argv[index]); -} - -double cell_volume(const ModuleBase::Matrix3& latvec) -{ - const double cx = latvec.e22 * latvec.e33 - latvec.e23 * latvec.e32; - const double cy = latvec.e23 * latvec.e31 - latvec.e21 * latvec.e33; - const double cz = latvec.e21 * latvec.e32 - latvec.e22 * latvec.e31; - return std::abs(latvec.e11 * cx + latvec.e12 * cy + latvec.e13 * cz); -} - -ModuleBase::Matrix3 make_simple_lattice_latvec(int nx, int ny, int nz, double spacing, double skew) -{ - ModuleBase::Matrix3 latvec; - latvec.e11 = nx * spacing; - latvec.e12 = 0.0; - latvec.e13 = 0.0; - latvec.e21 = skew * ny * spacing; - latvec.e22 = ny * spacing; - latvec.e23 = 0.0; - latvec.e31 = 0.25 * skew * nz * spacing; - latvec.e32 = 0.5 * skew * nz * spacing; - latvec.e33 = nz * spacing; - return latvec; -} - -ModuleBase::Vector3 direct_to_cartesian(const ModuleBase::Matrix3& latvec, - double fx, - double fy, - double fz) -{ - return ModuleBase::Vector3(fx * latvec.e11 + fy * latvec.e21 + fz * latvec.e31, - fx * latvec.e12 + fy * latvec.e22 + fz * latvec.e32, - fx * latvec.e13 + fy * latvec.e23 + fz * latvec.e33); -} - -UnitCellLite make_simple_lattice_ucell(int nx, int ny, int nz, double spacing, double skew) -{ - const ModuleBase::Matrix3 latvec = make_simple_lattice_latvec(nx, ny, nz, spacing, skew); - - std::vector> tau; - tau.reserve(static_cast(nx) * ny * nz); - for (int ix = 0; ix < nx; ++ix) - { - for (int iy = 0; iy < ny; ++iy) - { - for (int iz = 0; iz < nz; ++iz) - { - tau.push_back(direct_to_cartesian(latvec, - static_cast(ix) / nx, - static_cast(iy) / ny, - static_cast(iz) / nz)); - } - } - } - - UnitCellLite ucell; - const double omega = cell_volume(latvec); - ucell.set_lattice(1.0, omega, latvec); - ucell.set_atoms(1, {static_cast(tau.size())}, tau); - return ucell; -} - -long long checked_lattice_atom_count(int nx, int ny, int nz) -{ - const long long lx = nx; - const long long ly = ny; - const long long lz = nz; - if (lx > std::numeric_limits::max() / ly || - lx * ly > std::numeric_limits::max() / lz) - { - throw std::overflow_error("benchmark lattice atom count overflows."); - } - return lx * ly * lz; -} - -long long owner_begin_index(long long n, int coord, int dims) -{ - return (static_cast(coord) * n + dims - 1) / dims; -} - -long long owner_end_index(long long n, int coord, int dims) -{ - return (static_cast(coord + 1) * n + dims - 1) / dims; -} - -void generate_owned_atoms_from_lattice(const DomainDecomposition& decomp, - const ModuleBase::Matrix3& latvec, - int nx, - int ny, - int nz, - std::vector& owned_atoms) -{ - owned_atoms.clear(); - - const auto& coords = decomp.coords(); - const auto& dims = decomp.dims(); - - const long long ix_begin = owner_begin_index(nx, coords[0], dims[0]); - const long long ix_end = owner_end_index(nx, coords[0], dims[0]); - const long long iy_begin = owner_begin_index(ny, coords[1], dims[1]); - const long long iy_end = owner_end_index(ny, coords[1], dims[1]); - const long long iz_begin = owner_begin_index(nz, coords[2], dims[2]); - const long long iz_end = owner_end_index(nz, coords[2], dims[2]); - - const std::size_t local_count - = ModuleNeighList::checked_size_product( - static_cast(ix_end - ix_begin), - ModuleNeighList::checked_size_product(static_cast(iy_end - iy_begin), - static_cast(iz_end - iz_begin), - "benchmark local atom count"), - "benchmark local atom count"); - owned_atoms.reserve(local_count); - - for (long long ix = ix_begin; ix < ix_end; ++ix) - { - for (long long iy = iy_begin; iy < iy_end; ++iy) - { - for (long long iz = iz_begin; iz < iz_end; ++iz) - { - const double fx = static_cast(ix) / nx; - const double fy = static_cast(iy) / ny; - const double fz = static_cast(iz) / nz; - const ModuleBase::Vector3 frac(fx, fy, fz); - const ModuleBase::Vector3 cart = direct_to_cartesian(latvec, fx, fy, fz); - const ModuleNeighList::GlobalAtomId global_id - = static_cast((ix * ny + iy) * nz + iz); - - owned_atoms.push_back(LocalAtom(cart, - frac, - 0, - 0, - global_id, - decomp.rank(), - false)); - } - } - } -} - -long long count_neighbor_pairs(const NeighborList& list) -{ - long long pairs = 0; - for (int local_i = 0; local_i < list.get_nlocal(); ++local_i) - { - pairs += list.get_numneigh(local_i); - } - return pairs; -} - -long long square_sum(long long n) -{ - const __int128 value = static_cast<__int128>(n) * (n - 1) * (2 * n - 1) / 6; - if (value > std::numeric_limits::max()) - { - throw std::overflow_error("benchmark square sum exceeds long long range."); - } - return static_cast(value); -} -} // namespace - -int main(int argc, char** argv) -{ - MPI_Init(&argc, &argv); - - int mpi_rank = 0; - int mpi_size = 1; - MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); - MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); - - if (argc > 1 && std::string(argv[1]) == "--help") - { - if (mpi_rank == 0) - { - std::cout << "Usage: neighbor_search_mpi_benchmark [nx ny nz repeat cutoff spacing skew check_serial]\n" - << "Defaults: nx=16 ny=16 nz=16 repeat=5 cutoff=1.75 spacing=1.0 skew=0.0 check_serial=1\n"; - } - MPI_Finalize(); - return 0; - } - - const int nx = read_int_arg(argc, argv, 1, 16); - const int ny = read_int_arg(argc, argv, 2, 16); - const int nz = read_int_arg(argc, argv, 3, 16); - const int repeat = read_int_arg(argc, argv, 4, 5); - const double cutoff = read_double_arg(argc, argv, 5, 1.75); - const double spacing = read_double_arg(argc, argv, 6, 1.0); - const double skew = read_double_arg(argc, argv, 7, 0.0); - const int check_serial = read_int_arg(argc, argv, 8, 1); - - if (nx <= 0 || ny <= 0 || nz <= 0 || repeat <= 0 || cutoff <= 0.0 || spacing <= 0.0) - { - if (mpi_rank == 0) - { - std::cerr << "All dimensions, repeat, cutoff, and spacing must be positive.\n"; - } - MPI_Finalize(); - return 2; - } - - const ModuleBase::Matrix3 latvec = make_simple_lattice_latvec(nx, ny, nz, spacing, skew); - const double lat0 = 1.0; - const long long nat = checked_lattice_atom_count(nx, ny, nz); - - long long serial_all_atoms = -1; - long long serial_neighbor_pairs = -1; - double serial_init_time = 0.0; - double serial_build_time = 0.0; - if (mpi_rank == 0 && check_serial) - { - UnitCellLite ucell = make_simple_lattice_ucell(nx, ny, nz, spacing, skew); - NeighborSearch serial; - const double t0 = MPI_Wtime(); - serial.init(ucell, cutoff); - const double t1 = MPI_Wtime(); - serial.build_neighbors(); - const double t2 = MPI_Wtime(); - serial_all_atoms = static_cast(serial.get_all_atoms().size()); - serial_neighbor_pairs = count_neighbor_pairs(serial.get_neighbor_list()); - serial_init_time = t1 - t0; - serial_build_time = t2 - t1; - } - MPI_Bcast(&serial_all_atoms, 1, MPI_LONG_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast(&serial_neighbor_pairs, 1, MPI_LONG_LONG, 0, MPI_COMM_WORLD); - - double init_time = 0.0; - double build_time = 0.0; - double total_time = 0.0; - long long last_inside = 0; - long long last_ghost = 0; - long long last_all = 0; - long long last_pairs = 0; - long long inside_index_sum = 0; - long long inside_index_square_sum = 0; - int local_failure = 0; - - for (int i = 0; i < repeat; ++i) - { - MPI_Barrier(MPI_COMM_WORLD); - const double t0 = MPI_Wtime(); - DomainDecomposition decomp; - std::vector owned_atoms; - std::vector ghost_atoms; - NeighborSearch ns; - decomp.init(MPI_COMM_WORLD, latvec, lat0, cutoff, 0.0); - generate_owned_atoms_from_lattice(decomp, latvec, nx, ny, nz, owned_atoms); - decomp.exchange_ghost_atoms(owned_atoms, ghost_atoms); - ns.init_distributed(owned_atoms, ghost_atoms, cutoff, lat0); - const double t1 = MPI_Wtime(); - ns.build_neighbors(); - const double t2 = MPI_Wtime(); - - init_time += t1 - t0; - build_time += t2 - t1; - total_time += t2 - t0; - - if (i == repeat - 1) - { - const auto& inside_atoms = ns.get_inside_atoms(); - const auto& ghost_atoms = ns.get_ghost_atoms(); - const auto& all_atoms = ns.get_all_atoms(); - const auto& list = ns.get_neighbor_list(); - - last_inside = static_cast(inside_atoms.size()); - last_ghost = static_cast(ghost_atoms.size()); - last_all = static_cast(all_atoms.size()); - last_pairs = 0; - inside_index_sum = 0; - inside_index_square_sum = 0; - - for (size_t atom_id = 0; atom_id < all_atoms.size(); ++atom_id) - { - if (all_atoms[atom_id].atom_id != - ModuleNeighList::checked_local_atom_index(atom_id, "benchmark atom id")) - { - local_failure = 1; - } - } - - for (const NeighborAtom& atom : inside_atoms) - { - inside_index_sum += atom.global_id; - inside_index_square_sum += static_cast(atom.global_id) * atom.global_id; - } - - for (int local_i = 0; local_i < list.get_nlocal(); ++local_i) - { - last_pairs += list.get_numneigh(local_i); - for (int ad = 0; ad < list.get_numneigh(local_i); ++ad) - { - const int neighbor_id = list.get_firstneigh(local_i)[ad]; - if (neighbor_id < 0 || static_cast(neighbor_id) >= all_atoms.size()) - { - local_failure = 1; - } - } - } - } - } - - long long global_inside = 0; - long long global_ghost = 0; - long long global_all = 0; - long long global_pairs = 0; - long long global_index_sum = 0; - long long global_index_square_sum = 0; - long long min_all = 0; - long long max_all = 0; - long long min_inside = 0; - long long max_inside = 0; - long long min_ghost = 0; - long long max_ghost = 0; - long long min_pairs = 0; - long long max_pairs = 0; - int global_failure = 0; - MPI_Allreduce(&last_inside, &global_inside, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&last_ghost, &global_ghost, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&last_all, &global_all, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&last_pairs, &global_pairs, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&inside_index_sum, &global_index_sum, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&inside_index_square_sum, &global_index_square_sum, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&last_all, &min_all, 1, MPI_LONG_LONG, MPI_MIN, MPI_COMM_WORLD); - MPI_Allreduce(&last_all, &max_all, 1, MPI_LONG_LONG, MPI_MAX, MPI_COMM_WORLD); - MPI_Allreduce(&last_inside, &min_inside, 1, MPI_LONG_LONG, MPI_MIN, MPI_COMM_WORLD); - MPI_Allreduce(&last_inside, &max_inside, 1, MPI_LONG_LONG, MPI_MAX, MPI_COMM_WORLD); - MPI_Allreduce(&last_ghost, &min_ghost, 1, MPI_LONG_LONG, MPI_MIN, MPI_COMM_WORLD); - MPI_Allreduce(&last_ghost, &max_ghost, 1, MPI_LONG_LONG, MPI_MAX, MPI_COMM_WORLD); - MPI_Allreduce(&last_pairs, &min_pairs, 1, MPI_LONG_LONG, MPI_MIN, MPI_COMM_WORLD); - MPI_Allreduce(&last_pairs, &max_pairs, 1, MPI_LONG_LONG, MPI_MAX, MPI_COMM_WORLD); - MPI_Allreduce(&local_failure, &global_failure, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD); - - double max_init_time = 0.0; - double max_build_time = 0.0; - double max_total_time = 0.0; - MPI_Reduce(&init_time, &max_init_time, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); - MPI_Reduce(&build_time, &max_build_time, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); - MPI_Reduce(&total_time, &max_total_time, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); - - const bool ownership_ok = global_inside == nat && - global_index_sum == nat * (nat - 1) / 2 && - global_index_square_sum == square_sum(nat); - const bool neighbor_pairs_ok = !check_serial || global_pairs == serial_neighbor_pairs; - const bool all_ok = ownership_ok && global_failure == 0 && neighbor_pairs_ok; - - if (mpi_rank == 0) - { - std::cout << "NeighborSearch MPI halo benchmark\n" - << "algorithm fractional_halo_bins\n" - << "np " << mpi_size << "\n" - << "atoms " << nat << "\n" - << "grid " << nx << " " << ny << " " << nz << "\n" - << "repeat " << repeat << "\n" - << "cutoff " << cutoff << "\n" - << "spacing " << spacing << "\n" - << "skew " << skew << "\n" - << "check_serial " << check_serial << "\n" - << "serial_all_atoms " << serial_all_atoms << "\n" - << "serial_neighbor_pairs " << serial_neighbor_pairs << "\n" - << "inside_sum " << global_inside << "\n" - << "inside_min " << min_inside << "\n" - << "inside_max " << max_inside << "\n" - << "ghost_sum " << global_ghost << "\n" - << "ghost_min " << min_ghost << "\n" - << "ghost_max " << max_ghost << "\n" - << "all_atoms_sum " << global_all << "\n" - << "all_atoms_min " << min_all << "\n" - << "all_atoms_max " << max_all << "\n" - << "neighbor_pairs_sum " << global_pairs << "\n" - << "neighbor_pairs_min " << min_pairs << "\n" - << "neighbor_pairs_max " << max_pairs << "\n" - << "time_serial_ref_init " << serial_init_time << "\n" - << "time_serial_ref_build " << serial_build_time << "\n" - << "time_serial_ref_total " << serial_init_time + serial_build_time << "\n" - << "time_init_max_total " << max_init_time << "\n" - << "time_build_max_total " << max_build_time << "\n" - << "time_total_max_total " << max_total_time << "\n" - << "time_init_max_avg " << max_init_time / repeat << "\n" - << "time_build_max_avg " << max_build_time / repeat << "\n" - << "time_total_max_avg " << max_total_time / repeat << "\n" - << "ownership_ok " << (ownership_ok ? 1 : 0) << "\n" - << "neighbor_pairs_ok " << (neighbor_pairs_ok ? 1 : 0) << "\n" - << "neighbor_ids_ok " << (global_failure == 0 ? 1 : 0) << "\n"; - } - - MPI_Finalize(); - return all_ok ? 0 : 1; -} diff --git a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp index c689019691c..761c37f46e5 100644 --- a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp +++ b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp @@ -1,25 +1,34 @@ #include -#include "source_cell/module_neighlist/local_atom.h" -#include "source_cell/module_neighlist/neighbor_search.h" -#include "source_cell/module_neighlist/unitcell_lite.h" +#include "source_cell/md_cell.h" +#include "../neighbor_search.h" +#include + +#include #include +#include #include namespace { -UnitCellLite make_test_ucell(double lat0, - double omega, - const ModuleBase::Matrix3& latvec, - int ntype, - const std::vector& na, - const std::vector>& tau) +void ensure_mpi_initialized() { - UnitCellLite ucell; - ucell.set_lattice(lat0, omega, latvec); - ucell.set_atoms(ntype, na, tau); - return ucell; + int initialized = 0; + MPI_Initialized(&initialized); + if (!initialized) + { + int provided = 0; + MPI_Init_thread(NULL, NULL, MPI_THREAD_SINGLE, &provided); + } +} + +double cell_volume(const ModuleBase::Matrix3& latvec) +{ + const double cx = latvec.e22 * latvec.e33 - latvec.e23 * latvec.e32; + const double cy = latvec.e23 * latvec.e31 - latvec.e21 * latvec.e33; + const double cz = latvec.e21 * latvec.e32 - latvec.e22 * latvec.e31; + return std::abs(latvec.e11 * cx + latvec.e12 * cy + latvec.e13 * cz); } ModuleBase::Matrix3 identity_lattice() @@ -37,6 +46,41 @@ ModuleBase::Matrix3 identity_lattice() return latvec; } +MdCell make_mdcell(const ModuleBase::Matrix3& latvec, + const std::vector >& positions, + double cutoff) +{ + int rank = 0; + MPI_Comm_rank(MPI_COMM_SELF, &rank); + + const ModuleBase::Matrix3 gt = latvec.Inverse(); + std::vector owned_atoms; + owned_atoms.reserve(positions.size()); + for (std::size_t iat = 0; iat < positions.size(); ++iat) + { + owned_atoms.push_back(LocalAtom(positions[iat], + positions[iat] * gt, + ModuleBase::Vector3(0.0, 0.0, 0.0), + ModuleBase::Vector3(0.0, 0.0, 0.0), + ModuleBase::Vector3(1, 1, 1), + 1.0, + 0, + static_cast(iat), + rank, + false)); + } + return MdCell(latvec, + gt, + 1.0, + cell_volume(latvec), + static_cast(positions.size()), + owned_atoms, + std::vector(1, "X"), + std::vector(1, 1.0), + cutoff, + 0.0); +} + std::size_t count_pairs(const NeighborList& list) { std::size_t pairs = 0; @@ -48,17 +92,16 @@ std::size_t count_pairs(const NeighborList& list) } } // namespace -TEST(NeighborSearchTest, TwoAtomsNeighbor) +TEST(NeighborSearchTest, MdCellTwoAtomsNeighbor) { - UnitCellLite ucell = make_test_ucell(1.0, - 1.0, - identity_lattice(), - 1, - {2}, - {{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}); + ensure_mpi_initialized(); + + const ModuleBase::Matrix3 latvec = identity_lattice(); + const std::vector > positions{{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}; + MdCell mdcell = make_mdcell(latvec, positions, 1.0); NeighborSearch ns; - ns.init(ucell, 1.0); + ns.init(mdcell, 1.0); ns.build_neighbors(); const NeighborList& list = ns.get_neighbor_list(); @@ -67,17 +110,16 @@ TEST(NeighborSearchTest, TwoAtomsNeighbor) EXPECT_EQ(list.get_numneigh(1), 8); } -TEST(NeighborSearchTest, NoNeighbor) +TEST(NeighborSearchTest, MdCellNoNeighbor) { - UnitCellLite ucell = make_test_ucell(1.0, - 1.0, - identity_lattice(), - 1, - {2}, - {{0.0, 0.0, 0.0}, {5.0, 0.0, 0.0}}); + ensure_mpi_initialized(); + + const ModuleBase::Matrix3 latvec = identity_lattice(); + const std::vector > positions{{0.0, 0.0, 0.0}, {0.49, 0.0, 0.0}}; + MdCell mdcell = make_mdcell(latvec, positions, 0.1); NeighborSearch ns; - ns.init(ucell, 0.1); + ns.init(mdcell, 0.1); ns.build_neighbors(); const NeighborList& list = ns.get_neighbor_list(); @@ -86,97 +128,46 @@ TEST(NeighborSearchTest, NoNeighbor) EXPECT_EQ(list.get_numneigh(1), 0); } -TEST(NeighborSearchTest, SerialInitOwnsCentralAtomsAndBuildsImages) +TEST(NeighborSearchTest, MdCellInitBuildsOwnedAndGhostAtoms) { - UnitCellLite ucell = make_test_ucell(1.0, - 1.0, - identity_lattice(), - 1, - {2}, - {{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}); + ensure_mpi_initialized(); + + const ModuleBase::Matrix3 latvec = identity_lattice(); + const std::vector > positions{{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}; + MdCell mdcell = make_mdcell(latvec, positions, 1.0); NeighborSearch ns; - ns.init(ucell, 1.0); + ns.init(mdcell, 1.0); + EXPECT_EQ(mdcell.mpi_size(), 1); EXPECT_EQ(ns.get_inside_atoms().size(), 2U); + EXPECT_GT(ns.get_ghost_atoms().size(), 0U); + EXPECT_GT(ns.get_all_atoms().size(), ns.get_inside_atoms().size()); EXPECT_EQ(ns.get_neighbor_list().get_nlocal(), 2); - EXPECT_EQ(ns.get_all_atoms().size(), 54U); - - const std::vector& all_atoms = ns.get_all_atoms(); - for (std::size_t i = 0; i < all_atoms.size(); ++i) - { - EXPECT_EQ(all_atoms[i].atom_id, - ModuleNeighList::checked_local_atom_index(i, "test atom id")); - } } -TEST(NeighborSearchTest, DistributedInputUsesOwnedCentersAndGhostNeighbors) +TEST(NeighborSearchTest, MdCellNeighborIdsStayLocalToAllAtoms) { - std::vector owned_atoms; - std::vector ghost_atoms; - owned_atoms.push_back(LocalAtom(ModuleBase::Vector3(0.0, 0.0, 0.0), - ModuleBase::Vector3(0.0, 0.0, 0.0), - 0, - 0, - 0, - 0, - false)); - ghost_atoms.push_back(LocalAtom(ModuleBase::Vector3(0.5, 0.0, 0.0), - ModuleBase::Vector3(0.5, 0.0, 0.0), - 0, - 1, - 1, - 1, - true)); - - NeighborSearch ns; - ns.init_distributed(owned_atoms, ghost_atoms, 1.0, 1.0); - ns.build_neighbors(); - - const NeighborList& list = ns.get_neighbor_list(); - ASSERT_EQ(list.get_nlocal(), 1); - ASSERT_EQ(list.get_numneigh(0), 1); - - const int neighbor_id = list.get_firstneigh(0)[0]; - ASSERT_GE(neighbor_id, 0); - ASSERT_LT(static_cast(neighbor_id), ns.get_all_atoms().size()); - EXPECT_EQ(ns.get_all_atoms()[neighbor_id].global_id, 1); - EXPECT_EQ(ns.get_all_atoms()[neighbor_id].owner_rank, 1); -} + ensure_mpi_initialized(); -TEST(NeighborSearchTest, DistributedNeighborIdsStayLocalToAllAtoms) -{ - std::vector owned_atoms; - std::vector ghost_atoms; - owned_atoms.push_back(LocalAtom(ModuleBase::Vector3(0.0, 0.0, 0.0), - ModuleBase::Vector3(0.0, 0.0, 0.0), - 0, - 10, - 0, - 0, - false)); - owned_atoms.push_back(LocalAtom(ModuleBase::Vector3(2.0, 0.0, 0.0), - ModuleBase::Vector3(2.0, 0.0, 0.0), - 0, - 11, - 1, - 0, - false)); - ghost_atoms.push_back(LocalAtom(ModuleBase::Vector3(0.5, 0.0, 0.0), - ModuleBase::Vector3(0.5, 0.0, 0.0), - 0, - 20, - 2, - 1, - true)); + const ModuleBase::Matrix3 latvec = identity_lattice(); + const std::vector > positions{{0.0, 0.0, 0.0}, + {0.5, 0.0, 0.0}, + {0.0, 0.5, 0.0}}; + MdCell mdcell = make_mdcell(latvec, positions, 0.75); NeighborSearch ns; - ns.init_distributed(owned_atoms, ghost_atoms, 0.75, 1.0); + ns.init(mdcell, 0.75); ns.build_neighbors(); const NeighborList& list = ns.get_neighbor_list(); const std::vector& all_atoms = ns.get_all_atoms(); - EXPECT_EQ(count_pairs(list), 1U); + EXPECT_GT(count_pairs(list), 0U); + for (std::size_t atom_id = 0; atom_id < all_atoms.size(); ++atom_id) + { + EXPECT_EQ(all_atoms[atom_id].atom_id, + ModuleNeighList::checked_local_atom_index(atom_id, "test atom id")); + } for (int local_i = 0; local_i < list.get_nlocal(); ++local_i) { for (int ad = 0; ad < list.get_numneigh(local_i); ++ad) @@ -187,3 +178,83 @@ TEST(NeighborSearchTest, DistributedNeighborIdsStayLocalToAllAtoms) } } } + +TEST(NeighborSearchTest, MdCellPreservesMdAtomStateAcrossOwnedAndGhostStorage) +{ + ensure_mpi_initialized(); + + const ModuleBase::Matrix3 latvec = identity_lattice(); + const std::vector > positions{{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}; + MdCell mdcell = make_mdcell(latvec, positions, 1.0); + + ASSERT_EQ(mdcell.nlocal(), 2); + std::vector& owned_atoms = mdcell.mutable_owned_atoms(); + owned_atoms[0].vel.set(1.0, 2.0, 3.0); + owned_atoms[0].mbl.set(1, 0, 1); + owned_atoms[0].mass = 7.5; + owned_atoms[1].vel.set(-1.0, -2.0, -3.0); + owned_atoms[1].mbl.set(0, 1, 1); + owned_atoms[1].mass = 8.5; + + mdcell.exchange_ghost_atoms(); + + ASSERT_GT(mdcell.nghost(), 0); + const std::vector& ghost_atoms = mdcell.ghost_atoms(); + EXPECT_DOUBLE_EQ(ghost_atoms[0].force.x, 0.0); + EXPECT_DOUBLE_EQ(ghost_atoms[0].force.y, 0.0); + EXPECT_DOUBLE_EQ(ghost_atoms[0].force.z, 0.0); + + bool found_first = false; + bool found_second = false; + for (std::size_t i = 0; i < ghost_atoms.size(); ++i) + { + if (ghost_atoms[i].type == owned_atoms[0].type && + ghost_atoms[i].type_index == owned_atoms[0].type_index) + { + found_first = true; + EXPECT_DOUBLE_EQ(ghost_atoms[i].vel.x, 1.0); + EXPECT_DOUBLE_EQ(ghost_atoms[i].vel.y, 2.0); + EXPECT_DOUBLE_EQ(ghost_atoms[i].vel.z, 3.0); + EXPECT_EQ(ghost_atoms[i].mbl.x, 1); + EXPECT_EQ(ghost_atoms[i].mbl.y, 0); + EXPECT_EQ(ghost_atoms[i].mbl.z, 1); + EXPECT_DOUBLE_EQ(ghost_atoms[i].mass, 7.5); + } + if (ghost_atoms[i].type == owned_atoms[1].type && + ghost_atoms[i].type_index == owned_atoms[1].type_index) + { + found_second = true; + EXPECT_DOUBLE_EQ(ghost_atoms[i].vel.x, -1.0); + EXPECT_DOUBLE_EQ(ghost_atoms[i].vel.y, -2.0); + EXPECT_DOUBLE_EQ(ghost_atoms[i].vel.z, -3.0); + EXPECT_EQ(ghost_atoms[i].mbl.x, 0); + EXPECT_EQ(ghost_atoms[i].mbl.y, 1); + EXPECT_EQ(ghost_atoms[i].mbl.z, 1); + EXPECT_DOUBLE_EQ(ghost_atoms[i].mass, 8.5); + } + } + + EXPECT_TRUE(found_first); + EXPECT_TRUE(found_second); +} + +TEST(NeighborSearchTest, MdCellMigrateOwnedAtomsReassignsOwnership) +{ + ensure_mpi_initialized(); + + const ModuleBase::Matrix3 latvec = identity_lattice(); + const std::vector > positions{{0.1, 0.1, 0.1}}; + MdCell mdcell = make_mdcell(latvec, positions, 0.2); + + ASSERT_EQ(mdcell.nlocal(), 1); + mdcell.mutable_owned_atoms()[0].cart.set(1.2, -0.1, 0.1); + mdcell.migrate_owned_atoms(); + + ASSERT_EQ(mdcell.nlocal(), 1); + EXPECT_DOUBLE_EQ(mdcell.owned_atoms()[0].frac.x, 0.2); + EXPECT_DOUBLE_EQ(mdcell.owned_atoms()[0].frac.y, 0.9); + EXPECT_DOUBLE_EQ(mdcell.owned_atoms()[0].frac.z, 0.1); + EXPECT_DOUBLE_EQ(mdcell.owned_atoms()[0].cart.x, 0.2); + EXPECT_DOUBLE_EQ(mdcell.owned_atoms()[0].cart.y, 0.9); + EXPECT_DOUBLE_EQ(mdcell.owned_atoms()[0].cart.z, 0.1); +} diff --git a/source/source_cell/module_neighlist/unitcell_lite.cpp b/source/source_cell/module_neighlist/unitcell_lite.cpp deleted file mode 100644 index 877d214497d..00000000000 --- a/source/source_cell/module_neighlist/unitcell_lite.cpp +++ /dev/null @@ -1,95 +0,0 @@ -#include "unitcell_lite.h" -#include "source_cell/module_neighlist/neighbor_types.h" - -#include - -// === AtomProvider interface implementation === - -double UnitCellLite::get_lat0() const { - return lat0_; -} - -double UnitCellLite::get_omega() const { - return omega_; -} - -const ModuleBase::Matrix3& UnitCellLite::get_latvec() const { - return latvec_; -} - -int UnitCellLite::get_natom() const { - return nat_; -} - -int UnitCellLite::get_na(int i) const { - assert(i >= 0 && i < ntype_); - return na_[i]; -} - -int UnitCellLite::get_ntype() const { - return ntype_; -} - -ModuleBase::Vector3 UnitCellLite::get_tau(int i, int j) const { - assert(i >= 0 && i < ntype_); - assert(j >= 0 && j < na_[i]); - if (i == 0) { - return tau_[j]; - } - return tau_[naa_[i - 1] + j]; -} - -// === Setter methods === - -void UnitCellLite::set_lat0(double lat0) { - lat0_ = lat0; -} - -void UnitCellLite::set_omega(double omega) { - omega_ = omega; -} - -void UnitCellLite::set_latvec(const ModuleBase::Matrix3& latvec) { - latvec_ = latvec; -} - -void UnitCellLite::set_lattice(double lat0, double omega, const ModuleBase::Matrix3& latvec) { - lat0_ = lat0; - omega_ = omega; - latvec_ = latvec; -} - -void UnitCellLite::set_atoms(int ntype, - const std::vector& na, - const std::vector>& tau) { - assert(ntype >= 0); - assert(na.size() == static_cast(ntype)); - - ntype_ = ntype; - na_ = na; - tau_ = tau; - - // compute total number of atoms - std::size_t nat = 0; - for (int i = 0; i < ntype_; ++i) { - assert(na_[i] >= 0); - nat += static_cast(na_[i]); - } - nat_ = ModuleNeighList::checked_int_size(nat, "UnitCellLite atom count"); - assert(tau_.size() == static_cast(nat_)); - - // compute cumulative counts - compute_naa_(); -} - -// === Internal methods === - -void UnitCellLite::compute_naa_() { - naa_.resize(na_.size()); - if (naa_.size() > 0) { - naa_[0] = na_[0]; - } - for (size_t i = 1; i < naa_.size(); ++i) { - naa_[i] = naa_[i - 1] + na_[i]; - } -} diff --git a/source/source_cell/module_neighlist/unitcell_lite.h b/source/source_cell/module_neighlist/unitcell_lite.h deleted file mode 100644 index e951945799f..00000000000 --- a/source/source_cell/module_neighlist/unitcell_lite.h +++ /dev/null @@ -1,169 +0,0 @@ -#ifndef UNITCELL_LITE_H -#define UNITCELL_LITE_H - -#include "source_cell/module_neighlist/atom_provider.h" -#include - -/** - * @brief A lightweight unit cell class for molecular dynamics simulations. - * - * This class provides a minimal set of unit cell information needed for - * large-scale molecular dynamics simulations (e.g., billion-atom simulations). - * It implements the AtomProvider interface and stores only essential data: - * lattice parameters and atomic coordinates. - * - * Compared to the full UnitCell class, UnitCellLite has significantly lower - * memory overhead by omitting electronic structure-related data such as - * pseudopotentials, orbitals, magnetism, and symmetry information. - * - * @see AtomProvider - * @see UnitCell - */ -class UnitCellLite : public AtomProvider -{ -public: - /** - * @brief Default constructor. - * - * Initializes all data members to zero/empty state. - */ - UnitCellLite() = default; - - /** - * @brief Default destructor. - */ - ~UnitCellLite() = default; - - // ========== AtomProvider interface implementation ========== - - /** - * @brief Get the lattice constant in Bohr. - * @return Lattice constant lat0. - */ - double get_lat0() const override; - - /** - * @brief Get the unit cell volume. - * @return Cell volume omega in Bohr^3. - */ - double get_omega() const override; - - /** - * @brief Get the lattice vectors. - * @return Reference to the 3x3 matrix of lattice vectors. - */ - const ModuleBase::Matrix3& get_latvec() const override; - - /** - * @brief Get the total number of atoms. - * @return Total atom count nat. - */ - int get_natom() const override; - - /** - * @brief Get the number of atoms for a given type. - * @param i Atom type index (0-based). - * @return Number of atoms of type i. - * @note Asserts that i is in valid range [0, ntype_). - */ - int get_na(int i) const override; - - /** - * @brief Get the number of atom types. - * @return Number of atom types ntype. - */ - int get_ntype() const override; - - /** - * @brief Get the coordinate of atom (type i, index j). - * @param i Atom type index (0-based). - * @param j Atom index within type i (0-based). - * @return Cartesian coordinate of the atom in Bohr. - * @note Asserts that i and j are in valid ranges. - */ - ModuleBase::Vector3 get_tau(int i, int j) const override; - - // ========== Setter methods ========== - - /** - * @brief Set the lattice constant. - * @param lat0 Lattice constant in Bohr. - */ - void set_lat0(double lat0); - - /** - * @brief Set the unit cell volume. - * @param omega Cell volume in Bohr^3. - */ - void set_omega(double omega); - - /** - * @brief Set the lattice vectors. - * @param latvec 3x3 matrix of lattice vectors. - */ - void set_latvec(const ModuleBase::Matrix3& latvec); - - /** - * @brief Set all lattice parameters together. - * @param lat0 Lattice constant in Bohr. - * @param omega Cell volume in Bohr^3. - * @param latvec 3x3 matrix of lattice vectors. - */ - void set_lattice(double lat0, double omega, const ModuleBase::Matrix3& latvec); - - /** - * @brief Set atom information for all types. - * - * This method sets the number of atom types, the count of atoms per type, - * and all atomic coordinates. It automatically computes the total atom - * count (nat_) and the cumulative atom counts (naa_). - * - * @param ntype Number of atom types. - * @param na Vector of atom counts for each type [ntype]. - * @param tau Vector of all atomic coordinates [nat]. - * - * @note Asserts that na.size() == ntype and tau.size() == sum(na). - */ - void set_atoms(int ntype, - const std::vector& na, - const std::vector>& tau); - -private: - // ========== Data members ========== - - /// Lattice constant in Bohr - double lat0_ = 0.0; - - /// Unit cell volume in Bohr^3 - double omega_ = 0.0; - - /// Total number of atoms - int nat_ = 0; - - /// Number of atom types - int ntype_ = 0; - - /// Lattice vectors (3x3 matrix) - ModuleBase::Matrix3 latvec_; - - /// Number of atoms for each type [ntype] - std::vector na_; - - /// Cumulative sum of na: naa_[i] = na_[0] + na_[1] + ... + na_[i] - std::vector naa_; - - /// Atomic coordinates in Cartesian (Bohr) [nat] - std::vector> tau_; - - // ========== Internal methods ========== - - /** - * @brief Compute cumulative atom counts from na_. - * - * Updates naa_ such that naa_[i] = sum of na_[0] to na_[i]. - * Called internally by set_atoms(). - */ - void compute_naa_(); -}; - -#endif // UNITCELL_LITE_H \ No newline at end of file diff --git a/source/source_cell/unitcell.h b/source/source_cell/unitcell.h index 2be50a98777..bc992f16280 100644 --- a/source/source_cell/unitcell.h +++ b/source/source_cell/unitcell.h @@ -6,49 +6,32 @@ #include "source_cell/sep_cell.h" #include "source_cell/magnetism.h" #include "module_symmetry/symmetry.h" -#include "source_cell/module_neighlist/atom_provider.h" #include "source_cell/base_cell.h" #include "source_cell/nonlocal_info_base.h" /** * @brief Provide the basic information about unitcell. */ -class UnitCell : public AtomProvider, public BaseCell { +class UnitCell : public BaseCell { public: UnitCell(); ~UnitCell(); - /// @name BaseCell / AtomProvider interface overrides - /// @{ - double get_lat0() const override { + double get_lat0() const override + { return lat0; } - double get_omega() const override { + double get_omega() const override + { return omega; } - const ModuleBase::Matrix3& get_latvec() const override { + 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_tau(int i, int j) const override { - return atoms[i].tau[j]; - } - /// @} - /// @brief Initialize basic cell parameters (latname, ntype, lmaxmax, init_vel) /// from INPUT and parse fixed_axes into lat_axis_free flags. void setup_from_input(const std::string& latname_in, diff --git a/source/source_esolver/esolver_lj.cpp b/source/source_esolver/esolver_lj.cpp index 30b8ff78f5a..3b499d1adde 100644 --- a/source/source_esolver/esolver_lj.cpp +++ b/source/source_esolver/esolver_lj.cpp @@ -21,27 +21,6 @@ namespace ModuleESolver { - UnitCellLite ESolver_LJ::change_from_ucell_to_ucell_lite(const UnitCell& ucell) - { - UnitCellLite ucell_lite; - - // Set lattice parameters - ucell_lite.set_lattice(ucell.lat0, ucell.omega, ucell.latvec); - - // Build atom information - std::vector na; - std::vector> tau; - for (int i = 0; i < ucell.ntype; i++) { - na.push_back(ucell.atoms[i].na); - for (int j = 0; j < ucell.atoms[i].na; j++) { - tau.push_back(ucell.atoms[i].tau[j]); - } - } - ucell_lite.set_atoms(ucell.ntype, na, tau); - - return ucell_lite; - } - void ESolver_LJ::before_all_runners(BaseCell& cell, const Input_para& inp) { cell.require_kind(BaseCell::Kind::unit_cell, __FUNCTION__); @@ -65,7 +44,6 @@ void ESolver_LJ::runner(BaseCell& cell, const int istep) static_cast(istep); cell.require_kind(BaseCell::Kind::unit_cell, __FUNCTION__); UnitCell& ucell = static_cast(cell); - UnitCellLite ucell_lite = change_from_ucell_to_ucell_lite(ucell); NeighborSearch neighbor_search; // Important! potential, force, virial must be zero per step @@ -81,12 +59,12 @@ void ESolver_LJ::runner(BaseCell& cell, const int istep) ModuleBase::timer::start("ESolverLJ", "mpi_total"); ModuleBase::timer::start("ESolverLJ", "neigh_init"); DomainDecomposition decomp; - decomp.init(MPI_COMM_WORLD, ucell_lite.get_latvec(), ucell_lite.get_lat0(), search_radius, 0.0); + decomp.init(MPI_COMM_WORLD, ucell.latvec, ucell.lat0, search_radius, 0.0); std::vector owned_atoms; std::vector ghost_atoms; - decomp.split_owned_atoms_from_ucell(ucell_lite, owned_atoms); + decomp.split_owned_atoms_from_ucell(ucell, owned_atoms); decomp.exchange_ghost_atoms(owned_atoms, ghost_atoms); - neighbor_search.init_distributed(owned_atoms, ghost_atoms, search_radius, ucell_lite.get_lat0()); + neighbor_search.init_distributed(owned_atoms, ghost_atoms, search_radius, ucell.lat0); ModuleBase::timer::end("ESolverLJ", "neigh_init"); ModuleBase::timer::start("ESolverLJ", "neigh_bld"); neighbor_search.build_neighbors(); @@ -181,7 +159,7 @@ void ESolver_LJ::runner(BaseCell& cell, const int istep) { ModuleBase::timer::start("ESolverLJ", "serial_tot"); ModuleBase::timer::start("ESolverLJ", "ser_neigh"); - neighbor_search.init(ucell_lite, search_radius); + neighbor_search.init(ucell, search_radius); neighbor_search.build_neighbors(); ModuleBase::timer::end("ESolverLJ", "ser_neigh"); diff --git a/source/source_esolver/esolver_lj.h b/source/source_esolver/esolver_lj.h index 42ed6cfcc71..fd0b390e1ed 100644 --- a/source/source_esolver/esolver_lj.h +++ b/source/source_esolver/esolver_lj.h @@ -2,7 +2,6 @@ #define ESOLVER_LJ_H #include "esolver.h" -#include "source_cell/module_neighlist/unitcell_lite.h" namespace ModuleESolver { @@ -15,8 +14,6 @@ namespace ModuleESolver classname = "ESolver_LJ"; } - UnitCellLite change_from_ucell_to_ucell_lite(const UnitCell& ucell); - void before_all_runners(BaseCell& cell, const Input_para& inp) override; void runner(BaseCell& cell, const int istep) override; diff --git a/source/source_md/test/CMakeLists.txt b/source/source_md/test/CMakeLists.txt index 9f74a519a76..cb72982a3fe 100644 --- a/source/source_md/test/CMakeLists.txt +++ b/source/source_md/test/CMakeLists.txt @@ -4,6 +4,7 @@ abacus_add_local_feature_definitions(__NORMAL) list(APPEND depend_files ../md_func.cpp ../../source_cell/base_cell.cpp + ../../source_cell/md_cell.cpp ../../source_cell/unitcell.cpp ../../source_cell/update_cell.cpp ../../source_cell/bcast_cell.cpp @@ -50,7 +51,6 @@ list(APPEND depend_files ../../source_cell/module_neighlist/neighbor_search.cpp ../../source_cell/module_neighlist/bin_manager.cpp ../../source_cell/module_neighlist/page_allocator.cpp - ../../source_cell/module_neighlist/unitcell_lite.cpp ../../source_base/output.cpp ../../source_io/module_output/output_log.cpp ../../source_io/module_output/print_info.cpp From 84714ba2f91cc4c182623f945414b5778fb46cb6 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Mon, 3 Aug 2026 19:43:41 +0800 Subject: [PATCH 02/10] fix: update NeighborAtom test construction --- source/source_cell/module_neighlist/test/bin_manager_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/source_cell/module_neighlist/test/bin_manager_test.cpp b/source/source_cell/module_neighlist/test/bin_manager_test.cpp index 07e34c488ac..274aefe2089 100644 --- a/source/source_cell/module_neighlist/test/bin_manager_test.cpp +++ b/source/source_cell/module_neighlist/test/bin_manager_test.cpp @@ -153,7 +153,7 @@ TEST(BinManagerUnit, GhostAtomsAreCounted) 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, 3, 1); + ghost.emplace_back(0.4, 0.0, 0.0, 0, 1, 1, 1); BinManager bm; std::vector all_atoms = inside; From dac05c738b427786c64a0b75afbf9902cee0984c Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Mon, 10 Aug 2026 19:47:34 +0800 Subject: [PATCH 03/10] refactor: generalize distributed MDCell reader --- source/source_cell/CMakeLists.txt | 2 - .../source_cell/distributed_mdcell_reader.cpp | 109 +++++++---- .../source_cell/distributed_mdcell_reader.h | 12 +- source/source_cell/md_cell.cpp | 169 ++++++++++-------- source/source_cell/md_cell.h | 17 +- .../module_neighlist/domain_decomposition.cpp | 92 ++++++++++ .../module_neighlist/domain_decomposition.h | 9 + .../module_neighlist/neighbor_search.cpp | 58 +++--- .../module_neighlist/neighbor_search.h | 9 +- .../module_neighlist/test/CMakeLists.txt | 2 +- .../test/distributed_mdcell_reader_test.cpp | 9 +- .../test/md_cell_migrate_mpi_test.cpp | 2 +- .../test/neighbor_search_test.cpp | 17 +- 13 files changed, 346 insertions(+), 161 deletions(-) diff --git a/source/source_cell/CMakeLists.txt b/source/source_cell/CMakeLists.txt index 3aa11e61dfa..5351fa730b9 100644 --- a/source/source_cell/CMakeLists.txt +++ b/source/source_cell/CMakeLists.txt @@ -24,7 +24,6 @@ add_library( klist.cpp parallel_kpoints.cpp cell_index.cpp - cell_tools.cpp check_atomic_stru.cpp update_cell.cpp magnetism.cpp @@ -40,7 +39,6 @@ add_library( cal_nelec_nband.cpp read_pseudo.cpp cal_wfc.cpp - cal_ux.cpp ) if(ENABLE_COVERAGE) diff --git a/source/source_cell/distributed_mdcell_reader.cpp b/source/source_cell/distributed_mdcell_reader.cpp index 314ac05eea7..797f6dbb16e 100644 --- a/source/source_cell/distributed_mdcell_reader.cpp +++ b/source/source_cell/distributed_mdcell_reader.cpp @@ -9,6 +9,8 @@ #endif #include +#include +#include #include #include #include @@ -180,24 +182,41 @@ StruMetadata parse_stru_metadata(std::ifstream& ifs) std::vector read_owned_atoms(std::ifstream& ifs, StruMetadata& metadata, - double cutoff_bohr, - double skin_bohr, + const ModuleBase::Matrix3& primitive_latvec, + const ModuleBase::Matrix3& primitive_gt, + const std::vector& replicate, + double cutoff, + double skin, int& nat) { int rank = 0; #ifdef __MPI DomainDecomposition decomposition; - decomposition.init(MPI_COMM_WORLD, metadata.latvec, metadata.lat0, cutoff_bohr, skin_bohr); + decomposition.init(MPI_COMM_WORLD, metadata.latvec, metadata.lat0, cutoff, skin); MPI_Comm_rank(MPI_COMM_WORLD, &rank); #endif + int begin[3] = {0, 0, 0}; + int end[3] = {replicate[0], replicate[1], replicate[2]}; +#ifdef __MPI + const std::array& dims = decomposition.dims(); + const std::array& coords = decomposition.coords(); + for (int idim = 0; idim < 3; ++idim) + { + begin[idim] = std::max(0, static_cast(std::floor( + static_cast(coords[idim]) * replicate[idim] / dims[idim])) - 1); + end[idim] = std::min(replicate[idim], static_cast(std::ceil( + static_cast(coords[idim] + 1) * replicate[idim] / dims[idim])) + 1); + } +#endif + expect_keyword(ifs, "ATOMIC_POSITIONS"); const std::string coord_type = next_data_line(ifs, "ATOMIC_POSITIONS type"); const bool is_cartesian = coord_type == "Cartesian"; const bool is_direct = coord_type == "Direct"; if (!is_cartesian && !is_direct) { - throw std::runtime_error("Only Direct and Cartesian ATOMIC_POSITIONS are supported for LJ MD."); + throw std::runtime_error("Only Direct and Cartesian ATOMIC_POSITIONS are supported for MD."); } std::vector owned_atoms; @@ -231,13 +250,13 @@ std::vector read_owned_atoms(std::ifstream& ifs, if (is_cartesian) { cart.set(c1, c2, c3); - frac = wrap_fractional(cart * metadata.gt); - cart = frac * metadata.latvec; + frac = wrap_fractional(cart * primitive_gt); + cart = frac * primitive_latvec; } else { frac = wrap_fractional(ModuleBase::Vector3(c1, c2, c3)); - cart = frac * metadata.latvec; + cart = frac * primitive_latvec; } ModuleBase::Vector3 mbl(1, 1, 1); @@ -267,37 +286,52 @@ std::vector read_owned_atoms(std::ifstream& ifs, } } - int owner = 0; + for (int ix = begin[0]; ix < end[0]; ++ix) + { + for (int iy = begin[1]; iy < end[1]; ++iy) + { + for (int iz = begin[2]; iz < end[2]; ++iz) + { + ModuleBase::Vector3 final_frac( + (ix + frac.x) / replicate[0], + (iy + frac.y) / replicate[1], + (iz + frac.z) / replicate[2]); + int owner = 0; #ifdef __MPI - owner = decomposition.owner_rank_from_frac(frac); + owner = decomposition.owner_rank_from_frac(final_frac); #endif - if (owner == rank) - { - owned_atoms.push_back(LocalAtom(cart, - frac, - vel, - ModuleBase::Vector3(0.0, 0.0, 0.0), - mbl, - metadata.masses[it] / ModuleBase::AU_to_MASS, - static_cast(it), - ia, - owner, - false)); + if (owner == rank) + { + owned_atoms.push_back(LocalAtom(final_frac * metadata.latvec, + final_frac, + vel, + ModuleBase::Vector3(0.0, 0.0, 0.0), + mbl, + metadata.masses[it] / ModuleBase::AU_to_MASS, + static_cast(it), + ((ix * replicate[1] + iy) * replicate[2] + iz) * nat_type + ia, + owner, + false)); + } + } + } } - ++nat; } + metadata.stru_metadata.species[it].atom_count = nat_type * replicate[0] * replicate[1] * replicate[2]; + nat += metadata.stru_metadata.species[it].atom_count; } return owned_atoms; } } // namespace -MdCell DistributedMdCellReader::read_lj_stru(const std::string& stru_file, - double cutoff_bohr, - double skin_bohr) +MDCell DistributedMDCellReader::read_stru(const std::string& stru_file, + const std::vector& replicate, + double cutoff, + double skin) { - if (cutoff_bohr <= 0.0) + if (cutoff <= 0.0) { - throw std::runtime_error("MdCell requires a positive LJ cutoff from Parameter."); + throw std::runtime_error("MDCell requires a positive cutoff."); } std::ifstream ifs(stru_file.c_str(), std::ios::in); @@ -306,10 +340,22 @@ MdCell DistributedMdCellReader::read_lj_stru(const std::string& stru_file, throw std::runtime_error("Failed to open STRU file: " + stru_file); } + if (replicate.size() != 3 || replicate[0] <= 0 || replicate[1] <= 0 || replicate[2] <= 0) + { + throw std::runtime_error("replicate requires three positive integers."); + } StruMetadata metadata = parse_stru_metadata(ifs); + const ModuleBase::Matrix3 primitive_latvec = metadata.latvec; + const ModuleBase::Matrix3 primitive_gt = metadata.gt; + metadata.latvec.e11 *= replicate[0]; metadata.latvec.e12 *= replicate[0]; metadata.latvec.e13 *= replicate[0]; + metadata.latvec.e21 *= replicate[1]; metadata.latvec.e22 *= replicate[1]; metadata.latvec.e23 *= replicate[1]; + metadata.latvec.e31 *= replicate[2]; metadata.latvec.e32 *= replicate[2]; metadata.latvec.e33 *= replicate[2]; + metadata.gt = metadata.latvec.Inverse(); + metadata.omega = std::abs(metadata.latvec.Det()) * metadata.lat0 * metadata.lat0 * metadata.lat0; int nat = 0; - const std::vector owned_atoms = read_owned_atoms(ifs, metadata, cutoff_bohr, skin_bohr, nat); - MdCell mdcell(metadata.latvec, + const std::vector owned_atoms = read_owned_atoms(ifs, metadata, primitive_latvec, primitive_gt, + replicate, cutoff, skin, nat); + MDCell mdcell(metadata.latvec, metadata.gt, metadata.lat0, metadata.omega, @@ -317,8 +363,9 @@ MdCell DistributedMdCellReader::read_lj_stru(const std::string& stru_file, owned_atoms, metadata.labels, metadata.masses, - cutoff_bohr, - skin_bohr); + cutoff, + skin); mdcell.set_stru_metadata(metadata.stru_metadata); + mdcell.set_uses_replicated_stru(replicate[0] != 1 || replicate[1] != 1 || replicate[2] != 1); return mdcell; } diff --git a/source/source_cell/distributed_mdcell_reader.h b/source/source_cell/distributed_mdcell_reader.h index 1af6fa66708..bd2e8177e32 100644 --- a/source/source_cell/distributed_mdcell_reader.h +++ b/source/source_cell/distributed_mdcell_reader.h @@ -2,15 +2,17 @@ #define DISTRIBUTED_MDCELL_READER_H #include +#include -class MdCell; +class MDCell; -class DistributedMdCellReader +class DistributedMDCellReader { public: - static MdCell read_lj_stru(const std::string& stru_file, - double cutoff_bohr, - double skin_bohr); + static MDCell read_stru(const std::string& stru_file, + const std::vector& replicate, + double cutoff, + double skin); }; #endif diff --git a/source/source_cell/md_cell.cpp b/source/source_cell/md_cell.cpp index 6ae4fd138c8..babb85b46d3 100644 --- a/source/source_cell/md_cell.cpp +++ b/source/source_cell/md_cell.cpp @@ -1,13 +1,11 @@ #include "source_cell/md_cell.h" #include "source_cell/unitcell.h" -#include "source_io/module_parameter/parameter.h" -#include #include #include -double MdCell::wrap_fractional_(double value) +double MDCell::wrap_fractional_(double value) { value -= std::floor(value); if (value >= 1.0 - 1.0e-12 || value < 1.0e-12) @@ -17,18 +15,7 @@ double MdCell::wrap_fractional_(double value) return value; } -double MdCell::infer_cutoff_from_parameter_(const Parameter& param) -{ - double cutoff = 0.0; - const std::vector& lj_rcut = param.inp.mdp.lj_rcut; - for (std::size_t i = 0; i < lj_rcut.size(); ++i) - { - cutoff = std::max(cutoff, lj_rcut[i] * ModuleBase::ANGSTROM_AU); - } - return cutoff; -} - -void MdCell::clear_forces_(std::vector& atoms) +void MDCell::clear_forces_(std::vector& atoms) { for (std::size_t i = 0; i < atoms.size(); ++i) { @@ -36,7 +23,7 @@ void MdCell::clear_forces_(std::vector& atoms) } } -void MdCell::sync_backing_unitcell_geometry_() +void MDCell::sync_backing_unitcell_geometry_() { if (backing_unitcell_ == nullptr) { @@ -58,7 +45,7 @@ void MdCell::sync_backing_unitcell_geometry_() backing_unitcell_->cell_parameter_updated = true; } -void MdCell::sync_backing_unitcell_owned_atoms_() +void MDCell::sync_backing_unitcell_owned_atoms_() { if (backing_unitcell_ == nullptr) { @@ -75,7 +62,7 @@ void MdCell::sync_backing_unitcell_owned_atoms_() } } -void MdCell::initialize_from_ucell_serial_(UnitCell& ucell, double cutoff, double skin) +void MDCell::initialize_from_ucell_serial_(UnitCell& ucell, double cutoff, double skin) { backing_unitcell_ = &ucell; nat_ = ucell.nat; @@ -83,16 +70,18 @@ void MdCell::initialize_from_ucell_serial_(UnitCell& ucell, double cutoff, doubl omega_ = ucell.omega; latvec_ = ucell.latvec; gt_ = ucell.GT; - type_labels_.resize(static_cast(ucell.ntype)); - type_masses_.resize(static_cast(ucell.ntype)); + type_labels_.clear(); + type_masses_.clear(); + type_labels_.reserve(static_cast(ucell.ntype)); + type_masses_.reserve(static_cast(ucell.ntype)); stru_metadata_.species.resize(static_cast(ucell.ntype)); for (int it = 0; it < ucell.ntype; ++it) { MdStruSpecies& species = stru_metadata_.species[static_cast(it)]; species.label = ucell.atoms[it].label; species.mass = ucell.atoms[it].mass; - type_labels_[static_cast(it)] = species.label; - type_masses_[static_cast(it)] = species.mass; + type_labels_.push_back(species.label); + type_masses_.push_back(species.mass); if (static_cast(it) < ucell.pseudo_fn.size()) species.pseudo_file = ucell.pseudo_fn[it]; if (static_cast(it) < ucell.pseudo_type.size()) species.pseudo_type = ucell.pseudo_type[it]; if (static_cast(it) < ucell.orbital_fn.size()) species.orbital_file = ucell.orbital_fn[it]; @@ -126,7 +115,7 @@ void MdCell::initialize_from_ucell_serial_(UnitCell& ucell, double cutoff, doubl } #ifdef __MPI -void MdCell::initialize_from_ucell_(UnitCell& ucell, MPI_Comm comm, double cutoff, double skin) +void MDCell::initialize_from_ucell_(UnitCell& ucell, MPI_Comm comm, double cutoff, double skin) { backing_unitcell_ = &ucell; nat_ = ucell.nat; @@ -134,16 +123,18 @@ void MdCell::initialize_from_ucell_(UnitCell& ucell, MPI_Comm comm, double cutof omega_ = ucell.omega; latvec_ = ucell.latvec; gt_ = ucell.GT; - type_labels_.resize(static_cast(ucell.ntype)); - type_masses_.resize(static_cast(ucell.ntype)); + type_labels_.clear(); + type_masses_.clear(); + type_labels_.reserve(static_cast(ucell.ntype)); + type_masses_.reserve(static_cast(ucell.ntype)); stru_metadata_.species.resize(static_cast(ucell.ntype)); for (int it = 0; it < ucell.ntype; ++it) { MdStruSpecies& species = stru_metadata_.species[static_cast(it)]; species.label = ucell.atoms[it].label; species.mass = ucell.atoms[it].mass; - type_labels_[static_cast(it)] = species.label; - type_masses_[static_cast(it)] = species.mass; + type_labels_.push_back(species.label); + type_masses_.push_back(species.mass); if (static_cast(it) < ucell.pseudo_fn.size()) species.pseudo_file = ucell.pseudo_fn[it]; if (static_cast(it) < ucell.pseudo_type.size()) species.pseudo_type = ucell.pseudo_type[it]; if (static_cast(it) < ucell.orbital_fn.size()) species.orbital_file = ucell.orbital_fn[it]; @@ -167,7 +158,7 @@ void MdCell::initialize_from_ucell_(UnitCell& ucell, MPI_Comm comm, double cutof exchange_ghost_atoms(); } -void MdCell::initialize_from_owned_atoms_(MPI_Comm comm, double cutoff, double skin) +void MDCell::initialize_from_owned_atoms_(MPI_Comm comm, double cutoff, double skin) { comm_ = comm; cutoff_ = cutoff; @@ -182,7 +173,7 @@ void MdCell::initialize_from_owned_atoms_(MPI_Comm comm, double cutoff, double s #endif #ifndef __MPI -void MdCell::initialize_from_owned_atoms_(double cutoff, double skin) +void MDCell::initialize_from_owned_atoms_(double cutoff, double skin) { cutoff_ = cutoff; skin_ = skin; @@ -191,17 +182,16 @@ void MdCell::initialize_from_owned_atoms_(double cutoff, double skin) } #endif -MdCell::MdCell(UnitCell& ucell, const Parameter& param) +MDCell::MDCell(UnitCell& ucell, double cutoff, double skin) { - const double cutoff = infer_cutoff_from_parameter_(param); #ifdef __MPI - initialize_from_ucell_(ucell, MPI_COMM_WORLD, cutoff, 0.0); + initialize_from_ucell_(ucell, MPI_COMM_WORLD, cutoff, skin); #else - initialize_from_ucell_serial_(ucell, cutoff, 0.0); + initialize_from_ucell_serial_(ucell, cutoff, skin); #endif } -MdCell::MdCell(const ModuleBase::Matrix3& latvec, +MDCell::MDCell(const ModuleBase::Matrix3& latvec, const ModuleBase::Matrix3& gt, double lat0, double omega, @@ -229,7 +219,7 @@ MdCell::MdCell(const ModuleBase::Matrix3& latvec, } #ifdef __MPI -MdCell::MdCell(const ModuleBase::Matrix3& latvec, +MDCell::MDCell(const ModuleBase::Matrix3& latvec, const ModuleBase::Matrix3& gt, double lat0, double omega, @@ -253,33 +243,33 @@ MdCell::MdCell(const ModuleBase::Matrix3& latvec, initialize_from_owned_atoms_(comm, cutoff, skin); } -MdCell::MdCell(UnitCell& ucell, MPI_Comm comm, double cutoff, double skin) +MDCell::MDCell(UnitCell& ucell, MPI_Comm comm, double cutoff, double skin) { initialize_from_ucell_(ucell, comm, cutoff, skin); } -int MdCell::mpi_rank() const +int MDCell::mpi_rank() const { return rank_; } -int MdCell::mpi_size() const +int MDCell::mpi_size() const { return size_; } -MPI_Comm MdCell::communicator() const +MPI_Comm MDCell::communicator() const { return comm_; } -const DomainDecomposition& MdCell::decomposition() const +const DomainDecomposition& MDCell::decomposition() const { return decomp_; } #endif -void MdCell::exchange_ghost_atoms() +void MDCell::exchange_ghost_atoms() { #ifdef __MPI decomp_.exchange_ghost_atoms(owned_atoms_, ghost_atoms_); @@ -309,7 +299,7 @@ void MdCell::exchange_ghost_atoms() const double volume = std::abs(a1.x * a2xa3.x + a1.y * a2xa3.y + a1.z * a2xa3.z); if (volume <= 0.0) { - throw std::runtime_error("MdCell requires a nonzero cell volume for periodic ghosts."); + throw std::runtime_error("MDCell requires a nonzero cell volume for periodic ghosts."); } const double search_radius = (cutoff_ + skin_) / lat0_; @@ -344,7 +334,28 @@ void MdCell::exchange_ghost_atoms() } } -void MdCell::migrate_owned_atoms() +void MDCell::accumulate_ghost_forces() +{ +#ifdef __MPI + decomp_.accumulate_ghost_forces(owned_atoms_, ghost_atoms_); +#else + for (std::size_t ighost = 0; ighost < ghost_atoms_.size(); ++ighost) + { + const LocalAtom& ghost = ghost_atoms_[ighost]; + for (std::size_t iowned = 0; iowned < owned_atoms_.size(); ++iowned) + { + LocalAtom& owned = owned_atoms_[iowned]; + if (owned.type == ghost.type && owned.type_index == ghost.type_index) + { + owned.force += ghost.force; + break; + } + } + } +#endif +} + +void MDCell::migrate_owned_atoms() { #ifdef __MPI decomp_.migrate_owned_atoms(owned_atoms_); @@ -361,9 +372,10 @@ void MdCell::migrate_owned_atoms() atom.cart = atom.frac * latvec_; } sync_backing_unitcell_owned_atoms_(); + exchange_ghost_atoms(); } -void MdCell::set_lattice_vectors(const ModuleBase::Matrix3& latvec) +void MDCell::set_lattice_vectors(const ModuleBase::Matrix3& latvec) { latvec_ = latvec; gt_ = latvec_.Inverse(); @@ -377,7 +389,7 @@ void MdCell::set_lattice_vectors(const ModuleBase::Matrix3& latvec) sync_backing_unitcell_geometry_(); } -void MdCell::refresh_cart_from_frac() +void MDCell::refresh_cart_from_frac() { for (std::size_t i = 0; i < owned_atoms_.size(); ++i) { @@ -390,94 +402,111 @@ void MdCell::refresh_cart_from_frac() exchange_ghost_atoms(); } -const std::vector& MdCell::owned_atoms() const +const std::vector& MDCell::owned_atoms() const { return owned_atoms_; } -const std::vector& MdCell::ghost_atoms() const +const std::vector& MDCell::ghost_atoms() const { return ghost_atoms_; } -const std::vector& MdCell::type_labels() const +const std::vector& MDCell::type_labels() const { return type_labels_; } -const std::vector& MdCell::type_masses() const +const std::vector& MDCell::type_masses() const { return type_masses_; } -const MdStruMetadata& MdCell::stru_metadata() const +const MdStruMetadata& MDCell::stru_metadata() const { return stru_metadata_; } -void MdCell::set_stru_metadata(const MdStruMetadata& metadata) +void MDCell::set_stru_metadata(const MdStruMetadata& metadata) { stru_metadata_ = metadata; } -std::vector& MdCell::mutable_owned_atoms() +std::vector& MDCell::mutable_owned_atoms() { return owned_atoms_; } -std::vector& MdCell::mutable_ghost_atoms() +std::vector& MDCell::mutable_ghost_atoms() { return ghost_atoms_; } -int MdCell::nlocal() const +void MDCell::replace_owned_atoms_for_restart(const std::vector& owned_atoms) +{ + owned_atoms_ = owned_atoms; + clear_forces_(owned_atoms_); + exchange_ghost_atoms(); +} + +int MDCell::nlocal() const { return static_cast(owned_atoms_.size()); } -int MdCell::nghost() const +int MDCell::nghost() const { return static_cast(ghost_atoms_.size()); } -bool MdCell::init_vel() const +bool MDCell::init_vel() const { return init_vel_; } -void MdCell::set_init_vel(bool init_vel) +void MDCell::set_init_vel(bool init_vel) { init_vel_ = init_vel; } -double MdCell::cutoff() const +double MDCell::cutoff() const { return cutoff_; } -double MdCell::skin() const +double MDCell::skin() const { return skin_; } -bool MdCell::has_backing_unitcell() const +bool MDCell::has_backing_unitcell() const { return backing_unitcell_ != nullptr; } -UnitCell& MdCell::backing_unitcell() +bool MDCell::uses_replicated_stru() const +{ + return uses_replicated_stru_; +} + +void MDCell::set_uses_replicated_stru(const bool uses_replicated_stru) +{ + uses_replicated_stru_ = uses_replicated_stru; +} + +UnitCell& MDCell::backing_unitcell() { assert(backing_unitcell_ != nullptr); return *backing_unitcell_; } -const UnitCell& MdCell::backing_unitcell() const +const UnitCell& MDCell::backing_unitcell() const { assert(backing_unitcell_ != nullptr); return *backing_unitcell_; } -void MdCell::sync_backing_unitcell() +void MDCell::sync_backing_unitcell() { if (backing_unitcell_ == nullptr) { @@ -533,7 +562,7 @@ void MdCell::sync_backing_unitcell() const int iat = type_offset[it] + ia; if (owner[iat] != 1) { - throw std::runtime_error("MdCell backing UnitCell atom ownership is invalid."); + throw std::runtime_error("MDCell backing UnitCell atom ownership is invalid."); } backing_unitcell_->atoms[it].tau[ia].set(cart[3 * iat], cart[3 * iat + 1], cart[3 * iat + 2]); backing_unitcell_->atoms[it].taud[ia].set(frac[3 * iat], frac[3 * iat + 1], frac[3 * iat + 2]); @@ -548,32 +577,32 @@ void MdCell::sync_backing_unitcell() sync_backing_unitcell_owned_atoms_(); } -BaseCell::Kind MdCell::get_kind() const +BaseCell::Kind MDCell::get_kind() const { return Kind::md_cell; } -int MdCell::get_nat() const +int MDCell::get_nat() const { return nat_; } -double MdCell::get_lat0() const +double MDCell::get_lat0() const { return lat0_; } -double MdCell::get_omega() const +double MDCell::get_omega() const { return omega_; } -const ModuleBase::Matrix3& MdCell::get_latvec() const +const ModuleBase::Matrix3& MDCell::get_latvec() const { return latvec_; } -const ModuleBase::Matrix3& MdCell::get_GT() const +const ModuleBase::Matrix3& MDCell::get_GT() const { return gt_; } diff --git a/source/source_cell/md_cell.h b/source/source_cell/md_cell.h index 27fe453a22d..48643ad1c08 100644 --- a/source/source_cell/md_cell.h +++ b/source/source_cell/md_cell.h @@ -11,7 +11,6 @@ #include #include -class Parameter; class UnitCell; struct MdStruSpecies @@ -31,11 +30,11 @@ struct MdStruMetadata std::string descriptor_file; }; -class MdCell : public BaseCell +class MDCell : public BaseCell { public: - MdCell(UnitCell& ucell, const Parameter& param); - MdCell(const ModuleBase::Matrix3& latvec, + MDCell(UnitCell& ucell, double cutoff, double skin); + MDCell(const ModuleBase::Matrix3& latvec, const ModuleBase::Matrix3& gt, double lat0, double omega, @@ -47,7 +46,7 @@ class MdCell : public BaseCell double skin); #ifdef __MPI - MdCell(const ModuleBase::Matrix3& latvec, + MDCell(const ModuleBase::Matrix3& latvec, const ModuleBase::Matrix3& gt, double lat0, double omega, @@ -59,7 +58,7 @@ class MdCell : public BaseCell double cutoff, double skin); - MdCell(UnitCell& ucell, MPI_Comm comm, double cutoff, double skin); + MDCell(UnitCell& ucell, MPI_Comm comm, double cutoff, double skin); int mpi_rank() const; int mpi_size() const; @@ -69,6 +68,7 @@ class MdCell : public BaseCell #endif void exchange_ghost_atoms(); + void accumulate_ghost_forces(); void migrate_owned_atoms(); void set_lattice_vectors(const ModuleBase::Matrix3& latvec); void refresh_cart_from_frac(); @@ -81,6 +81,7 @@ class MdCell : public BaseCell void set_stru_metadata(const MdStruMetadata& metadata); std::vector& mutable_owned_atoms(); std::vector& mutable_ghost_atoms(); + void replace_owned_atoms_for_restart(const std::vector& owned_atoms); int nlocal() const; int nghost() const; @@ -89,6 +90,8 @@ class MdCell : public BaseCell double cutoff() const; double skin() const; bool has_backing_unitcell() const; + bool uses_replicated_stru() const; + void set_uses_replicated_stru(bool uses_replicated_stru); UnitCell& backing_unitcell(); const UnitCell& backing_unitcell() const; void sync_backing_unitcell(); @@ -101,7 +104,6 @@ class MdCell : public BaseCell const ModuleBase::Matrix3& get_latvec() const override; const ModuleBase::Matrix3& get_GT() const override; - static double infer_cutoff_from_parameter_(const Parameter& param); #ifdef __MPI void initialize_from_ucell_(UnitCell& ucell, MPI_Comm comm, double cutoff, double skin); #endif @@ -130,6 +132,7 @@ class MdCell : public BaseCell double cutoff_ = 0.0; double skin_ = 0.0; UnitCell* backing_unitcell_ = nullptr; + bool uses_replicated_stru_ = false; #ifdef __MPI MPI_Comm comm_ = MPI_COMM_NULL; diff --git a/source/source_cell/module_neighlist/domain_decomposition.cpp b/source/source_cell/module_neighlist/domain_decomposition.cpp index 443a7d6313b..bcd5b003295 100644 --- a/source/source_cell/module_neighlist/domain_decomposition.cpp +++ b/source/source_cell/module_neighlist/domain_decomposition.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include DomainDecomposition::DomainDecomposition() @@ -536,6 +537,97 @@ void DomainDecomposition::exchange_ghost_atoms(const std::vector& own } } +void DomainDecomposition::accumulate_ghost_forces(std::vector& owned_atoms, + std::vector& ghost_atoms) const +{ + std::map, std::size_t> owned_lookup; + for (std::size_t iat = 0; iat < owned_atoms.size(); ++iat) + { + const LocalAtom& atom = owned_atoms[iat]; + owned_lookup[std::make_pair(atom.type, atom.type_index)] = iat; + } + + std::vector > send_buffers(static_cast(size_)); + for (std::size_t iat = 0; iat < ghost_atoms.size(); ++iat) + { + LocalAtom& atom = ghost_atoms[iat]; + if (atom.owner_rank == rank_) + { + const std::map, std::size_t>::const_iterator found + = owned_lookup.find(std::make_pair(atom.type, atom.type_index)); + if (found == owned_lookup.end()) + { + throw std::runtime_error("Cannot match a local ghost force to an owned atom."); + } + owned_atoms[found->second].force += atom.force; + } + else + { + ForceRecord record; + record.type = atom.type; + record.type_index = atom.type_index; + record.force[0] = atom.force.x; + record.force[1] = atom.force.y; + record.force[2] = atom.force.z; + send_buffers[static_cast(atom.owner_rank)].push_back(record); + } + atom.force.set(0.0, 0.0, 0.0); + } + + std::vector send_counts(static_cast(size_), 0); + std::vector recv_counts(static_cast(size_), 0); + for (int irank = 0; irank < size_; ++irank) + { + const std::size_t bytes = send_buffers[static_cast(irank)].size() * sizeof(ForceRecord); + if (bytes > static_cast(std::numeric_limits::max())) + { + throw std::overflow_error("DomainDecomposition ghost force message exceeds MPI int range."); + } + send_counts[static_cast(irank)] = static_cast(bytes); + } + MPI_Alltoall(&send_counts[0], 1, MPI_INT, &recv_counts[0], 1, MPI_INT, comm_); + + std::vector send_displs(static_cast(size_), 0); + std::vector recv_displs(static_cast(size_), 0); + int total_send_bytes = 0; + int total_recv_bytes = 0; + for (int irank = 0; irank < size_; ++irank) + { + send_displs[static_cast(irank)] = total_send_bytes; + recv_displs[static_cast(irank)] = total_recv_bytes; + total_send_bytes += send_counts[static_cast(irank)]; + total_recv_bytes += recv_counts[static_cast(irank)]; + } + + std::vector send_records; + send_records.reserve(static_cast(total_send_bytes / static_cast(sizeof(ForceRecord)))); + for (int irank = 0; irank < size_; ++irank) + { + const std::vector& records = send_buffers[static_cast(irank)]; + send_records.insert(send_records.end(), records.begin(), records.end()); + } + std::vector recv_records(static_cast(total_recv_bytes / static_cast(sizeof(ForceRecord)))); + MPI_Alltoallv(send_records.empty() ? NULL : reinterpret_cast(&send_records[0]), + &send_counts[0], &send_displs[0], MPI_BYTE, + recv_records.empty() ? NULL : reinterpret_cast(&recv_records[0]), + &recv_counts[0], &recv_displs[0], MPI_BYTE, comm_); + + for (std::size_t irecord = 0; irecord < recv_records.size(); ++irecord) + { + const ForceRecord& record = recv_records[irecord]; + const std::map, std::size_t>::const_iterator found + = owned_lookup.find(std::make_pair(record.type, record.type_index)); + if (found == owned_lookup.end()) + { + throw std::runtime_error("Cannot match a received ghost force to an owned atom."); + } + LocalAtom& atom = owned_atoms[found->second]; + atom.force.x += record.force[0]; + atom.force.y += record.force[1]; + atom.force.z += record.force[2]; + } +} + void DomainDecomposition::migrate_owned_atoms(std::vector& owned_atoms) const { std::vector > send_atoms(static_cast(size_)); diff --git a/source/source_cell/module_neighlist/domain_decomposition.h b/source/source_cell/module_neighlist/domain_decomposition.h index 3b00c3d56fe..30b3165d175 100644 --- a/source/source_cell/module_neighlist/domain_decomposition.h +++ b/source/source_cell/module_neighlist/domain_decomposition.h @@ -40,6 +40,8 @@ class DomainDecomposition void exchange_ghost_atoms(const std::vector& owned_atoms, std::vector& ghost_atoms) const; + void accumulate_ghost_forces(std::vector& owned_atoms, + std::vector& ghost_atoms) const; void migrate_owned_atoms(std::vector& owned_atoms) const; const std::array& dims() const; @@ -71,6 +73,13 @@ class DomainDecomposition int recv_rank; }; + struct ForceRecord + { + int type; + int type_index; + double force[3]; + }; + MPI_Comm comm_; MPI_Comm cart_comm_; bool owns_cart_comm_; diff --git a/source/source_cell/module_neighlist/neighbor_search.cpp b/source/source_cell/module_neighlist/neighbor_search.cpp index cd5c5125ab4..5ca3e292243 100644 --- a/source/source_cell/module_neighlist/neighbor_search.cpp +++ b/source/source_cell/module_neighlist/neighbor_search.cpp @@ -39,53 +39,61 @@ const NeighborList& NeighborSearch::get_neighbor_list() const { // ========== Main public interface ========== -void NeighborSearch::init_distributed(const std::vector& owned_atoms, - const std::vector& ghost_atoms, - double sr, - double lat0) +void NeighborSearch::init_from_mdcell_(const MDCell& cell, double sr) { inside_atoms_.clear(); ghost_atoms_.clear(); all_atoms_.clear(); bin_manager_.clear(); - search_radius_ = sr / lat0; - const std::size_t total_atoms = ModuleNeighList::checked_size_sum(owned_atoms.size(), - ghost_atoms.size(), + search_radius_ = sr / cell.lat0(); + + const std::size_t total_atoms = ModuleNeighList::checked_size_sum(cell.owned_atoms().size(), + cell.ghost_atoms().size(), "NeighborSearch distributed atom count"); if (total_atoms > static_cast(std::numeric_limits::max())) { throw std::overflow_error("NeighborSearch distributed atom count exceeds local atom index range."); } + all_atoms_.reserve(total_atoms); - inside_atoms_.reserve(owned_atoms.size()); - ghost_atoms_.reserve(ghost_atoms.size()); - for (std::size_t iat = 0; iat < owned_atoms.size(); ++iat) + inside_atoms_.reserve(cell.owned_atoms().size()); + ghost_atoms_.reserve(cell.ghost_atoms().size()); + + for (size_t iat = 0; iat < cell.owned_atoms().size(); ++iat) { - const LocalAtom& local = owned_atoms[iat]; - NeighborAtom atom(local.cart.x, local.cart.y, local.cart.z, local.type, local.type_index, - ModuleNeighList::checked_local_atom_index(all_atoms_.size(), "NeighborSearch owned atom id"), + const LocalAtom& local = cell.owned_atoms()[iat]; + NeighborAtom atom(local.cart.x, + local.cart.y, + local.cart.z, + local.type, + local.type_index, + ModuleNeighList::checked_local_atom_index(all_atoms_.size(), + "NeighborSearch owned atom id"), local.owner_rank); all_atoms_.push_back(atom); inside_atoms_.push_back(atom); } - for (std::size_t iat = 0; iat < ghost_atoms.size(); ++iat) + + for (size_t iat = 0; iat < cell.ghost_atoms().size(); ++iat) { - const LocalAtom& local = ghost_atoms[iat]; - NeighborAtom atom(local.cart.x, local.cart.y, local.cart.z, local.type, local.type_index, - ModuleNeighList::checked_local_atom_index(all_atoms_.size(), "NeighborSearch ghost atom id"), + const LocalAtom& local = cell.ghost_atoms()[iat]; + NeighborAtom atom(local.cart.x, + local.cart.y, + local.cart.z, + local.type, + local.type_index, + ModuleNeighList::checked_local_atom_index(all_atoms_.size(), + "NeighborSearch ghost atom id"), local.owner_rank); all_atoms_.push_back(atom); ghost_atoms_.push_back(atom); } - neighbor_list_.initialize(inside_atoms_.size(), - ModuleNeighList::checked_size_product(all_atoms_.size(), neighbor_reserve_factor, - "NeighborSearch page size")); -} -void NeighborSearch::init_from_mdcell_(const MdCell& cell, double sr) -{ - init_distributed(cell.owned_atoms(), cell.ghost_atoms(), sr, cell.lat0()); + const std::size_t page_size = ModuleNeighList::checked_size_product(all_atoms_.size(), + neighbor_reserve_factor, + "NeighborSearch page size"); + neighbor_list_.initialize(inside_atoms_.size(), page_size); } void NeighborSearch::init_from_unitcell_(const UnitCell& ucell, double sr) @@ -138,7 +146,7 @@ void NeighborSearch::init(BaseCell& cell, double sr) { if (cell.kind() == BaseCell::Kind::md_cell) { - MdCell& md_cell = static_cast(cell); + MDCell& md_cell = static_cast(cell); init_from_mdcell_(md_cell, sr); return; } diff --git a/source/source_cell/module_neighlist/neighbor_search.h b/source/source_cell/module_neighlist/neighbor_search.h index fae73246d83..b132b4815c9 100644 --- a/source/source_cell/module_neighlist/neighbor_search.h +++ b/source/source_cell/module_neighlist/neighbor_search.h @@ -7,7 +7,7 @@ #include "source_cell/module_neighlist/local_atom.h" #include "source_cell/base_cell.h" -class MdCell; +class MDCell; class UnitCell; /** @@ -48,11 +48,6 @@ class NeighborSearch */ void init(BaseCell& cell, double sr); - void init_distributed(const std::vector& owned_atoms, - const std::vector& ghost_atoms, - double sr, - double lat0); - /** * @brief Build the neighbor list for all inside atoms. * @@ -115,7 +110,7 @@ class NeighborSearch */ void init_from_unitcell_(const UnitCell& ucell, double sr); - void init_from_mdcell_(const MdCell& cell, double sr); + void init_from_mdcell_(const MDCell& cell, double sr); void check_expand_condition(const UnitCell& ucell, int& glayerX_minus, int& glayerX, int& glayerY_minus, int& glayerY, int& glayerZ_minus, int& glayerZ); diff --git a/source/source_cell/module_neighlist/test/CMakeLists.txt b/source/source_cell/module_neighlist/test/CMakeLists.txt index ae0b535c9f4..58008609493 100644 --- a/source/source_cell/module_neighlist/test/CMakeLists.txt +++ b/source/source_cell/module_neighlist/test/CMakeLists.txt @@ -43,7 +43,7 @@ if(ENABLE_MPI) ) target_link_libraries(MODULE_CELL_NEIGHBOR_mdcell_migrate_mpi PRIVATE - parameter base device MPI::MPI_CXX GTest::gtest_main GTest::gmock_main abacus::linalg_libs + parameter base device GTest::gtest_main GTest::gmock_main abacus::linalg_libs ) install(TARGETS MODULE_CELL_NEIGHBOR_mdcell_migrate_mpi DESTINATION ${CMAKE_BINARY_DIR}/tests) add_test(NAME MODULE_CELL_NEIGHBOR_mdcell_migrate_mpi diff --git a/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp b/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp index 7691e9d19ef..aaac2dc9fe9 100644 --- a/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp +++ b/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp @@ -50,7 +50,7 @@ ModuleBase::Matrix3 make_lattice() } } // namespace -TEST(DistributedMdCellReaderTest, ReadOwnedAtomsFromSTRUWithoutUnitCell) +TEST(DistributedMDCellReaderTest, ReadOwnedAtomsFromSTRUWithoutUnitCell) { int rank = 0; MPI_Comm_rank(MPI_COMM_WORLD, &rank); @@ -62,9 +62,10 @@ TEST(DistributedMdCellReaderTest, ReadOwnedAtomsFromSTRUWithoutUnitCell) } MPI_Barrier(MPI_COMM_WORLD); - MdCell mdcell = DistributedMdCellReader::read_lj_stru(stru_file, - 1.0 * ModuleBase::ANGSTROM_AU, - 0.0); + MDCell mdcell = DistributedMDCellReader::read_stru(stru_file, + std::vector{1, 1, 1}, + 1.0 * ModuleBase::ANGSTROM_AU, + 0.0); EXPECT_EQ(mdcell.type_labels().size(), 1U); EXPECT_EQ(mdcell.type_labels()[0], "He"); diff --git a/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp b/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp index 406b88e02e8..a7dc0be0a56 100644 --- a/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp +++ b/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp @@ -54,7 +54,7 @@ TEST(MdCellMigrateMpiTest, AtomCrossingDomainMigratesToNewOwner) rank, false)); } - MdCell mdcell(latvec, + MDCell mdcell(latvec, latvec.Inverse(), 1.0, 1.0, diff --git a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp index 761c37f46e5..846bb6dc8b8 100644 --- a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp +++ b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp @@ -46,7 +46,7 @@ ModuleBase::Matrix3 identity_lattice() return latvec; } -MdCell make_mdcell(const ModuleBase::Matrix3& latvec, +MDCell make_mdcell(const ModuleBase::Matrix3& latvec, const std::vector >& positions, double cutoff) { @@ -69,7 +69,7 @@ MdCell make_mdcell(const ModuleBase::Matrix3& latvec, rank, false)); } - return MdCell(latvec, + return MDCell(latvec, gt, 1.0, cell_volume(latvec), @@ -77,6 +77,7 @@ MdCell make_mdcell(const ModuleBase::Matrix3& latvec, owned_atoms, std::vector(1, "X"), std::vector(1, 1.0), + MPI_COMM_SELF, cutoff, 0.0); } @@ -98,7 +99,7 @@ TEST(NeighborSearchTest, MdCellTwoAtomsNeighbor) const ModuleBase::Matrix3 latvec = identity_lattice(); const std::vector > positions{{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}; - MdCell mdcell = make_mdcell(latvec, positions, 1.0); + MDCell mdcell = make_mdcell(latvec, positions, 1.0); NeighborSearch ns; ns.init(mdcell, 1.0); @@ -116,7 +117,7 @@ TEST(NeighborSearchTest, MdCellNoNeighbor) const ModuleBase::Matrix3 latvec = identity_lattice(); const std::vector > positions{{0.0, 0.0, 0.0}, {0.49, 0.0, 0.0}}; - MdCell mdcell = make_mdcell(latvec, positions, 0.1); + MDCell mdcell = make_mdcell(latvec, positions, 0.1); NeighborSearch ns; ns.init(mdcell, 0.1); @@ -134,7 +135,7 @@ TEST(NeighborSearchTest, MdCellInitBuildsOwnedAndGhostAtoms) const ModuleBase::Matrix3 latvec = identity_lattice(); const std::vector > positions{{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}; - MdCell mdcell = make_mdcell(latvec, positions, 1.0); + MDCell mdcell = make_mdcell(latvec, positions, 1.0); NeighborSearch ns; ns.init(mdcell, 1.0); @@ -154,7 +155,7 @@ TEST(NeighborSearchTest, MdCellNeighborIdsStayLocalToAllAtoms) const std::vector > positions{{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}, {0.0, 0.5, 0.0}}; - MdCell mdcell = make_mdcell(latvec, positions, 0.75); + MDCell mdcell = make_mdcell(latvec, positions, 0.75); NeighborSearch ns; ns.init(mdcell, 0.75); @@ -185,7 +186,7 @@ TEST(NeighborSearchTest, MdCellPreservesMdAtomStateAcrossOwnedAndGhostStorage) const ModuleBase::Matrix3 latvec = identity_lattice(); const std::vector > positions{{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}; - MdCell mdcell = make_mdcell(latvec, positions, 1.0); + MDCell mdcell = make_mdcell(latvec, positions, 1.0); ASSERT_EQ(mdcell.nlocal(), 2); std::vector& owned_atoms = mdcell.mutable_owned_atoms(); @@ -244,7 +245,7 @@ TEST(NeighborSearchTest, MdCellMigrateOwnedAtomsReassignsOwnership) const ModuleBase::Matrix3 latvec = identity_lattice(); const std::vector > positions{{0.1, 0.1, 0.1}}; - MdCell mdcell = make_mdcell(latvec, positions, 0.2); + MDCell mdcell = make_mdcell(latvec, positions, 0.2); ASSERT_EQ(mdcell.nlocal(), 1); mdcell.mutable_owned_atoms()[0].cart.set(1.2, -0.1, 0.1); From 87ffe33b2d08791fb298b08abfac5107299ef4a2 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Mon, 10 Aug 2026 20:13:17 +0800 Subject: [PATCH 04/10] fix: restore distributed MDCell build integration --- source/source_cell/CMakeLists.txt | 2 ++ .../source_cell/distributed_mdcell_reader.cpp | 1 - source/source_cell/md_cell.cpp | 17 ----------------- source/source_cell/md_cell.h | 4 ---- .../module_neighlist/test/CMakeLists.txt | 2 +- source/source_esolver/esolver_lj.cpp | 11 +++-------- 6 files changed, 6 insertions(+), 31 deletions(-) diff --git a/source/source_cell/CMakeLists.txt b/source/source_cell/CMakeLists.txt index 5351fa730b9..3aa11e61dfa 100644 --- a/source/source_cell/CMakeLists.txt +++ b/source/source_cell/CMakeLists.txt @@ -24,6 +24,7 @@ add_library( klist.cpp parallel_kpoints.cpp cell_index.cpp + cell_tools.cpp check_atomic_stru.cpp update_cell.cpp magnetism.cpp @@ -39,6 +40,7 @@ add_library( cal_nelec_nband.cpp read_pseudo.cpp cal_wfc.cpp + cal_ux.cpp ) if(ENABLE_COVERAGE) diff --git a/source/source_cell/distributed_mdcell_reader.cpp b/source/source_cell/distributed_mdcell_reader.cpp index 797f6dbb16e..cbf92d5b6c4 100644 --- a/source/source_cell/distributed_mdcell_reader.cpp +++ b/source/source_cell/distributed_mdcell_reader.cpp @@ -366,6 +366,5 @@ MDCell DistributedMDCellReader::read_stru(const std::string& stru_file, cutoff, skin); mdcell.set_stru_metadata(metadata.stru_metadata); - mdcell.set_uses_replicated_stru(replicate[0] != 1 || replicate[1] != 1 || replicate[2] != 1); return mdcell; } diff --git a/source/source_cell/md_cell.cpp b/source/source_cell/md_cell.cpp index babb85b46d3..18fcc66bc2e 100644 --- a/source/source_cell/md_cell.cpp +++ b/source/source_cell/md_cell.cpp @@ -442,13 +442,6 @@ std::vector& MDCell::mutable_ghost_atoms() return ghost_atoms_; } -void MDCell::replace_owned_atoms_for_restart(const std::vector& owned_atoms) -{ - owned_atoms_ = owned_atoms; - clear_forces_(owned_atoms_); - exchange_ghost_atoms(); -} - int MDCell::nlocal() const { return static_cast(owned_atoms_.size()); @@ -484,16 +477,6 @@ bool MDCell::has_backing_unitcell() const return backing_unitcell_ != nullptr; } -bool MDCell::uses_replicated_stru() const -{ - return uses_replicated_stru_; -} - -void MDCell::set_uses_replicated_stru(const bool uses_replicated_stru) -{ - uses_replicated_stru_ = uses_replicated_stru; -} - UnitCell& MDCell::backing_unitcell() { assert(backing_unitcell_ != nullptr); diff --git a/source/source_cell/md_cell.h b/source/source_cell/md_cell.h index 48643ad1c08..09d91a23170 100644 --- a/source/source_cell/md_cell.h +++ b/source/source_cell/md_cell.h @@ -81,7 +81,6 @@ class MDCell : public BaseCell void set_stru_metadata(const MdStruMetadata& metadata); std::vector& mutable_owned_atoms(); std::vector& mutable_ghost_atoms(); - void replace_owned_atoms_for_restart(const std::vector& owned_atoms); int nlocal() const; int nghost() const; @@ -90,8 +89,6 @@ class MDCell : public BaseCell double cutoff() const; double skin() const; bool has_backing_unitcell() const; - bool uses_replicated_stru() const; - void set_uses_replicated_stru(bool uses_replicated_stru); UnitCell& backing_unitcell(); const UnitCell& backing_unitcell() const; void sync_backing_unitcell(); @@ -132,7 +129,6 @@ class MDCell : public BaseCell double cutoff_ = 0.0; double skin_ = 0.0; UnitCell* backing_unitcell_ = nullptr; - bool uses_replicated_stru_ = false; #ifdef __MPI MPI_Comm comm_ = MPI_COMM_NULL; diff --git a/source/source_cell/module_neighlist/test/CMakeLists.txt b/source/source_cell/module_neighlist/test/CMakeLists.txt index 58008609493..ae0b535c9f4 100644 --- a/source/source_cell/module_neighlist/test/CMakeLists.txt +++ b/source/source_cell/module_neighlist/test/CMakeLists.txt @@ -43,7 +43,7 @@ if(ENABLE_MPI) ) target_link_libraries(MODULE_CELL_NEIGHBOR_mdcell_migrate_mpi PRIVATE - parameter base device GTest::gtest_main GTest::gmock_main abacus::linalg_libs + parameter base device MPI::MPI_CXX GTest::gtest_main GTest::gmock_main abacus::linalg_libs ) install(TARGETS MODULE_CELL_NEIGHBOR_mdcell_migrate_mpi DESTINATION ${CMAKE_BINARY_DIR}/tests) add_test(NAME MODULE_CELL_NEIGHBOR_mdcell_migrate_mpi diff --git a/source/source_esolver/esolver_lj.cpp b/source/source_esolver/esolver_lj.cpp index 3b499d1adde..7d6abe9ab5c 100644 --- a/source/source_esolver/esolver_lj.cpp +++ b/source/source_esolver/esolver_lj.cpp @@ -5,10 +5,10 @@ #include "source_io/module_output/output_log.h" #include "source_cell/module_neighlist/neighbor_types.h" #include "source_cell/module_neighlist/neighbor_search.h" +#include "source_cell/md_cell.h" #include "source_base/global_variable.h" #include "source_base/timer.h" #ifdef __MPI -#include "source_cell/module_neighlist/domain_decomposition.h" #include "source_base/parallel_reduce.h" #endif @@ -58,13 +58,8 @@ void ESolver_LJ::runner(BaseCell& cell, const int istep) { ModuleBase::timer::start("ESolverLJ", "mpi_total"); ModuleBase::timer::start("ESolverLJ", "neigh_init"); - DomainDecomposition decomp; - decomp.init(MPI_COMM_WORLD, ucell.latvec, ucell.lat0, search_radius, 0.0); - std::vector owned_atoms; - std::vector ghost_atoms; - decomp.split_owned_atoms_from_ucell(ucell, owned_atoms); - decomp.exchange_ghost_atoms(owned_atoms, ghost_atoms); - neighbor_search.init_distributed(owned_atoms, ghost_atoms, search_radius, ucell.lat0); + MDCell mdcell(ucell, MPI_COMM_WORLD, search_radius, 0.0); + neighbor_search.init(mdcell, search_radius); ModuleBase::timer::end("ESolverLJ", "neigh_init"); ModuleBase::timer::start("ESolverLJ", "neigh_bld"); neighbor_search.build_neighbors(); From 380660279a41ce92d7fe5c73c05ef99e6b54c972 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Mon, 10 Aug 2026 20:24:14 +0800 Subject: [PATCH 05/10] fix: preserve MDCell build and MPI ownership --- source/Makefile.Objects | 2 + source/source_cell/md_cell.h | 5 ++ .../module_neighlist/domain_decomposition.cpp | 49 +++++++++++++++++++ .../module_neighlist/domain_decomposition.h | 4 ++ .../test/distributed_mdcell_reader_test.cpp | 4 ++ 5 files changed, 64 insertions(+) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index c5c503b33bb..0fca0e5b2a0 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -211,6 +211,8 @@ OBJS_CELL=atom_pseudo.o\ read_pseudo.o\ cal_wfc.o\ cal_ux.o\ + distributed_mdcell_reader.o\ + md_cell.o\ OBJS_DEEPKS=LCAO_deepks.o\ deepks_basic.o\ diff --git a/source/source_cell/md_cell.h b/source/source_cell/md_cell.h index 09d91a23170..0f764e43e74 100644 --- a/source/source_cell/md_cell.h +++ b/source/source_cell/md_cell.h @@ -34,6 +34,11 @@ class MDCell : public BaseCell { public: MDCell(UnitCell& ucell, double cutoff, double skin); + MDCell(const MDCell&) = delete; + MDCell& operator=(const MDCell&) = delete; + MDCell(MDCell&&) = default; + MDCell& operator=(MDCell&&) = default; + MDCell(const ModuleBase::Matrix3& latvec, const ModuleBase::Matrix3& gt, double lat0, diff --git a/source/source_cell/module_neighlist/domain_decomposition.cpp b/source/source_cell/module_neighlist/domain_decomposition.cpp index bcd5b003295..8472d7d0565 100644 --- a/source/source_cell/module_neighlist/domain_decomposition.cpp +++ b/source/source_cell/module_neighlist/domain_decomposition.cpp @@ -40,6 +40,55 @@ DomainDecomposition::~DomainDecomposition() } } +DomainDecomposition::DomainDecomposition(DomainDecomposition&& other) noexcept + : comm_(other.comm_), + cart_comm_(other.cart_comm_), + owns_cart_comm_(other.owns_cart_comm_), + rank_(other.rank_), + size_(other.size_), + dims_(other.dims_), + coords_(other.coords_), + margin_(other.margin_), + latvec_(other.latvec_), + inv_latvec_(other.inv_latvec_), + lat0_(other.lat0_), + cutoff_(other.cutoff_), + skin_(other.skin_) +{ + other.comm_ = MPI_COMM_NULL; + other.cart_comm_ = MPI_COMM_NULL; + other.owns_cart_comm_ = false; +} + +DomainDecomposition& DomainDecomposition::operator=(DomainDecomposition&& other) noexcept +{ + if (this != &other) + { + if (owns_cart_comm_ && cart_comm_ != MPI_COMM_NULL) + { + MPI_Comm_free(&cart_comm_); + } + comm_ = other.comm_; + cart_comm_ = other.cart_comm_; + owns_cart_comm_ = other.owns_cart_comm_; + rank_ = other.rank_; + size_ = other.size_; + dims_ = other.dims_; + coords_ = other.coords_; + margin_ = other.margin_; + latvec_ = other.latvec_; + inv_latvec_ = other.inv_latvec_; + lat0_ = other.lat0_; + cutoff_ = other.cutoff_; + skin_ = other.skin_; + + other.comm_ = MPI_COMM_NULL; + other.cart_comm_ = MPI_COMM_NULL; + other.owns_cart_comm_ = false; + } + return *this; +} + double DomainDecomposition::wrap_fractional(double value) { value -= std::floor(value); diff --git a/source/source_cell/module_neighlist/domain_decomposition.h b/source/source_cell/module_neighlist/domain_decomposition.h index 30b3165d175..cd6d010881b 100644 --- a/source/source_cell/module_neighlist/domain_decomposition.h +++ b/source/source_cell/module_neighlist/domain_decomposition.h @@ -26,6 +26,10 @@ class DomainDecomposition public: DomainDecomposition(); ~DomainDecomposition(); + DomainDecomposition(const DomainDecomposition&) = delete; + DomainDecomposition& operator=(const DomainDecomposition&) = delete; + DomainDecomposition(DomainDecomposition&& other) noexcept; + DomainDecomposition& operator=(DomainDecomposition&& other) noexcept; void init(MPI_Comm comm, const ModuleBase::Matrix3& latvec, diff --git a/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp b/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp index aaac2dc9fe9..9a80826183b 100644 --- a/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp +++ b/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp @@ -9,6 +9,10 @@ #include #include #include +#include + +static_assert(!std::is_copy_constructible::value, + "MDCell must not copy MPI communicator ownership."); namespace { From 40314c99dee8f3343cca0f586a3a7e53a4ab4657 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Mon, 10 Aug 2026 23:21:02 +0800 Subject: [PATCH 06/10] refactor: simplify MDCell construction --- .../source_cell/distributed_mdcell_reader.cpp | 20 ++---- source/source_cell/md_cell.cpp | 67 ++----------------- source/source_cell/md_cell.h | 35 ---------- .../test/md_cell_migrate_mpi_test.cpp | 1 - .../test/neighbor_search_test.cpp | 51 ++++++-------- 5 files changed, 30 insertions(+), 144 deletions(-) diff --git a/source/source_cell/distributed_mdcell_reader.cpp b/source/source_cell/distributed_mdcell_reader.cpp index cbf92d5b6c4..85c7f220a86 100644 --- a/source/source_cell/distributed_mdcell_reader.cpp +++ b/source/source_cell/distributed_mdcell_reader.cpp @@ -27,7 +27,6 @@ struct StruMetadata ModuleBase::Matrix3 gt; std::vector labels; std::vector masses; - MdStruMetadata stru_metadata; }; std::string trim_copy(const std::string& value) @@ -126,15 +125,15 @@ StruMetadata parse_stru_metadata(std::ifstream& ifs) } if (line == "NUMERICAL_ORBITAL") { - for (std::size_t it = 0; it < metadata.stru_metadata.species.size(); ++it) + for (std::size_t it = 0; it < metadata.labels.size(); ++it) { - metadata.stru_metadata.species[it].orbital_file = next_data_line(ifs, "NUMERICAL_ORBITAL body"); + next_data_line(ifs, "NUMERICAL_ORBITAL body"); } continue; } if (line == "NUMERICAL_DESCRIPTOR") { - metadata.stru_metadata.descriptor_file = next_data_line(ifs, "NUMERICAL_DESCRIPTOR body"); + next_data_line(ifs, "NUMERICAL_DESCRIPTOR body"); continue; } @@ -149,11 +148,6 @@ StruMetadata parse_stru_metadata(std::ifstream& ifs) metadata.labels.push_back(label); metadata.masses.push_back(parse_double(mass_token, "atomic mass")); - MdStruSpecies species; - species.label = label; - species.mass = metadata.masses.back(); - iss >> species.pseudo_file >> species.pseudo_type; - metadata.stru_metadata.species.push_back(species); } expect_keyword(ifs, "LATTICE_CONSTANT"); @@ -228,10 +222,8 @@ std::vector read_owned_atoms(std::ifstream& ifs, { throw std::runtime_error("ATOMIC_POSITIONS label order does not match ATOMIC_SPECIES."); } - std::istringstream magnetism(next_data_line(ifs, "magnetism")); - magnetism >> metadata.stru_metadata.species[it].start_mag; + next_data_line(ifs, "magnetism"); const int nat_type = parse_int(next_data_line(ifs, "atom count"), "atom count"); - metadata.stru_metadata.species[it].atom_count = nat_type; for (int ia = 0; ia < nat_type; ++ia) { @@ -317,8 +309,7 @@ std::vector read_owned_atoms(std::ifstream& ifs, } } } - metadata.stru_metadata.species[it].atom_count = nat_type * replicate[0] * replicate[1] * replicate[2]; - nat += metadata.stru_metadata.species[it].atom_count; + nat += nat_type * replicate[0] * replicate[1] * replicate[2]; } return owned_atoms; } @@ -365,6 +356,5 @@ MDCell DistributedMDCellReader::read_stru(const std::string& stru_file, metadata.masses, cutoff, skin); - mdcell.set_stru_metadata(metadata.stru_metadata); return mdcell; } diff --git a/source/source_cell/md_cell.cpp b/source/source_cell/md_cell.cpp index 18fcc66bc2e..b2c7f1f8839 100644 --- a/source/source_cell/md_cell.cpp +++ b/source/source_cell/md_cell.cpp @@ -74,21 +74,11 @@ void MDCell::initialize_from_ucell_serial_(UnitCell& ucell, double cutoff, doubl type_masses_.clear(); type_labels_.reserve(static_cast(ucell.ntype)); type_masses_.reserve(static_cast(ucell.ntype)); - stru_metadata_.species.resize(static_cast(ucell.ntype)); for (int it = 0; it < ucell.ntype; ++it) { - MdStruSpecies& species = stru_metadata_.species[static_cast(it)]; - species.label = ucell.atoms[it].label; - species.mass = ucell.atoms[it].mass; - type_labels_.push_back(species.label); - type_masses_.push_back(species.mass); - if (static_cast(it) < ucell.pseudo_fn.size()) species.pseudo_file = ucell.pseudo_fn[it]; - if (static_cast(it) < ucell.pseudo_type.size()) species.pseudo_type = ucell.pseudo_type[it]; - if (static_cast(it) < ucell.orbital_fn.size()) species.orbital_file = ucell.orbital_fn[it]; - if (static_cast(it) < ucell.magnet.start_mag.size()) species.start_mag = ucell.magnet.start_mag[it]; - species.atom_count = ucell.atoms[it].na; + type_labels_.push_back(ucell.atoms[it].label); + type_masses_.push_back(ucell.atoms[it].mass); } - stru_metadata_.descriptor_file = ucell.descriptor_file; init_vel_ = ucell.init_vel; cutoff_ = cutoff; skin_ = skin; @@ -127,21 +117,11 @@ void MDCell::initialize_from_ucell_(UnitCell& ucell, MPI_Comm comm, double cutof type_masses_.clear(); type_labels_.reserve(static_cast(ucell.ntype)); type_masses_.reserve(static_cast(ucell.ntype)); - stru_metadata_.species.resize(static_cast(ucell.ntype)); for (int it = 0; it < ucell.ntype; ++it) { - MdStruSpecies& species = stru_metadata_.species[static_cast(it)]; - species.label = ucell.atoms[it].label; - species.mass = ucell.atoms[it].mass; - type_labels_.push_back(species.label); - type_masses_.push_back(species.mass); - if (static_cast(it) < ucell.pseudo_fn.size()) species.pseudo_file = ucell.pseudo_fn[it]; - if (static_cast(it) < ucell.pseudo_type.size()) species.pseudo_type = ucell.pseudo_type[it]; - if (static_cast(it) < ucell.orbital_fn.size()) species.orbital_file = ucell.orbital_fn[it]; - if (static_cast(it) < ucell.magnet.start_mag.size()) species.start_mag = ucell.magnet.start_mag[it]; - species.atom_count = ucell.atoms[it].na; + type_labels_.push_back(ucell.atoms[it].label); + type_masses_.push_back(ucell.atoms[it].mass); } - stru_metadata_.descriptor_file = ucell.descriptor_file; init_vel_ = ucell.init_vel; comm_ = comm; cutoff_ = cutoff; @@ -219,35 +199,6 @@ MDCell::MDCell(const ModuleBase::Matrix3& latvec, } #ifdef __MPI -MDCell::MDCell(const ModuleBase::Matrix3& latvec, - const ModuleBase::Matrix3& gt, - double lat0, - double omega, - int nat, - const std::vector& owned_atoms, - const std::vector& type_labels, - const std::vector& type_masses, - MPI_Comm comm, - double cutoff, - double skin) -{ - latvec_ = latvec; - gt_ = gt; - lat0_ = lat0; - omega_ = omega; - nat_ = nat; - owned_atoms_ = owned_atoms; - type_labels_ = type_labels; - type_masses_ = type_masses; - init_vel_ = true; - initialize_from_owned_atoms_(comm, cutoff, skin); -} - -MDCell::MDCell(UnitCell& ucell, MPI_Comm comm, double cutoff, double skin) -{ - initialize_from_ucell_(ucell, comm, cutoff, skin); -} - int MDCell::mpi_rank() const { return rank_; @@ -422,16 +373,6 @@ const std::vector& MDCell::type_masses() const return type_masses_; } -const MdStruMetadata& MDCell::stru_metadata() const -{ - return stru_metadata_; -} - -void MDCell::set_stru_metadata(const MdStruMetadata& metadata) -{ - stru_metadata_ = metadata; -} - std::vector& MDCell::mutable_owned_atoms() { return owned_atoms_; diff --git a/source/source_cell/md_cell.h b/source/source_cell/md_cell.h index 0f764e43e74..1ce89e7dfec 100644 --- a/source/source_cell/md_cell.h +++ b/source/source_cell/md_cell.h @@ -12,24 +12,6 @@ #include class UnitCell; - -struct MdStruSpecies -{ - std::string label; - double mass = 0.0; - std::string pseudo_file; - std::string pseudo_type; - std::string orbital_file; - double start_mag = 0.0; - int atom_count = 0; -}; - -struct MdStruMetadata -{ - std::vector species; - std::string descriptor_file; -}; - class MDCell : public BaseCell { public: @@ -51,20 +33,6 @@ class MDCell : public BaseCell double skin); #ifdef __MPI - MDCell(const ModuleBase::Matrix3& latvec, - const ModuleBase::Matrix3& gt, - double lat0, - double omega, - int nat, - const std::vector& owned_atoms, - const std::vector& type_labels, - const std::vector& type_masses, - MPI_Comm comm, - double cutoff, - double skin); - - MDCell(UnitCell& ucell, MPI_Comm comm, double cutoff, double skin); - int mpi_rank() const; int mpi_size() const; MPI_Comm communicator() const; @@ -82,8 +50,6 @@ class MDCell : public BaseCell const std::vector& ghost_atoms() const; const std::vector& type_labels() const; const std::vector& type_masses() const; - const MdStruMetadata& stru_metadata() const; - void set_stru_metadata(const MdStruMetadata& metadata); std::vector& mutable_owned_atoms(); std::vector& mutable_ghost_atoms(); @@ -129,7 +95,6 @@ class MDCell : public BaseCell std::vector ghost_atoms_; std::vector type_labels_; std::vector type_masses_; - MdStruMetadata stru_metadata_; bool init_vel_ = false; double cutoff_ = 0.0; double skin_ = 0.0; diff --git a/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp b/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp index a7dc0be0a56..406d14c0bd7 100644 --- a/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp +++ b/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp @@ -62,7 +62,6 @@ TEST(MdCellMigrateMpiTest, AtomCrossingDomainMigratesToNewOwner) owned_atoms, std::vector(1, "X"), std::vector(1, 1.0), - MPI_COMM_WORLD, 0.1, 0.0); diff --git a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp index 846bb6dc8b8..64b56b53bdd 100644 --- a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp +++ b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp @@ -3,7 +3,9 @@ #include "source_cell/md_cell.h" #include "../neighbor_search.h" +#ifdef __MPI #include +#endif #include #include @@ -14,6 +16,7 @@ namespace { void ensure_mpi_initialized() { +#ifdef __MPI int initialized = 0; MPI_Initialized(&initialized); if (!initialized) @@ -21,14 +24,7 @@ void ensure_mpi_initialized() int provided = 0; MPI_Init_thread(NULL, NULL, MPI_THREAD_SINGLE, &provided); } -} - -double cell_volume(const ModuleBase::Matrix3& latvec) -{ - const double cx = latvec.e22 * latvec.e33 - latvec.e23 * latvec.e32; - const double cy = latvec.e23 * latvec.e31 - latvec.e21 * latvec.e33; - const double cz = latvec.e21 * latvec.e32 - latvec.e22 * latvec.e31; - return std::abs(latvec.e11 * cx + latvec.e12 * cy + latvec.e13 * cz); +#endif } ModuleBase::Matrix3 identity_lattice() @@ -46,40 +42,36 @@ ModuleBase::Matrix3 identity_lattice() return latvec; } +double cell_volume(const ModuleBase::Matrix3& latvec) +{ + const double cx = latvec.e22 * latvec.e33 - latvec.e23 * latvec.e32; + const double cy = latvec.e23 * latvec.e31 - latvec.e21 * latvec.e33; + const double cz = latvec.e21 * latvec.e32 - latvec.e22 * latvec.e31; + return std::abs(latvec.e11 * cx + latvec.e12 * cy + latvec.e13 * cz); +} + MDCell make_mdcell(const ModuleBase::Matrix3& latvec, const std::vector >& positions, double cutoff) { int rank = 0; - MPI_Comm_rank(MPI_COMM_SELF, &rank); - +#ifdef __MPI + MPI_Comm_rank(MPI_COMM_WORLD, &rank); +#endif const ModuleBase::Matrix3 gt = latvec.Inverse(); std::vector owned_atoms; owned_atoms.reserve(positions.size()); for (std::size_t iat = 0; iat < positions.size(); ++iat) { - owned_atoms.push_back(LocalAtom(positions[iat], - positions[iat] * gt, + owned_atoms.push_back(LocalAtom(positions[iat], positions[iat] * gt, ModuleBase::Vector3(0.0, 0.0, 0.0), ModuleBase::Vector3(0.0, 0.0, 0.0), - ModuleBase::Vector3(1, 1, 1), - 1.0, - 0, - static_cast(iat), - rank, - false)); + ModuleBase::Vector3(1, 1, 1), 1.0, 0, + static_cast(iat), rank, false)); } - return MDCell(latvec, - gt, - 1.0, - cell_volume(latvec), - static_cast(positions.size()), - owned_atoms, - std::vector(1, "X"), - std::vector(1, 1.0), - MPI_COMM_SELF, - cutoff, - 0.0); + return MDCell(latvec, gt, 1.0, cell_volume(latvec), static_cast(positions.size()), + owned_atoms, std::vector(1, "X"), + std::vector(1, 1.0), cutoff, 0.0); } std::size_t count_pairs(const NeighborList& list) @@ -140,7 +132,6 @@ TEST(NeighborSearchTest, MdCellInitBuildsOwnedAndGhostAtoms) NeighborSearch ns; ns.init(mdcell, 1.0); - EXPECT_EQ(mdcell.mpi_size(), 1); EXPECT_EQ(ns.get_inside_atoms().size(), 2U); EXPECT_GT(ns.get_ghost_atoms().size(), 0U); EXPECT_GT(ns.get_all_atoms().size(), ns.get_inside_atoms().size()); From 2f50429807e24f554d3969f7ec9c53f286e2866f Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Tue, 11 Aug 2026 00:17:52 +0800 Subject: [PATCH 07/10] refactor: pass MD communication domain --- source/source_base/CMakeLists.txt | 1 + source/source_base/communication_domain.cpp | 43 +++++++++++++++++++ source/source_base/communication_domain.h | 32 ++++++++++++++ .../source_cell/distributed_mdcell_reader.cpp | 16 ++++--- .../source_cell/distributed_mdcell_reader.h | 7 ++- source/source_cell/md_cell.cpp | 15 +++++-- source/source_cell/md_cell.h | 12 +++++- .../module_neighlist/test/CMakeLists.txt | 1 + .../test/distributed_mdcell_reader_test.cpp | 24 +++++++---- .../test/md_cell_migrate_mpi_test.cpp | 4 +- .../test/neighbor_search_test.cpp | 4 +- 11 files changed, 136 insertions(+), 23 deletions(-) create mode 100644 source/source_base/communication_domain.cpp create mode 100644 source/source_base/communication_domain.h diff --git a/source/source_base/CMakeLists.txt b/source/source_base/CMakeLists.txt index 596b3f3b463..a7b96c3f014 100644 --- a/source/source_base/CMakeLists.txt +++ b/source/source_base/CMakeLists.txt @@ -51,6 +51,7 @@ add_library( tool_title.cpp ylm.cpp parallel_common.cpp + communication_domain.cpp parallel_global.cpp parallel_comm.cpp parallel_reduce.cpp diff --git a/source/source_base/communication_domain.cpp b/source/source_base/communication_domain.cpp new file mode 100644 index 00000000000..946222bc3eb --- /dev/null +++ b/source/source_base/communication_domain.cpp @@ -0,0 +1,43 @@ +#include "source_base/communication_domain.h" + +namespace ModuleBase +{ +CommunicationDomain::CommunicationDomain() +{ +} + +#ifdef __MPI +CommunicationDomain::CommunicationDomain(MPI_Comm communicator) : communicator_(communicator) +{ + if (communicator_ != MPI_COMM_NULL) + { + MPI_Comm_rank(communicator_, &rank_); + MPI_Comm_size(communicator_, &size_); + } +} + +MPI_Comm CommunicationDomain::communicator() const +{ + return communicator_; +} +#endif + +int CommunicationDomain::rank() const +{ + return rank_; +} + +int CommunicationDomain::size() const +{ + return size_; +} + +CommunicationDomain world_communication_domain() +{ +#ifdef __MPI + return CommunicationDomain(MPI_COMM_WORLD); +#else + return CommunicationDomain(); +#endif +} +} // namespace ModuleBase diff --git a/source/source_base/communication_domain.h b/source/source_base/communication_domain.h new file mode 100644 index 00000000000..b91817db27e --- /dev/null +++ b/source/source_base/communication_domain.h @@ -0,0 +1,32 @@ +#ifndef COMMUNICATION_DOMAIN_H +#define COMMUNICATION_DOMAIN_H + +#ifdef __MPI +#include +#endif + +namespace ModuleBase +{ +class CommunicationDomain +{ +public: + CommunicationDomain(); +#ifdef __MPI + explicit CommunicationDomain(MPI_Comm communicator); + MPI_Comm communicator() const; +#endif + int rank() const; + int size() const; + +private: +#ifdef __MPI + MPI_Comm communicator_ = MPI_COMM_NULL; +#endif + int rank_ = 0; + int size_ = 1; +}; + +CommunicationDomain world_communication_domain(); +} // namespace ModuleBase + +#endif diff --git a/source/source_cell/distributed_mdcell_reader.cpp b/source/source_cell/distributed_mdcell_reader.cpp index 85c7f220a86..5640c54038c 100644 --- a/source/source_cell/distributed_mdcell_reader.cpp +++ b/source/source_cell/distributed_mdcell_reader.cpp @@ -1,6 +1,7 @@ #include "source_cell/distributed_mdcell_reader.h" #include "source_base/constants.h" +#include "source_base/communication_domain.h" #include "source_base/vector3.h" #include "source_cell/md_cell.h" @@ -181,13 +182,14 @@ std::vector read_owned_atoms(std::ifstream& ifs, const std::vector& replicate, double cutoff, double skin, - int& nat) + int& nat, + const ModuleBase::CommunicationDomain& communication_domain) { int rank = 0; #ifdef __MPI DomainDecomposition decomposition; - decomposition.init(MPI_COMM_WORLD, metadata.latvec, metadata.lat0, cutoff, skin); - MPI_Comm_rank(MPI_COMM_WORLD, &rank); + decomposition.init(communication_domain.communicator(), metadata.latvec, metadata.lat0, cutoff, skin); + rank = communication_domain.rank(); #endif int begin[3] = {0, 0, 0}; @@ -318,7 +320,8 @@ std::vector read_owned_atoms(std::ifstream& ifs, MDCell DistributedMDCellReader::read_stru(const std::string& stru_file, const std::vector& replicate, double cutoff, - double skin) + double skin, + const ModuleBase::CommunicationDomain& communication_domain) { if (cutoff <= 0.0) { @@ -345,7 +348,7 @@ MDCell DistributedMDCellReader::read_stru(const std::string& stru_file, metadata.omega = std::abs(metadata.latvec.Det()) * metadata.lat0 * metadata.lat0 * metadata.lat0; int nat = 0; const std::vector owned_atoms = read_owned_atoms(ifs, metadata, primitive_latvec, primitive_gt, - replicate, cutoff, skin, nat); + replicate, cutoff, skin, nat, communication_domain); MDCell mdcell(metadata.latvec, metadata.gt, metadata.lat0, @@ -355,6 +358,7 @@ MDCell DistributedMDCellReader::read_stru(const std::string& stru_file, metadata.labels, metadata.masses, cutoff, - skin); + skin, + communication_domain); return mdcell; } diff --git a/source/source_cell/distributed_mdcell_reader.h b/source/source_cell/distributed_mdcell_reader.h index bd2e8177e32..1d74c71ca88 100644 --- a/source/source_cell/distributed_mdcell_reader.h +++ b/source/source_cell/distributed_mdcell_reader.h @@ -5,6 +5,10 @@ #include class MDCell; +namespace ModuleBase +{ +class CommunicationDomain; +} class DistributedMDCellReader { @@ -12,7 +16,8 @@ class DistributedMDCellReader static MDCell read_stru(const std::string& stru_file, const std::vector& replicate, double cutoff, - double skin); + double skin, + const ModuleBase::CommunicationDomain& communication_domain); }; #endif diff --git a/source/source_cell/md_cell.cpp b/source/source_cell/md_cell.cpp index b2c7f1f8839..496e87d9b01 100644 --- a/source/source_cell/md_cell.cpp +++ b/source/source_cell/md_cell.cpp @@ -1,5 +1,6 @@ #include "source_cell/md_cell.h" +#include "source_base/communication_domain.h" #include "source_cell/unitcell.h" #include @@ -162,11 +163,15 @@ void MDCell::initialize_from_owned_atoms_(double cutoff, double skin) } #endif -MDCell::MDCell(UnitCell& ucell, double cutoff, double skin) +MDCell::MDCell(UnitCell& ucell, + double cutoff, + double skin, + const ModuleBase::CommunicationDomain& communication_domain) { #ifdef __MPI - initialize_from_ucell_(ucell, MPI_COMM_WORLD, cutoff, skin); + initialize_from_ucell_(ucell, communication_domain.communicator(), cutoff, skin); #else + static_cast(communication_domain); initialize_from_ucell_serial_(ucell, cutoff, skin); #endif } @@ -180,7 +185,8 @@ MDCell::MDCell(const ModuleBase::Matrix3& latvec, const std::vector& type_labels, const std::vector& type_masses, double cutoff, - double skin) + double skin, + const ModuleBase::CommunicationDomain& communication_domain) { latvec_ = latvec; gt_ = gt; @@ -192,8 +198,9 @@ MDCell::MDCell(const ModuleBase::Matrix3& latvec, type_masses_ = type_masses; init_vel_ = true; #ifdef __MPI - initialize_from_owned_atoms_(MPI_COMM_WORLD, cutoff, skin); + initialize_from_owned_atoms_(communication_domain.communicator(), cutoff, skin); #else + static_cast(communication_domain); initialize_from_owned_atoms_(cutoff, skin); #endif } diff --git a/source/source_cell/md_cell.h b/source/source_cell/md_cell.h index 1ce89e7dfec..0224ce6ccc9 100644 --- a/source/source_cell/md_cell.h +++ b/source/source_cell/md_cell.h @@ -12,10 +12,17 @@ #include class UnitCell; +namespace ModuleBase +{ +class CommunicationDomain; +} class MDCell : public BaseCell { public: - MDCell(UnitCell& ucell, double cutoff, double skin); + MDCell(UnitCell& ucell, + double cutoff, + double skin, + const ModuleBase::CommunicationDomain& communication_domain); MDCell(const MDCell&) = delete; MDCell& operator=(const MDCell&) = delete; MDCell(MDCell&&) = default; @@ -30,7 +37,8 @@ class MDCell : public BaseCell const std::vector& type_labels, const std::vector& type_masses, double cutoff, - double skin); + double skin, + const ModuleBase::CommunicationDomain& communication_domain); #ifdef __MPI int mpi_rank() const; diff --git a/source/source_cell/module_neighlist/test/CMakeLists.txt b/source/source_cell/module_neighlist/test/CMakeLists.txt index ae0b535c9f4..5c4f0336f9b 100644 --- a/source/source_cell/module_neighlist/test/CMakeLists.txt +++ b/source/source_cell/module_neighlist/test/CMakeLists.txt @@ -58,6 +58,7 @@ if(ENABLE_MPI) ../../md_cell.cpp ../domain_decomposition.cpp ../../../source_base/global_variable.cpp + ../../../source_base/communication_domain.cpp ../../../source_base/matrix.cpp ../../../source_base/matrix3.cpp ../../../source_base/tool_quit.cpp diff --git a/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp b/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp index 9a80826183b..6a4b9be6cab 100644 --- a/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp +++ b/source/source_cell/module_neighlist/test/distributed_mdcell_reader_test.cpp @@ -3,6 +3,7 @@ #include "source_cell/distributed_mdcell_reader.h" #include "source_cell/md_cell.h" #include "source_base/constants.h" +#include "source_base/communication_domain.h" #include "source_cell/module_neighlist/domain_decomposition.h" #include @@ -56,20 +57,25 @@ ModuleBase::Matrix3 make_lattice() TEST(DistributedMDCellReaderTest, ReadOwnedAtomsFromSTRUWithoutUnitCell) { - int rank = 0; - MPI_Comm_rank(MPI_COMM_WORLD, &rank); + int world_rank = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &world_rank); const std::string stru_file = "distributed_mdcell_reader_cartesian.STRU"; - if (rank == 0) + if (world_rank == 0) { write_cartesian_stru_case(stru_file); } MPI_Barrier(MPI_COMM_WORLD); + MPI_Comm md_comm = MPI_COMM_NULL; + MPI_Comm_split(MPI_COMM_WORLD, world_rank % 2, world_rank, &md_comm); + const ModuleBase::CommunicationDomain communication_domain(md_comm); + MDCell mdcell = DistributedMDCellReader::read_stru(stru_file, std::vector{1, 1, 1}, 1.0 * ModuleBase::ANGSTROM_AU, - 0.0); + 0.0, + communication_domain); EXPECT_EQ(mdcell.type_labels().size(), 1U); EXPECT_EQ(mdcell.type_labels()[0], "He"); @@ -78,18 +84,18 @@ TEST(DistributedMDCellReaderTest, ReadOwnedAtomsFromSTRUWithoutUnitCell) EXPECT_EQ(mdcell.nat(), 4); DomainDecomposition decomp; - decomp.init(MPI_COMM_WORLD, make_lattice(), 1.0, 1.0 * ModuleBase::ANGSTROM_AU, 0.0); + decomp.init(md_comm, make_lattice(), 1.0, 1.0 * ModuleBase::ANGSTROM_AU, 0.0); long long local_count = static_cast(mdcell.owned_atoms().size()); long long global_count = 0; - MPI_Allreduce(&local_count, &global_count, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&local_count, &global_count, 1, MPI_LONG_LONG, MPI_SUM, md_comm); EXPECT_EQ(global_count, 4); std::set > local_ids; for (std::size_t iat = 0; iat < mdcell.owned_atoms().size(); ++iat) { const LocalAtom& atom = mdcell.owned_atoms()[iat]; - EXPECT_EQ(decomp.owner_rank_from_frac(atom.frac), rank); + EXPECT_EQ(decomp.owner_rank_from_frac(atom.frac), communication_domain.rank()); local_ids.insert(std::make_pair(atom.type, atom.type_index)); EXPECT_GE(atom.type, 0); EXPECT_DOUBLE_EQ(atom.force.x, 0.0); @@ -124,9 +130,11 @@ TEST(DistributedMDCellReaderTest, ReadOwnedAtomsFromSTRUWithoutUnitCell) const int saw_flags[2] = {saw_v01 ? 1 : 0, saw_v04 ? 1 : 0}; int reduced_flags[2] = {0, 0}; - MPI_Allreduce(saw_flags, reduced_flags, 2, MPI_INT, MPI_MAX, MPI_COMM_WORLD); + MPI_Allreduce(saw_flags, reduced_flags, 2, MPI_INT, MPI_MAX, md_comm); EXPECT_EQ(reduced_flags[0], 1); EXPECT_EQ(reduced_flags[1], 1); + + MPI_Comm_free(&md_comm); } int main(int argc, char** argv) diff --git a/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp b/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp index 406d14c0bd7..7c418ad73e9 100644 --- a/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp +++ b/source/source_cell/module_neighlist/test/md_cell_migrate_mpi_test.cpp @@ -1,6 +1,7 @@ #include #include "source_cell/md_cell.h" +#include "source_base/communication_domain.h" #include @@ -63,7 +64,8 @@ TEST(MdCellMigrateMpiTest, AtomCrossingDomainMigratesToNewOwner) std::vector(1, "X"), std::vector(1, 1.0), 0.1, - 0.0); + 0.0, + ModuleBase::world_communication_domain()); ASSERT_EQ(mdcell.mpi_size(), size); if (size == 2) diff --git a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp index 64b56b53bdd..db4269beded 100644 --- a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp +++ b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp @@ -2,6 +2,7 @@ #include "source_cell/md_cell.h" #include "../neighbor_search.h" +#include "source_base/communication_domain.h" #ifdef __MPI #include @@ -71,7 +72,8 @@ MDCell make_mdcell(const ModuleBase::Matrix3& latvec, } return MDCell(latvec, gt, 1.0, cell_volume(latvec), static_cast(positions.size()), owned_atoms, std::vector(1, "X"), - std::vector(1, 1.0), cutoff, 0.0); + std::vector(1, 1.0), cutoff, 0.0, + ModuleBase::world_communication_domain()); } std::size_t count_pairs(const NeighborList& list) From 5436605d596f7a00eda313f27298e956ca79f483 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Tue, 11 Aug 2026 00:36:10 +0800 Subject: [PATCH 08/10] fix: pass communication domain to LJ solver --- source/source_esolver/esolver_lj.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/source/source_esolver/esolver_lj.cpp b/source/source_esolver/esolver_lj.cpp index 7d6abe9ab5c..3ddd664e9b2 100644 --- a/source/source_esolver/esolver_lj.cpp +++ b/source/source_esolver/esolver_lj.cpp @@ -6,6 +6,7 @@ #include "source_cell/module_neighlist/neighbor_types.h" #include "source_cell/module_neighlist/neighbor_search.h" #include "source_cell/md_cell.h" +#include "source_base/communication_domain.h" #include "source_base/global_variable.h" #include "source_base/timer.h" #ifdef __MPI @@ -58,7 +59,8 @@ void ESolver_LJ::runner(BaseCell& cell, const int istep) { ModuleBase::timer::start("ESolverLJ", "mpi_total"); ModuleBase::timer::start("ESolverLJ", "neigh_init"); - MDCell mdcell(ucell, MPI_COMM_WORLD, search_radius, 0.0); + const ModuleBase::CommunicationDomain communication_domain = ModuleBase::world_communication_domain(); + MDCell mdcell(ucell, search_radius, 0.0, communication_domain); neighbor_search.init(mdcell, search_radius); ModuleBase::timer::end("ESolverLJ", "neigh_init"); ModuleBase::timer::start("ESolverLJ", "neigh_bld"); From 5c46d95d1bfadab68ea0ab9a895e364d4600d0f3 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Thu, 13 Aug 2026 12:56:53 +0800 Subject: [PATCH 09/10] fix: restore neighbor search serial initialization --- .../module_neighlist/neighbor_search.cpp | 29 +------------------ .../test/neighbor_search_test.cpp | 1 - 2 files changed, 1 insertion(+), 29 deletions(-) diff --git a/source/source_cell/module_neighlist/neighbor_search.cpp b/source/source_cell/module_neighlist/neighbor_search.cpp index d8f5b648c93..5ca3e292243 100644 --- a/source/source_cell/module_neighlist/neighbor_search.cpp +++ b/source/source_cell/module_neighlist/neighbor_search.cpp @@ -4,7 +4,6 @@ #include #include -#include #include #include #include @@ -111,7 +110,6 @@ void NeighborSearch::init_from_unitcell_(const UnitCell& ucell, double sr) { for (int j = 0; j < ucell.atoms[i].na; j++) { - const ModuleBase::Vector3 position = ucell.get_tau(i, j); const ModuleNeighList::LocalAtomIndex atom_count = ModuleNeighList::checked_local_atom_index(all_atoms_.size(), "NeighborSearch atom id"); @@ -207,32 +205,7 @@ void NeighborSearch::set_member_variables(const UnitCell& ucell, int glayerX_min ModuleBase::Vector3 vec2(ucell.latvec.e21, ucell.latvec.e22, ucell.latvec.e23); ModuleBase::Vector3 vec3(ucell.latvec.e31, ucell.latvec.e32, ucell.latvec.e33); - const ModuleBase::Matrix3& lattice = ucell.get_latvec(); - const ModuleBase::Vector3 vec1(lattice.e11, lattice.e12, lattice.e13); - const ModuleBase::Vector3 vec2(lattice.e21, lattice.e22, lattice.e23); - const ModuleBase::Vector3 vec3(lattice.e31, lattice.e32, lattice.e33); - - const std::size_t image_count_x - = ModuleNeighList::checked_size_sum(static_cast(glayerX_minus), - static_cast(glayerX), - "NeighborSearch x image count"); - const std::size_t image_count_y - = ModuleNeighList::checked_size_sum(static_cast(glayerY_minus), - static_cast(glayerY), - "NeighborSearch y image count"); - const std::size_t image_count_z - = ModuleNeighList::checked_size_sum(static_cast(glayerZ_minus), - static_cast(glayerZ), - "NeighborSearch z image count"); - const std::size_t image_count_yz - = ModuleNeighList::checked_size_product(image_count_y, - image_count_z, - "NeighborSearch yz image count"); - const std::size_t image_count - = ModuleNeighList::checked_size_product(image_count_x, - image_count_yz, - "NeighborSearch periodic image count"); - if (image_count == 0) + for (int ix = -glayerX_minus; ix < glayerX; ix++) { for (int iy = -glayerY_minus; iy < glayerY; iy++) { diff --git a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp index 5d99ec234b2..ea2271e9c17 100644 --- a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp +++ b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp @@ -103,7 +103,6 @@ void expect_same_atoms(const std::vector& lhs, EXPECT_EQ(lhs[i].atom_type, rhs[i].atom_type) << "atom " << i; EXPECT_EQ(lhs[i].atom_index, rhs[i].atom_index) << "atom " << i; EXPECT_EQ(lhs[i].atom_id, rhs[i].atom_id) << "atom " << i; - EXPECT_EQ(lhs[i].global_id, rhs[i].global_id) << "atom " << i; EXPECT_EQ(lhs[i].owner_rank, rhs[i].owner_rank) << "atom " << i; } } From 675a08543f626f8cae543a6958d352529aff6fa0 Mon Sep 17 00:00:00 2001 From: Fei Yang <2501213217@stu.pku.edu.cn> Date: Thu, 13 Aug 2026 13:17:21 +0800 Subject: [PATCH 10/10] refactor: restore serial neighbor list implementation --- .../module_neighlist/bin_manager.cpp | 166 ++++------- .../module_neighlist/bin_manager.h | 15 - .../module_neighlist/test/CMakeLists.txt | 3 +- .../test/bin_manager_test.cpp | 114 -------- .../test/neighbor_search_test.cpp | 274 ++++-------------- 5 files changed, 116 insertions(+), 456 deletions(-) diff --git a/source/source_cell/module_neighlist/bin_manager.cpp b/source/source_cell/module_neighlist/bin_manager.cpp index 8cf5f2bc928..0cae41420fe 100644 --- a/source/source_cell/module_neighlist/bin_manager.cpp +++ b/source/source_cell/module_neighlist/bin_manager.cpp @@ -5,15 +5,6 @@ #include #include "bin_manager.h" -#ifdef _OPENMP -#include -#endif - -namespace -{ -constexpr int neighbor_build_openmp_threshold = 256; -} - // ========== Bin class implementation ========== const std::vector& Bin::get_atom_indices() const { @@ -185,63 +176,6 @@ int BinManager::bin_index(int ix, int iy, int iz) const { return ix * nbiny_ * nbinz_ + iy * nbinz_ + iz; } -template -void BinManager::visit_neighbors(const NeighborAtom& atom, - const std::vector& binned_atoms, - double sradius2, - const Emit& emit) const -{ - const int ix = std::min( - std::max(int((atom.position_x - x_min_) / bin_sizex_), 0), - nbinx_ - 1 - ); - - const int iy = std::min( - std::max(int((atom.position_y - y_min_) / bin_sizey_), 0), - nbiny_ - 1 - ); - - const int iz = std::min( - std::max(int((atom.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++) - { - const int jx = ix + dx; - const int jy = iy + dy; - const int jz = iz + dz; - - if (jx < 0 || jx >= nbinx_ || - jy < 0 || jy >= nbiny_ || - jz < 0 || jz >= nbinz_) - { - continue; - } - - const int nidx = bin_index(jx, jy, jz); - for (const ModuleNeighList::LocalAtomIndex binned_atom_index : bins_[nidx].get_atom_indices()) - { - const NeighborAtom& natom = binned_atoms[static_cast(binned_atom_index)]; - const double delta_x = atom.position_x - natom.position_x; - const double delta_y = atom.position_y - natom.position_y; - const double delta_z = atom.position_z - natom.position_z; - const double dist2 = delta_x * delta_x + delta_y * delta_y + delta_z * delta_z; - - if (natom.atom_id != atom.atom_id && dist2 <= sradius2) - { - emit(natom.atom_id); - } - } - } - } - } -} - void BinManager::build_atom_neighbors( NeighborList& neighbor_list, const std::vector& atoms, @@ -250,63 +184,71 @@ void BinManager::build_atom_neighbors( { assert(atoms.size() == static_cast(neighbor_list.get_nlocal())); - const double sradius2 = sradius_ * sradius_; + double sradius2 = sradius_ * sradius_; neighbor_list.reset(); - const int nlocal = neighbor_list.get_nlocal(); + std::vector neigh_tmp; -#ifdef _OPENMP - const bool use_parallel = nlocal >= neighbor_build_openmp_threshold && omp_get_max_threads() > 1; - if (use_parallel) + const int nlocal = neighbor_list.get_nlocal(); + for (int i = 0; i < nlocal; i++) { - std::vector neighbor_counts(static_cast(nlocal), 0); + neigh_tmp.clear(); + const NeighborAtom& atom = atoms[i]; -#pragma omp parallel for schedule(static) - for (int i = 0; i < nlocal; i++) - { - std::size_t count = 0; - visit_neighbors(atoms[i], binned_atoms, sradius2, - [&count](ModuleNeighList::LocalAtomIndex) { ++count; }); - neighbor_counts[static_cast(i)] = count; - } + int ix = std::min( + std::max(int((atom.position_x - x_min_) / bin_sizex_), 0), + nbinx_ - 1 + ); - for (int i = 0; i < nlocal; i++) - { - const int n = ModuleNeighList::checked_int_size( - neighbor_counts[static_cast(i)], - "BinManager neighbor count" - ); - neighbor_list.firstneigh_[i] = neighbor_list.allocator_.allocate(n); - neighbor_list.numneigh_[i] = n; - } + 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 + ); -#pragma omp parallel for schedule(static) - for (int i = 0; i < nlocal; i++) + for (int dx = -1; dx <= 1; dx++) { - int* ptr = neighbor_list.firstneigh_[i]; - int k = 0; - visit_neighbors(atoms[i], binned_atoms, sradius2, - [&](ModuleNeighList::LocalAtomIndex atom_id) - { - assert(ptr != nullptr); - ptr[k++] = atom_id; - }); - assert(k == neighbor_list.numneigh_[i]); - } - return; - } -#endif + 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; - std::vector neigh_tmp; - for (int i = 0; i < nlocal; i++) - { - neigh_tmp.clear(); - visit_neighbors(atoms[i], binned_atoms, sradius2, - [&neigh_tmp](ModuleNeighList::LocalAtomIndex atom_id) + if (jx < 0 || jx >= nbinx_ || + jy < 0 || jy >= nbiny_ || + jz < 0 || jz >= nbinz_) + continue; + + int nidx = bin_index(jx, jy, jz); + + for (const ModuleNeighList::LocalAtomIndex binned_atom_index : bins_[nidx].get_atom_indices()) + { + const NeighborAtom& natom = binned_atoms[static_cast(binned_atom_index)]; + double dx = atom.position_x - natom.position_x; + double dy = atom.position_y - natom.position_y; + double dz = atom.position_z - natom.position_z; + + double dist2 = dx * dx + dy * dy + dz * dz; + + if (natom.atom_id == atom.atom_id) { - neigh_tmp.push_back(atom_id); - }); + continue; + } + if (dist2 <= sradius2) + { + neigh_tmp.push_back(natom.atom_id); + } + } + } + } + } const int n = ModuleNeighList::checked_int_size(neigh_tmp.size(), "BinManager neighbor count"); diff --git a/source/source_cell/module_neighlist/bin_manager.h b/source/source_cell/module_neighlist/bin_manager.h index 08eae906668..ffb470e8724 100644 --- a/source/source_cell/module_neighlist/bin_manager.h +++ b/source/source_cell/module_neighlist/bin_manager.h @@ -198,21 +198,6 @@ class BinManager * @return Flat index in the bins_ array. */ int bin_index(int ix, int iy, int iz) const; - - /** - * @brief Visit neighbors of one atom in the existing deterministic bin order. - * - * @tparam Emit Callable accepting a rank-local neighbor atom ID. - * @param atom Atom used as the neighbor-list center. - * @param binned_atoms All atoms assigned to bins by do_binning(). - * @param sradius2 Squared search radius. - * @param emit Callback invoked once for every accepted neighbor. - */ - template - void visit_neighbors(const NeighborAtom& atom, - const std::vector& binned_atoms, - double sradius2, - const Emit& emit) const; }; #endif // BIN_MANAGER_H diff --git a/source/source_cell/module_neighlist/test/CMakeLists.txt b/source/source_cell/module_neighlist/test/CMakeLists.txt index 5c4f0336f9b..08d5b030ced 100644 --- a/source/source_cell/module_neighlist/test/CMakeLists.txt +++ b/source/source_cell/module_neighlist/test/CMakeLists.txt @@ -8,11 +8,10 @@ abacus_disable_feature_definitions(__EXX) AddTest( TARGET MODULE_CELL_NEIGHBOR_neighbor_search - LIBS parameter base device + LIBS parameter base device cell symmetry SOURCES neighbor_search_test.cpp ../neighbor_search.cpp - ../../md_cell.cpp ../domain_decomposition.cpp ../bin_manager.cpp ../page_allocator.cpp diff --git a/source/source_cell/module_neighlist/test/bin_manager_test.cpp b/source/source_cell/module_neighlist/test/bin_manager_test.cpp index 0d75e505ed4..274aefe2089 100644 --- a/source/source_cell/module_neighlist/test/bin_manager_test.cpp +++ b/source/source_cell/module_neighlist/test/bin_manager_test.cpp @@ -2,30 +2,6 @@ #include "source_cell/module_neighlist/bin_manager.h" #include "source_cell/module_neighlist/neighbor_list.h" -#include - -#ifdef _OPENMP -#include -#endif - -namespace -{ -std::vector> snapshot_neighbors(const NeighborList& list) -{ - std::vector> result(static_cast(list.get_nlocal())); - for (int i = 0; i < list.get_nlocal(); ++i) - { - const int count = list.get_numneigh(i); - const int* first = list.get_firstneigh(i); - if (count > 0) - { - result[static_cast(i)].assign(first, first + count); - } - } - return result; -} -} // namespace - TEST(BinManagerUnit, InitAndBinning) { std::vector inside; @@ -245,93 +221,3 @@ TEST(BinManagerUnit, MultipleBinsNeighborSearch) int center_index = 13; EXPECT_EQ(nl.get_numneigh(center_index), 6); } - -#ifdef _OPENMP -TEST(BinManagerUnit, ParallelBuildPreservesSerialNeighborOrder) -{ - std::vector centers; - std::vector binned_atoms; - centers.reserve(300); - binned_atoms.reserve(600); - for (int i = 0; i < 300; ++i) - { - const double x = 0.2 * static_cast(i % 10); - const double y = 0.2 * static_cast((i / 10) % 10); - const double z = 0.2 * static_cast(i / 100); - centers.emplace_back(x, y, z, 0, i, i); - binned_atoms.push_back(centers.back()); - } - for (int i = 0; i < 300; ++i) - { - const NeighborAtom& center = centers[static_cast(i)]; - binned_atoms.emplace_back(center.position_x + 0.05, - center.position_y, - center.position_z, - 0, - i, - 300 + i, - 1000 + i, - 1); - } - - BinManager bm; - bm.init_bins(0.31, binned_atoms); - bm.do_binning(binned_atoms); - - const int previous_dynamic = omp_get_dynamic(); - const int previous_threads = omp_get_max_threads(); - omp_set_dynamic(0); - - NeighborList serial_list; - serial_list.initialize(centers.size(), centers.size() * 128); - omp_set_num_threads(1); - bm.build_atom_neighbors(serial_list, centers, binned_atoms); - const std::vector> serial = snapshot_neighbors(serial_list); - - NeighborList parallel_list; - parallel_list.initialize(centers.size(), centers.size() * 128); - omp_set_num_threads(4); - bm.build_atom_neighbors(parallel_list, centers, binned_atoms); - const std::vector> parallel = snapshot_neighbors(parallel_list); - - omp_set_num_threads(previous_threads); - omp_set_dynamic(previous_dynamic); - - EXPECT_EQ(parallel, serial); - EXPECT_NE(std::find_if(serial[0].begin(), serial[0].end(), - [](int atom_id) { return atom_id >= 300; }), - serial[0].end()); -} - -TEST(BinManagerUnit, ParallelBuildKeepsZeroNeighborPointersNull) -{ - std::vector atoms; - atoms.reserve(256); - for (int i = 0; i < 256; ++i) - { - atoms.emplace_back(static_cast(i), 0.0, 0.0, 0, i, i); - } - - BinManager bm; - bm.init_bins(0.4, atoms); - bm.do_binning(atoms); - - const int previous_dynamic = omp_get_dynamic(); - const int previous_threads = omp_get_max_threads(); - omp_set_dynamic(0); - omp_set_num_threads(4); - - NeighborList list; - list.initialize(atoms.size(), 1024); - bm.build_atom_neighbors(list, atoms, atoms); - - omp_set_num_threads(previous_threads); - omp_set_dynamic(previous_dynamic); - - for (int i = 0; i < list.get_nlocal(); ++i) - { - EXPECT_EQ(list.get_numneigh(i), 0); - EXPECT_EQ(list.get_firstneigh(i), nullptr); - } -} -#endif diff --git a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp index ea2271e9c17..e80e9e852c1 100644 --- a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp +++ b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp @@ -1,35 +1,41 @@ #include -#include "source_cell/md_cell.h" -#include "../neighbor_search.h" -#include "source_base/communication_domain.h" +#include "source_cell/module_neighlist/local_atom.h" +#include "source_cell/module_neighlist/neighbor_search.h" +#include "source_cell/unitcell.h" -#ifdef __MPI -#include -#endif - -#include #include -#include #include -#ifdef _OPENMP -#include -#endif - namespace { -void ensure_mpi_initialized() +void initialize_test_ucell(UnitCell& ucell, + double lat0, + double omega, + const ModuleBase::Matrix3& latvec, + int ntype, + const std::vector& na, + const std::vector>& tau) { -#ifdef __MPI - int initialized = 0; - MPI_Initialized(&initialized); - if (!initialized) + ucell.lat0 = lat0; + ucell.omega = omega; + ucell.latvec = latvec; + ucell.GT = latvec.Inverse(); + ucell.ntype = ntype; + ucell.nat = 0; + ucell.atoms = new Atom[ntype]; + std::size_t iat = 0; + for (int it = 0; it < ntype; ++it) { - int provided = 0; - MPI_Init_thread(NULL, NULL, MPI_THREAD_SINGLE, &provided); + ucell.atoms[it].type = it; + ucell.atoms[it].na = na[static_cast(it)]; + ucell.atoms[it].tau.resize(static_cast(ucell.atoms[it].na)); + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + ucell.atoms[it].tau[static_cast(ia)] = tau[iat++]; + } + ucell.nat += ucell.atoms[it].na; } -#endif } ModuleBase::Matrix3 identity_lattice() @@ -47,78 +53,21 @@ ModuleBase::Matrix3 identity_lattice() return latvec; } -double cell_volume(const ModuleBase::Matrix3& latvec) -{ - const double cx = latvec.e22 * latvec.e33 - latvec.e23 * latvec.e32; - const double cy = latvec.e23 * latvec.e31 - latvec.e21 * latvec.e33; - const double cz = latvec.e21 * latvec.e32 - latvec.e22 * latvec.e31; - return std::abs(latvec.e11 * cx + latvec.e12 * cy + latvec.e13 * cz); -} - -MDCell make_mdcell(const ModuleBase::Matrix3& latvec, - const std::vector >& positions, - double cutoff) -{ - int rank = 0; -#ifdef __MPI - MPI_Comm_rank(MPI_COMM_WORLD, &rank); -#endif - const ModuleBase::Matrix3 gt = latvec.Inverse(); - std::vector owned_atoms; - owned_atoms.reserve(positions.size()); - for (std::size_t iat = 0; iat < positions.size(); ++iat) - { - owned_atoms.push_back(LocalAtom(positions[iat], positions[iat] * gt, - ModuleBase::Vector3(0.0, 0.0, 0.0), - ModuleBase::Vector3(0.0, 0.0, 0.0), - ModuleBase::Vector3(1, 1, 1), 1.0, 0, - static_cast(iat), rank, false)); - } - return MDCell(latvec, gt, 1.0, cell_volume(latvec), static_cast(positions.size()), - owned_atoms, std::vector(1, "X"), - std::vector(1, 1.0), cutoff, 0.0, - ModuleBase::world_communication_domain()); -} - -std::size_t count_pairs(const NeighborList& list) -{ - std::size_t pairs = 0; - for (int local_i = 0; local_i < list.get_nlocal(); ++local_i) - { - pairs += static_cast(list.get_numneigh(local_i)); - } - return pairs; -} - -#ifdef _OPENMP -void expect_same_atoms(const std::vector& lhs, - const std::vector& rhs) -{ - ASSERT_EQ(lhs.size(), rhs.size()); - for (std::size_t i = 0; i < lhs.size(); ++i) - { - EXPECT_DOUBLE_EQ(lhs[i].position_x, rhs[i].position_x) << "atom " << i; - EXPECT_DOUBLE_EQ(lhs[i].position_y, rhs[i].position_y) << "atom " << i; - EXPECT_DOUBLE_EQ(lhs[i].position_z, rhs[i].position_z) << "atom " << i; - EXPECT_EQ(lhs[i].atom_type, rhs[i].atom_type) << "atom " << i; - EXPECT_EQ(lhs[i].atom_index, rhs[i].atom_index) << "atom " << i; - EXPECT_EQ(lhs[i].atom_id, rhs[i].atom_id) << "atom " << i; - EXPECT_EQ(lhs[i].owner_rank, rhs[i].owner_rank) << "atom " << i; - } -} -#endif } // namespace -TEST(NeighborSearchTest, MdCellTwoAtomsNeighbor) +TEST(NeighborSearchTest, TwoAtomsNeighbor) { - ensure_mpi_initialized(); - - const ModuleBase::Matrix3 latvec = identity_lattice(); - const std::vector > positions{{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}; - MDCell mdcell = make_mdcell(latvec, positions, 1.0); + UnitCell ucell; + initialize_test_ucell(ucell, + 1.0, + 1.0, + identity_lattice(), + 1, + {2}, + {{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}); NeighborSearch ns; - ns.init(mdcell, 1.0); + ns.init(ucell, 1.0); ns.build_neighbors(); const NeighborList& list = ns.get_neighbor_list(); @@ -127,16 +76,19 @@ TEST(NeighborSearchTest, MdCellTwoAtomsNeighbor) EXPECT_EQ(list.get_numneigh(1), 8); } -TEST(NeighborSearchTest, MdCellNoNeighbor) +TEST(NeighborSearchTest, NoNeighbor) { - ensure_mpi_initialized(); - - const ModuleBase::Matrix3 latvec = identity_lattice(); - const std::vector > positions{{0.0, 0.0, 0.0}, {0.49, 0.0, 0.0}}; - MDCell mdcell = make_mdcell(latvec, positions, 0.1); + UnitCell ucell; + initialize_test_ucell(ucell, + 1.0, + 1.0, + identity_lattice(), + 1, + {2}, + {{0.0, 0.0, 0.0}, {5.0, 0.0, 0.0}}); NeighborSearch ns; - ns.init(mdcell, 0.1); + ns.init(ucell, 0.1); ns.build_neighbors(); const NeighborList& list = ns.get_neighbor_list(); @@ -145,132 +97,28 @@ TEST(NeighborSearchTest, MdCellNoNeighbor) EXPECT_EQ(list.get_numneigh(1), 0); } -TEST(NeighborSearchTest, MdCellInitBuildsOwnedAndGhostAtoms) +TEST(NeighborSearchTest, SerialInitOwnsCentralAtomsAndBuildsImages) { - ensure_mpi_initialized(); - - const ModuleBase::Matrix3 latvec = identity_lattice(); - const std::vector > positions{{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}; - MDCell mdcell = make_mdcell(latvec, positions, 1.0); + UnitCell ucell; + initialize_test_ucell(ucell, + 1.0, + 1.0, + identity_lattice(), + 1, + {2}, + {{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}); NeighborSearch ns; - ns.init(mdcell, 1.0); + ns.init(ucell, 1.0); EXPECT_EQ(ns.get_inside_atoms().size(), 2U); - EXPECT_GT(ns.get_ghost_atoms().size(), 0U); - EXPECT_GT(ns.get_all_atoms().size(), ns.get_inside_atoms().size()); EXPECT_EQ(ns.get_neighbor_list().get_nlocal(), 2); -} + EXPECT_EQ(ns.get_all_atoms().size(), 54U); -TEST(NeighborSearchTest, MdCellNeighborIdsStayLocalToAllAtoms) -{ - ensure_mpi_initialized(); - - const ModuleBase::Matrix3 latvec = identity_lattice(); - const std::vector > positions{{0.0, 0.0, 0.0}, - {0.5, 0.0, 0.0}, - {0.0, 0.5, 0.0}}; - MDCell mdcell = make_mdcell(latvec, positions, 0.75); - - NeighborSearch ns; - ns.init(mdcell, 0.75); - ns.build_neighbors(); - - const NeighborList& list = ns.get_neighbor_list(); const std::vector& all_atoms = ns.get_all_atoms(); - EXPECT_GT(count_pairs(list), 0U); - for (std::size_t atom_id = 0; atom_id < all_atoms.size(); ++atom_id) - { - EXPECT_EQ(all_atoms[atom_id].atom_id, - ModuleNeighList::checked_local_atom_index(atom_id, "test atom id")); - } - for (int local_i = 0; local_i < list.get_nlocal(); ++local_i) - { - for (int ad = 0; ad < list.get_numneigh(local_i); ++ad) - { - const int neighbor_id = list.get_firstneigh(local_i)[ad]; - EXPECT_GE(neighbor_id, 0); - EXPECT_LT(static_cast(neighbor_id), all_atoms.size()); - } - } -} - -TEST(NeighborSearchTest, MdCellPreservesMdAtomStateAcrossOwnedAndGhostStorage) -{ - ensure_mpi_initialized(); - - const ModuleBase::Matrix3 latvec = identity_lattice(); - const std::vector > positions{{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}; - MDCell mdcell = make_mdcell(latvec, positions, 1.0); - - ASSERT_EQ(mdcell.nlocal(), 2); - std::vector& owned_atoms = mdcell.mutable_owned_atoms(); - owned_atoms[0].vel.set(1.0, 2.0, 3.0); - owned_atoms[0].mbl.set(1, 0, 1); - owned_atoms[0].mass = 7.5; - owned_atoms[1].vel.set(-1.0, -2.0, -3.0); - owned_atoms[1].mbl.set(0, 1, 1); - owned_atoms[1].mass = 8.5; - - mdcell.exchange_ghost_atoms(); - - ASSERT_GT(mdcell.nghost(), 0); - const std::vector& ghost_atoms = mdcell.ghost_atoms(); - EXPECT_DOUBLE_EQ(ghost_atoms[0].force.x, 0.0); - EXPECT_DOUBLE_EQ(ghost_atoms[0].force.y, 0.0); - EXPECT_DOUBLE_EQ(ghost_atoms[0].force.z, 0.0); - - bool found_first = false; - bool found_second = false; - for (std::size_t i = 0; i < ghost_atoms.size(); ++i) + for (std::size_t i = 0; i < all_atoms.size(); ++i) { - if (ghost_atoms[i].type == owned_atoms[0].type && - ghost_atoms[i].type_index == owned_atoms[0].type_index) - { - found_first = true; - EXPECT_DOUBLE_EQ(ghost_atoms[i].vel.x, 1.0); - EXPECT_DOUBLE_EQ(ghost_atoms[i].vel.y, 2.0); - EXPECT_DOUBLE_EQ(ghost_atoms[i].vel.z, 3.0); - EXPECT_EQ(ghost_atoms[i].mbl.x, 1); - EXPECT_EQ(ghost_atoms[i].mbl.y, 0); - EXPECT_EQ(ghost_atoms[i].mbl.z, 1); - EXPECT_DOUBLE_EQ(ghost_atoms[i].mass, 7.5); - } - if (ghost_atoms[i].type == owned_atoms[1].type && - ghost_atoms[i].type_index == owned_atoms[1].type_index) - { - found_second = true; - EXPECT_DOUBLE_EQ(ghost_atoms[i].vel.x, -1.0); - EXPECT_DOUBLE_EQ(ghost_atoms[i].vel.y, -2.0); - EXPECT_DOUBLE_EQ(ghost_atoms[i].vel.z, -3.0); - EXPECT_EQ(ghost_atoms[i].mbl.x, 0); - EXPECT_EQ(ghost_atoms[i].mbl.y, 1); - EXPECT_EQ(ghost_atoms[i].mbl.z, 1); - EXPECT_DOUBLE_EQ(ghost_atoms[i].mass, 8.5); - } + EXPECT_EQ(all_atoms[i].atom_id, + ModuleNeighList::checked_local_atom_index(i, "test atom id")); } - - EXPECT_TRUE(found_first); - EXPECT_TRUE(found_second); -} - -TEST(NeighborSearchTest, MdCellMigrateOwnedAtomsReassignsOwnership) -{ - ensure_mpi_initialized(); - - const ModuleBase::Matrix3 latvec = identity_lattice(); - const std::vector > positions{{0.1, 0.1, 0.1}}; - MDCell mdcell = make_mdcell(latvec, positions, 0.2); - - ASSERT_EQ(mdcell.nlocal(), 1); - mdcell.mutable_owned_atoms()[0].cart.set(1.2, -0.1, 0.1); - mdcell.migrate_owned_atoms(); - - ASSERT_EQ(mdcell.nlocal(), 1); - EXPECT_DOUBLE_EQ(mdcell.owned_atoms()[0].frac.x, 0.2); - EXPECT_DOUBLE_EQ(mdcell.owned_atoms()[0].frac.y, 0.9); - EXPECT_DOUBLE_EQ(mdcell.owned_atoms()[0].frac.z, 0.1); - EXPECT_DOUBLE_EQ(mdcell.owned_atoms()[0].cart.x, 0.2); - EXPECT_DOUBLE_EQ(mdcell.owned_atoms()[0].cart.y, 0.9); - EXPECT_DOUBLE_EQ(mdcell.owned_atoms()[0].cart.z, 0.1); }