diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 6aeaa0c75a..d62f2d325f 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -428,7 +428,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 ca64ccf60d..3aa11e61df 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 0000000000..314ac05eea --- /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 0000000000..1af6fa6670 --- /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 0000000000..6ae4fd138c --- /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 0000000000..27fe453a22 --- /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 dc3e1e7c50..3e6c282f02 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 3087148dac..0000000000 --- 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 abbf529a68..443a7d6313 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 9b74729abd..3b00c3d56f 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 f48a8da8f7..5969a8508e 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 3f62d30571..5e805d0e21 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 74e21cfac6..cd5c5125ab 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 b95ace5bc6..fae73246d8 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 a3a95aeb31..739f72445e 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 281b4fc5fc..ae0b535c9f 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/bin_manager_test.cpp b/source/source_cell/module_neighlist/test/bin_manager_test.cpp index 07e34c488a..274aefe208 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; 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 0000000000..7691e9d19e --- /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 0000000000..406b88e02e --- /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 0837d38ebe..0000000000 --- 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 c689019691..761c37f46e 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 877d214497..0000000000 --- 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 e951945799..0000000000 --- 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 2be50a9877..bc992f1628 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 30b8ff78f5..3b499d1add 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 42ed6cfcc7..fd0b390e1e 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 9f74a519a7..cb72982a3f 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