diff --git a/CMakeLists.txt b/CMakeLists.txt index 204fc5b1775..1de1439f365 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -580,6 +580,11 @@ if(ESPRESSO_BUILD_WITH_SHARED_MEMORY_PARALLELISM) set(CMAKE_SHARED_LIBRARY_PREFIX "lib") set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "") set(Kokkos_ENABLE_OPENMP ON CACHE BOOL "") + set(Kokkos_ENABLE_IMPL_VIEW_LEGACY ON CACHE BOOL "") + set(Kokkos_ENABLE_COMPLEX_ALIGN ON CACHE BOOL "") + set(Kokkos_ENABLE_AGGRESSIVE_VECTORIZATION ON CACHE BOOL "") + set(Kokkos_ENABLE_HWLOC ON CACHE BOOL "") + set(Kokkos_ARCH_NATIVE ON CACHE BOOL "") FetchContent_MakeAvailable(kokkos) set(BUILD_SHARED_LIBS ${ESPRESSO_BUILD_SHARED_LIBS_DEFAULT}) set(CMAKE_SHARED_LIBRARY_PREFIX "${ESPRESSO_SHARED_LIBRARY_PREFIX}") @@ -604,7 +609,7 @@ if(ESPRESSO_BUILD_WITH_SHARED_MEMORY_PARALLELISM) FetchContent_Declare( cabana GIT_REPOSITORY https://github.com/ECP-copa/Cabana.git - GIT_TAG ebfaa51 # 0.7.0 with patches + GIT_TAG e76c1a1 # 0.7.0 with patches PATCH_COMMAND patch -p0 < ${CMAKE_CURRENT_SOURCE_DIR}/cmake/cabana.patch ) # cmake-format: on @@ -628,6 +633,7 @@ if(ESPRESSO_BUILD_WITH_SHARED_MEMORY_PARALLELISM) list(APPEND SKIP_CLANG_TIDY_CHECKS "-modernize-use-bool-literals") list(APPEND SKIP_CLANG_TIDY_CHECKS "-modernize-use-equals-delete") list(APPEND SKIP_CLANG_TIDY_CHECKS "-modernize-use-equals-default") + list(APPEND SKIP_CLANG_TIDY_CHECKS "-modernize-use-nullptr") list(APPEND SKIP_CLANG_TIDY_CHECKS "-modernize-pass-by-value") list(APPEND SKIP_CLANG_TIDY_CHECKS "-modernize-loop-convert") list(APPEND SKIP_CLANG_TIDY_CHECKS "-modernize-return-braced-init-list") @@ -636,6 +642,10 @@ if(ESPRESSO_BUILD_WITH_SHARED_MEMORY_PARALLELISM) list(APPEND SKIP_CLANG_TIDY_CHECKS "-readability-simplify-boolean-expr") list(APPEND SKIP_CLANG_TIDY_CHECKS "-readability-avoid-const-params-in-decls") + list(APPEND SKIP_CLANG_TIDY_CHECKS + "-clang-analyzer-optin.performance.Padding") + list(APPEND SKIP_CLANG_TIDY_CHECKS + "-clang-analyzer-security.insecureAPI.strcpy") espresso_override_clang_tidy_checks( ESPRESSO_CABANA_CXX_CLANG_TIDY "${SKIP_CLANG_TIDY_CHECKS}" "${SKIP_CLANG_TIDY_CHECKS_CXX}") diff --git a/src/core/BoxGeometry.hpp b/src/core/BoxGeometry.hpp index b16fb3c48d8..fdddb16161d 100644 --- a/src/core/BoxGeometry.hpp +++ b/src/core/BoxGeometry.hpp @@ -225,6 +225,40 @@ class BoxGeometry { get_mi_coord(a[2], b[2], 2)}; } + /** + * @brief Get the minimum-image vector between two coordinates. + * + * @tparam T Floating point type. + * + * @param a0 x element of the terminal point. + * @param a1 y element of the terminal point. + * @param a2 z element of the terminal point. + * @param b0 x element of the initial point. + * @param b1 y element of the initial point. + * @param b2 z element of the initial point. + * @return Vector from @p b to @p a that minimizes the distance across + * periodic images, i.e. a - b. + */ + template + Utils::Vector get_mi_vector(T const &a0, T const &a1, T const &a2, + T const &b0, T const &b1, + T const &b2) const { + if (type() == BoxType::LEES_EDWARDS) { + auto const shear_plane_normal = lees_edwards_bc().shear_plane_normal; + auto a_tmp = Utils::Vector{a0, a1, a2}; + auto b_tmp = Utils::Vector{b0, b1, b2}; + a_tmp[shear_plane_normal] = Algorithm::periodic_fold( + a_tmp[shear_plane_normal], m_length[shear_plane_normal]); + b_tmp[shear_plane_normal] = Algorithm::periodic_fold( + b_tmp[shear_plane_normal], m_length[shear_plane_normal]); + return lees_edwards_bc().distance(a_tmp - b_tmp, m_length, m_length_half, + m_length_inv, m_periodic); + } + assert(type() == BoxType::CUBOID); + return {get_mi_coord(a0, b0, 0), get_mi_coord(a1, b1, 1), + get_mi_coord(a2, b2, 2)}; + } + BoxType type() const { return m_type; } void set_type(BoxType type) { m_type = type; } diff --git a/src/core/aosoa_pack.hpp b/src/core/aosoa_pack.hpp new file mode 100644 index 00000000000..3215fb01f88 --- /dev/null +++ b/src/core/aosoa_pack.hpp @@ -0,0 +1,41 @@ +/* + * Copyright (C) 2025 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#ifdef SHARED_MEMORY_PARALLELISM + +#include "cell_system/CellStructure.hpp" + +#include + +struct CellStructure::AoSoA_pack { + CellStructure::AoSoAType::member_slice_type<0> position; + CellStructure::AoSoAType::member_slice_type<1> charge; + CellStructure::AoSoAType::member_slice_type<2> id; + CellStructure::AoSoAType::member_slice_type<3> type; + + AoSoA_pack() = default; + + AoSoA_pack(CellStructure::AoSoAType &aosoa) + : position(Cabana::slice<0>(aosoa)), charge(Cabana::slice<1>(aosoa)), + id(Cabana::slice<2>(aosoa)), type(Cabana::slice<3>(aosoa)) {} +}; + +#endif // SHARED_MEMORY_PARALLELISM diff --git a/src/core/cell_system/CellStructure.cpp b/src/core/cell_system/CellStructure.cpp index 5b7b67842a9..2e477079ba3 100644 --- a/src/core/cell_system/CellStructure.cpp +++ b/src/core/cell_system/CellStructure.cpp @@ -32,11 +32,13 @@ #include "cell_system/CellStructureType.hpp" #include "communication.hpp" #include "lees_edwards/lees_edwards.hpp" +#include "particle_enumeration.hpp" #include "particle_reduction.hpp" #include "system/System.hpp" #include #include +#include #include #include @@ -56,9 +58,139 @@ #include #ifdef SHARED_MEMORY_PARALLELISM +#include "aosoa_pack.hpp" +#include "custom_verlet_list.hpp" +#include +#include #include #endif +CellStructure::~CellStructure() { +#ifdef SHARED_MEMORY_PARALLELISM + if (m_local_force) { + m_local_force.reset(); + } +#ifdef ROTATION + if (m_local_torque) { + m_local_torque.reset(); + } +#endif +#ifdef NPT + if (m_local_virial) { + m_local_virial.reset(); + } +#endif + if (m_aosoa) { + m_aosoa.reset(); + } + if (m_particle_storage) { + m_particle_storage.reset(); + } + if (m_verlet_list_cabana) { + m_verlet_list_cabana.reset(); + } + // Kokkos handle can be freed after all Cabana containers have been freed + m_kokkos_handle.reset(); +#endif +} + +#ifdef SHARED_MEMORY_PARALLELISM + +void CellStructure::set_kokkos_handle( + std::shared_ptr handle) { + m_kokkos_handle = std::move(handle); +} + +static auto estimate_max_counts(int max_prefactor, double pair_cutoff, + std::size_t number_of_unique_particles) { + if (std::isinf(pair_cutoff)) { + return number_of_unique_particles; + } + if (pair_cutoff < 0.) { + pair_cutoff = 0.; + } + auto const volume = Utils::int_pow<3>(pair_cutoff); + auto max_counts = static_cast( + std::ceil(static_cast(max_prefactor) * volume)); +#ifdef COLLISION_DETECTION + std::size_t constexpr threshold_num = 64; +#else + std::size_t constexpr threshold_num = 16; +#endif + if (max_counts < threshold_num) { + max_counts = std::min(threshold_num, number_of_unique_particles); + } + return max_counts; +} + +void CellStructure::rebuild_local_properties(std::size_t const num_threads, + double const pair_cutoff) { + assert(m_kokkos_handle); + auto const num_part = get_unique_particles().size(); + m_local_force = + std::make_unique("local_force", num_part, num_threads); +#ifdef ROTATION + m_local_torque = + std::make_unique("local_torque", num_part, num_threads); +#endif +#ifdef NPT + m_local_virial = std::make_unique("local_virial", num_threads); +#endif + m_particle_storage = std::make_unique("particles", num_part); + m_particle_storage->resize(num_part); + // particle properties are defined in aosoa_pack.hpp + m_aosoa = std::make_unique(*m_particle_storage); + + auto max_counts = estimate_max_counts(m_max_prefactor, pair_cutoff, num_part); + m_verlet_list_cabana = std::make_unique(0ul, num_part, max_counts); +} + +void CellStructure::reset_local_properties() { + Kokkos::deep_copy(get_local_force(), 0); +#ifdef ROTATION + Kokkos::deep_copy(get_local_torque(), 0); +#endif +#ifdef NPT + Kokkos::deep_copy(get_local_virial(), 0); +#endif +} + +void CellStructure::set_index_map() { + auto &unique_particles = m_unique_particles; + unique_particles.clear(); + unique_particles.resize(count_local_particles()); + std::unordered_set registered_index{}; + using execution_space = Kokkos::DefaultExecutionSpace; + int n_threads = execution_space().concurrency(); + std::vector max_ids(n_threads); + enumerate_local_particles( + *this, [&unique_particles, &max_ids](std::size_t index, Particle &p) { + unique_particles[index] = &p; + const int thread_num = omp_get_thread_num(); + max_ids[thread_num] = std::max(p.id(), max_ids[thread_num]); + }); + int max_id = *(std::max_element(max_ids.begin(), max_ids.end())); + for (auto &p : ghost_particles()) { + auto const *local_particle = get_local_particle(p.id()); + if (not local_particle) { + continue; + } + if (not local_particle->is_ghost()) { + continue; + } + if (registered_index.contains(p.id())) { + continue; + } + registered_index.insert(p.id()); + unique_particles.emplace_back(&p); + max_id = std::max(p.id(), max_id); + } + registered_index.clear(); + m_cached_max_local_particle_id = max_id; +} + +#endif // SHARED_MEMORY_PARALLELISM + CellStructure::CellStructure(BoxGeometry const &box) : m_decomposition{std::make_unique(box)} {} @@ -236,6 +368,7 @@ void CellStructure::resort_particles(bool global_flag) { auto const &lebc = get_system().box_geo->lees_edwards_bc(); m_rebuild_verlet_list = true; + m_rebuild_verlet_list_cabana = true; m_le_pos_offset_at_last_resort = lebc.pos_offset; #ifdef ADDITIONAL_CHECKS @@ -285,6 +418,7 @@ void CellStructure::set_verlet_skin(double value) { assert(value >= 0.); m_verlet_skin = value; m_verlet_skin_set = true; + m_rebuild_verlet_list_cabana = true; get_system().on_verlet_skin_change(); } diff --git a/src/core/cell_system/CellStructure.hpp b/src/core/cell_system/CellStructure.hpp index 31fe56822c7..32b8eaa5ae9 100644 --- a/src/core/cell_system/CellStructure.hpp +++ b/src/core/cell_system/CellStructure.hpp @@ -52,9 +52,42 @@ #include #include #include +#include #include #include +#ifdef CALIPER +#include +#endif + +// forward declarations +#ifdef SHARED_MEMORY_PARALLELISM +namespace Kokkos { +template class View; +class HostSpace; +struct LayoutRight; +template struct MemoryTraits; +} // namespace Kokkos +namespace Cabana { +class HalfNeighborTag; +struct VerletLayout2D; +class TeamVectorOpTag; +template struct MemberTypes; +template +class AoSoA; +} // namespace Cabana +namespace Communication { +struct KokkosHandle; +} // namespace Communication +template +class CustomVerletList; +#endif // SHARED_MEMORY_PARALLELISM + +template +concept ParticleCallback = requires(Callable c, Particle &p) { + { c(p) } -> std::same_as; +}; + using ParticleUnaryOp = std::function; namespace Cells { @@ -135,7 +168,23 @@ struct EuclidianDistance { * system which are not common between different cell systems have to * be stored in separate structures. */ -struct CellStructure : public System::Leaf { +class CellStructure : public System::Leaf { +#ifdef SHARED_MEMORY_PARALLELISM +public: + static constexpr auto vector_length = 1; + struct AoSoA_pack; + using ForceType = Kokkos::View; + using VirialType = Kokkos::View; + using data_types = Cabana::MemberTypes; + using memory_space = Kokkos::HostSpace; + using AoSoAType = Cabana::AoSoA>; + using ListAlgorithm = Cabana::HalfNeighborTag; + using ListType = + CustomVerletList; +#endif // SHARED_MEMORY_PARALLELISM + private: /** The local id-to-particle index */ std::vector m_particle_index; @@ -146,16 +195,37 @@ struct CellStructure : public System::Leaf { /** One of @ref Cells::Resort, announces the level of resort needed. */ unsigned m_resort_particles = Cells::RESORT_NONE; + bool m_verlet_skin_set = false; bool m_rebuild_verlet_list = true; + bool m_rebuild_verlet_list_cabana = true; std::vector> m_verlet_list; double m_le_pos_offset_at_last_resort = 0.; /** @brief Verlet list skin. */ double m_verlet_skin = 0.; - bool m_verlet_skin_set = false; double m_verlet_reuse = 0.; +#ifdef SHARED_MEMORY_PARALLELISM + int m_cached_max_local_particle_id = 0; + int m_max_prefactor = 8; + int m_max_id = 0; + std::unique_ptr m_local_force; +#ifdef ROTATION + std::unique_ptr m_local_torque; +#endif +#ifdef NPT + std::unique_ptr m_local_virial; +#endif + std::unique_ptr m_verlet_list_cabana; + std::unique_ptr m_particle_storage; + /** particle properties for Cabana */ + std::unique_ptr m_aosoa; + /** The local id-to-index for aosoa data */ + std::vector m_unique_particles; + std::shared_ptr m_kokkos_handle; +#endif // SHARED_MEMORY_PARALLELISM public: CellStructure(BoxGeometry const &box); + virtual ~CellStructure(); bool use_verlet_list = true; @@ -278,6 +348,15 @@ struct CellStructure : public System::Leaf { ParticleRange ghost_particles() const { return Cells::particles(decomposition().ghost_cells()); } + + std::size_t count_local_particles() const { + std::size_t count = 0; + for (auto const &cell : m_decomposition->local_cells()) { + count += cell->particles().size(); + } + return count; + } + /** @brief whether to use parallel version of @ref for_each_local_particle */ bool use_parallel_for_each_local_particle() const { #ifdef SHARED_MEMORY_PARALLELISM @@ -380,6 +459,11 @@ struct CellStructure : public System::Leaf { * this node, or -1 if there are no particles on this node. */ int get_max_local_particle_id() const; +#ifdef SHARED_MEMORY_PARALLELISM + int get_cached_max_local_particle_id() const { + return m_cached_max_local_particle_id; + } +#endif /** * @brief Remove all particles from the cell system. @@ -639,6 +723,89 @@ struct CellStructure : public System::Leaf { } } +#ifdef SHARED_MEMORY_PARALLELISM +public: + void set_max_prefactor(int value) { m_max_prefactor = value; } + auto get_max_id() const { return m_max_id; } + + void set_kokkos_handle(std::shared_ptr handle); + void rebuild_local_properties(std::size_t num_threads, double pair_cutoff); + void reset_local_properties(); + + auto &get_local_force() { return *m_local_force; } +#ifdef ROTATION + auto &get_local_torque() { return *m_local_torque; } +#endif +#ifdef NPT + auto &get_local_virial() { return *m_local_virial; } +#endif + auto &get_aosoa() { return *m_aosoa; } + auto const &get_unique_particles() const { return m_unique_particles; } + auto const &get_verlet_list_cabana() const { return *m_verlet_list_cabana; } + + [[nodiscard]] auto is_verlet_list_cabana_rebuild_needed() const { + return m_rebuild_verlet_list_cabana or (not use_verlet_list); + } + + /** + * @brief Reset local properties of the Verlet list. + * @param n_threads Number of threads. + * @param cutoff Pair interaction cutoff. + * @return True if a rebuild is needed. + */ + [[nodiscard]] auto prepare_verlet_list_cabana(int n_threads, double cutoff) { + auto const rebuild = is_verlet_list_cabana_rebuild_needed(); + if (rebuild) { + // If we have to rebuild, we need to count the particles + set_index_map(); // parallelized index_map + // Create essential variables for MD + rebuild_local_properties(n_threads, cutoff); + } else { + // If we do not rebuild we can use the saved map + reset_local_properties(); + } + return rebuild; + } + + void rebuild_verlet_list_cabana(auto &&kernel) { + assert(is_verlet_list_cabana_rebuild_needed()); + kernel(m_decomposition->local_cells(), m_decomposition->box(), + *m_verlet_list_cabana); + m_rebuild_verlet_list_cabana = false; + } + + void set_index_map(); + inline void set_index_map(ParticleRange const &particles, + ParticleRange const &ghost_particles) { + m_unique_particles.clear(); + m_max_id = 0; + std::unordered_set registered_index{}; + for (auto &p : particles) { + if (p.id() > m_max_id) + m_max_id = p.id(); + m_unique_particles.emplace_back(&p); + } + + for (auto &p : ghost_particles) { + if (not get_local_particle(p.id())) { + continue; + } + if (not get_local_particle(p.id())->is_ghost()) { + continue; + } + if (registered_index.contains(p.id())) { + continue; + } + if (p.id() > m_max_id) + m_max_id = p.id(); + registered_index.insert(p.id()); + m_unique_particles.emplace_back(&p); + } + registered_index.clear(); + } +#endif + +private: /** Non-bonded pair loop with verlet lists. * * @param pair_kernel Kernel to apply @@ -661,6 +828,7 @@ struct CellStructure : public System::Leaf { }); m_rebuild_verlet_list = false; + m_rebuild_verlet_list_cabana = true; } else { auto const maybe_box = decomposition().minimum_image_distance(); /* In this case the pair kernel is just run over the verlet list. */ diff --git a/src/core/cell_system/particle_enumeration.hpp b/src/core/cell_system/particle_enumeration.hpp new file mode 100644 index 00000000000..fa9e2adc575 --- /dev/null +++ b/src/core/cell_system/particle_enumeration.hpp @@ -0,0 +1,82 @@ +/* + * Copyright (C) 2025 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#include "Cell.hpp" +#include "CellStructure.hpp" +#include "config/config.hpp" + +#ifdef SHARED_MEMORY_PARALLELISM +#include +#endif + +#include +#include +#include + +/** + * @brief Run a kernel on all local particles with enumeration. + * The kernel is called with (index, particle) and is assumed to be thread-safe. + * + * @tparam Kernel Callable with signature void(int, Particle&) + * @param cs The CellStructure containing the particles + * @param kernel The kernel to apply to each particle with its index + */ +template +void enumerate_local_particles(CellStructure const &cs, Kernel &&kernel); + +template +inline void enumerate_local_particles(CellStructure const &cs, + Kernel &&kernel) { +#ifdef SHARED_MEMORY_PARALLELISM + if (cs.use_parallel_for_each_local_particle()) { + auto const local_cells = cs.decomposition().local_cells(); + + // Step 1: Calculate cell offsets + std::vector cell_offsets(local_cells.size() + 1, 0); + + // Calculate cumulative sum of particles per cell + for (std::size_t i = 0; i < local_cells.size(); ++i) { + cell_offsets[i + 1] = + cell_offsets[i] + local_cells[i]->particles().size(); + } + + // Step 2: Parallel loop over cells + Kokkos::parallel_for( + "enumerate_local_particles", local_cells.size(), [&](auto cell_idx) { + auto const base_offset = cell_offsets[cell_idx]; + auto &cell_particles = local_cells[cell_idx]->particles(); + + // Loop over particles in this cell + for (std::size_t part_idx = 0; part_idx < cell_particles.size(); + ++part_idx) { + int global_index = base_offset + part_idx; + kernel(global_index, *(cell_particles.begin() + part_idx)); + } + }); + return; + } +#endif // SHARED_MEMORY_PARALLELISM + // Sequential fallback + std::size_t index = 0; + for (auto &p : cs.local_particles()) { + kernel(index++, p); + } +} diff --git a/src/core/communication.cpp b/src/core/communication.cpp index e3181ce18ce..4fed2019173 100644 --- a/src/core/communication.cpp +++ b/src/core/communication.cpp @@ -53,8 +53,20 @@ #include #include +#ifdef SHARED_MEMORY_PARALLELISM +namespace Communication { +struct KokkosHandle { + KokkosHandle() { Kokkos::initialize(); } + ~KokkosHandle() { Kokkos::finalize(); } +}; +} // namespace Communication +#endif + boost::mpi::communicator comm_cart; Communicator communicator{}; +#ifdef SHARED_MEMORY_PARALLELISM +std::shared_ptr kokkos_handle{}; +#endif namespace Communication { static std::shared_ptr m_callbacks; @@ -106,7 +118,7 @@ void init(std::shared_ptr mpi_env) { #endif #ifdef SHARED_MEMORY_PARALLELISM - Kokkos::initialize(); + kokkos_handle = std::make_shared(); #endif } @@ -115,7 +127,7 @@ void deinit() { Communication::m_callbacks.reset(); #ifdef SHARED_MEMORY_PARALLELISM - Kokkos::finalize(); + kokkos_handle.reset(); #endif } } // namespace Communication diff --git a/src/core/communication.hpp b/src/core/communication.hpp index 905074ed002..bd2dbb21477 100644 --- a/src/core/communication.hpp +++ b/src/core/communication.hpp @@ -18,8 +18,9 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef CORE_COMMUNICATION_HPP -#define CORE_COMMUNICATION_HPP + +#pragma once + /** \file * This file contains the asynchronous MPI communication. * @@ -58,6 +59,12 @@ extern int this_node; /** The communicator */ extern boost::mpi::communicator comm_cart; +#ifdef SHARED_MEMORY_PARALLELISM +namespace Communication { +struct KokkosHandle; +} // namespace Communication +extern std::shared_ptr kokkos_handle; +#endif struct Communicator { boost::mpi::communicator &comm; @@ -108,5 +115,3 @@ namespace Communication { void init(std::shared_ptr mpi_env); void deinit(); } // namespace Communication - -#endif diff --git a/src/core/custom_verlet_list.hpp b/src/core/custom_verlet_list.hpp new file mode 100644 index 00000000000..d68e50e7909 --- /dev/null +++ b/src/core/custom_verlet_list.hpp @@ -0,0 +1,193 @@ +/* + * Copyright (C) 2025 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#ifdef SHARED_MEMORY_PARALLELISM + +#include + +#include +#include +#include + +// ONLY FOR 2D LAYOUT, OTHERWISE NEIGHBOR LIST INTERFACE IMPLEMENTATION WILL +// CAUSE PROBLEMS (NOT IMPLEMENTED) +template +class CustomVerletList : public Cabana::VerletList { +public: + CustomVerletList() = default; + CustomVerletList(std::size_t const begin, std::size_t const end, + std::size_t const max_neigh) { + initializeData(end - begin, max_neigh); + } + + Kokkos::View counts; + Kokkos::View neighbors; + + // Method to initialize _data without filling neighbors + KOKKOS_INLINE_FUNCTION + void initializeData(std::size_t const num_particles, + std::size_t const max_neigh) { + counts = Kokkos::View("num_neighbors", num_particles); + neighbors = Kokkos::View( + Kokkos::ViewAllocateWithoutInitializing("neighbors"), num_particles, + max_neigh); + } + + // Method to add a neighbor + KOKKOS_INLINE_FUNCTION + void addNeighborAtomicLB(int pid, int nid) { + auto count = counts(pid); + auto count_n = counts(nid); + + if (count > count_n) { + std::swap(pid, nid); + } + count = Kokkos::atomic_fetch_add(&counts(pid), 1); + assert(count < neighbors.extent(1)); + neighbors(pid, count) = nid; + } + + // Thread safe but non atomic method to add a neighbor + KOKKOS_INLINE_FUNCTION + void addNeighbor(int pid, int nid) { + auto const count = counts(pid); + + assert(count < neighbors.extent(1)); + neighbors(pid, count) = nid; + counts(pid) += 1; + } + + // Non atomic and load balancing method to add a neighbor + KOKKOS_INLINE_FUNCTION + void addNeighborLB(int pid, int nid) { + auto count = counts(pid); + auto count_n = counts(nid); + + if (count > count_n) { + std::swap(pid, nid); + count = counts(pid); + } + assert(count < neighbors.extent(1)); + neighbors(pid, count) = nid; + counts(pid) += 1; + } + + // Sorting a neighbor + KOKKOS_INLINE_FUNCTION + void sortNeighbors() { + Kokkos::parallel_for( + "custom_verlet_list::sort_neighbors", + Kokkos::RangePolicy(0, counts.size()), + [&](const int i) { + const int count = counts(i); + auto *ptr = &neighbors(i, 0); + std::sort(ptr, ptr + count); + }); + Kokkos::fence(); + } + + // Find max counts + KOKKOS_INLINE_FUNCTION + auto get_variance_max_counts(auto &ostream) { + auto max_counts = 0l; + auto ave_counts = 0l; + auto ave_sq_counts = 0l; + for (int pid = 0; pid < counts.extent(0); ++pid) { + auto const count = static_cast(counts(pid)); + if (max_counts < count) + max_counts = count; + ave_counts += count; + ave_sq_counts += count * count; + } + if (counts.extent(0) != 0) { + ave_counts /= static_cast(counts.extent(0)); + ave_sq_counts /= static_cast(counts.extent(0)); + ave_sq_counts -= ave_counts * ave_counts; + } + ostream << "max:" << max_counts << " ave:" << ave_counts + << " var:" << ave_sq_counts << std::endl; + return static_cast(max_counts); + } + + KOKKOS_INLINE_FUNCTION + auto get_max_counts() { + int max_counts; + Kokkos::Max max_reduce(max_counts); + Kokkos::parallel_reduce( + "custom_verlet_list::reduce_max", + Kokkos::RangePolicy(0, counts.size()), + [&](const int i, int &value) { + if (counts(i) > value) + value = counts(i); + }, + max_reduce); + Kokkos::fence(); + return max_counts; + } +}; + +template +class Cabana::NeighborList> { +public: + //! Kokkos memory space. + using memory_space = MemorySpace; + //! Neighbor list type. + using list_type = CustomVerletList; + + //! Get the total number of neighbors across all particles. + KOKKOS_INLINE_FUNCTION + static std::size_t totalNeighbor(list_type const &list) { + std::size_t const num_p = list.counts.size(); + std::size_t total_n = 0; + for (std::size_t i = 0; i < num_p; ++i) + total_n += list.counts(i); + return total_n; + } + + //! Get the maximum number of neighbors per particle. + KOKKOS_INLINE_FUNCTION + static std::size_t maxNeighbor(list_type const &list) { + // Stored during neighbor search. + return list.max_n; + } + + //! Get the number of neighbors for a given particle index. + KOKKOS_INLINE_FUNCTION + static std::size_t numNeighbor(list_type const &list, + std::size_t const particle_index) { + return list.counts(particle_index); + } + + //! Get the id for a neighbor for a given particle index and the index of + //! the neighbor relative to the particle. + KOKKOS_INLINE_FUNCTION + static std::size_t getNeighbor(list_type const &list, + std::size_t const particle_index, + std::size_t const count) { + return list.neighbors(particle_index, count); + } +}; + +#endif // SHARED_MEMORY_PARALLELISM diff --git a/src/core/exclusions.hpp b/src/core/exclusions.hpp index b126fcc732b..2ed62db6380 100644 --- a/src/core/exclusions.hpp +++ b/src/core/exclusions.hpp @@ -33,10 +33,14 @@ * calculated. */ inline bool do_nonbonded(Particle const &p1, Particle const &p2) { - /* check for particle 2 in particle 1's exclusion list. The exclusion list is - * symmetric, so this is sufficient. */ - return std::ranges::none_of( + /* check for particle 2 in particle 1's exclusion list. The exclusion list + * should be symmetric, so this is sufficient. */ + /* However. in present implementation, the exclusion list is not symmetric.*/ + bool p1_p2 = std::ranges::none_of( p1.exclusions(), [p2_id = p2.id()](int id) { return id == p2_id; }); + bool p2_p1 = std::ranges::none_of( + p2.exclusions(), [p1_id = p1.id()](int id) { return id == p1_id; }); + return (p1_p2 or p2_p1); } /** Remove exclusion from particle if possible */ diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 475825f5819..52791093212 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -58,6 +58,11 @@ #include #endif +#ifdef SHARED_MEMORY_PARALLELISM +#include "short_range_cabana.hpp" +#include +#endif + #include #include #include @@ -206,40 +211,145 @@ void System::System::calculate_forces() { auto const collision_detection_cutoff = INACTIVE_CUTOFF; #endif - short_range_loop( - [coulomb_kernel_ptr = get_ptr(coulomb_kernel), &bonded_ias = *bonded_ias, - &bond_breakage = *bond_breakage, &box_geo = *box_geo]( - Particle &p1, int bond_id, std::span partners) { - return add_bonded_force(p1, bond_id, partners, bonded_ias, - bond_breakage, box_geo, coulomb_kernel_ptr); - }, - [coulomb_kernel_ptr = get_ptr(coulomb_kernel), - dipoles_kernel_ptr = get_ptr(dipoles_kernel), - elc_kernel_ptr = get_ptr(elc_kernel), - coulomb_u_kernel_ptr = get_ptr(coulomb_u_kernel), - &nonbonded_ias = *nonbonded_ias, &thermostat = *thermostat, - &bonded_ias = *bonded_ias, + // interaction kernel is defined + auto bond_kernel = [coulomb_kernel_ptr = get_ptr(coulomb_kernel), + &bonded_ias = *bonded_ias, + &bond_breakage = *bond_breakage, + &box_geo = *box_geo](Particle &p1, int bond_id, + std::span partners) { + return add_bonded_force(p1, bond_id, partners, bonded_ias, bond_breakage, + box_geo, coulomb_kernel_ptr); + }; + + VerletCriterion<> const verlet_criterion{*this, + cell_structure->get_verlet_skin(), + get_interaction_range(), + coulomb_cutoff, + dipole_cutoff, + collision_detection_cutoff}; + +#ifdef SHARED_MEMORY_PARALLELISM +#ifdef CALIPER + CALI_MARK_BEGIN("parallel short range"); +#endif + using execution_space = Kokkos::DefaultExecutionSpace; + update_cabana_state(*cell_structure, verlet_criterion, + get_interaction_range()); + auto const &unique_particles = cell_structure->get_unique_particles(); + auto const &local_force = cell_structure->get_local_force(); +#ifdef ROTATION + auto const &local_torque = cell_structure->get_local_torque(); +#endif +#ifdef NPT + auto const &local_virial = cell_structure->get_local_virial(); +#endif + auto const &aosoa = cell_structure->get_aosoa(); + + ForcesKernel first_neighbor_kernel( + *bonded_ias, *nonbonded_ias, get_ptr(coulomb_kernel), +#if defined(LONG_RANGE_KERNELS) + get_ptr(dipoles_kernel), get_ptr(elc_kernel), get_ptr(coulomb_u_kernel), + *thermostat, +#endif + *box_geo, +#if defined(LONG_RANGE_KERNELS) or defined(EXCLUSIONS) + unique_particles, +#endif + local_force, +#ifdef ROTATION + local_torque, +#endif +#ifdef NPT + local_virial, +#endif + aosoa); + + cabana_short_range(bond_kernel, first_neighbor_kernel, *cell_structure, + get_interaction_range(), bonded_ias->maximal_cutoff()); + // Force and Torque reduction + int num_threads = execution_space().concurrency(); + Kokkos::RangePolicy policy(0, unique_particles.size()); + Kokkos::parallel_for("reduction", policy, + [&local_force, +#ifdef ROTATION + &local_torque, +#endif + &unique_particles, num_threads](const int i) { + Utils::Vector3d force{}; +#ifdef ROTATION + Utils::Vector3d torque{}; +#endif + for (int tid = 0; tid < num_threads; ++tid) { + force[0] += local_force(i, tid, 0); + force[1] += local_force(i, tid, 1); + force[2] += local_force(i, tid, 2); +#ifdef ROTATION + torque[0] += local_torque(i, tid, 0); + torque[1] += local_torque(i, tid, 1); + torque[2] += local_torque(i, tid, 2); +#endif + } + unique_particles.at(i)->force() += force; +#ifdef ROTATION + unique_particles.at(i)->torque() += torque; +#endif + }); + Kokkos::fence(); + +#ifdef NPT + Utils::Vector3d virial{}; + for (int tid = 0; tid < num_threads; ++tid) { + virial[0] += local_virial(tid, 0); + virial[1] += local_virial(tid, 1); + virial[2] += local_virial(tid, 2); + } + npt_add_virial_force_contribution(virial); +#endif + +#ifdef COLLISION_DETECTION + auto collision_kernel = [&collision_detection = *collision_detection]( + Particle const &p1, Particle const &p2, + Distance const &d) { + if (not collision_detection.is_off()) { + collision_detection.detect_collision(p1, p2, d.dist2); + } + }; + cell_structure->non_bonded_loop(collision_kernel, verlet_criterion); +#endif + +#ifdef CALIPER + CALI_MARK_END("parallel short range"); +#endif + +#else // SHARED_MEMORY_PARALLELISM + + auto pair_kernel = [coulomb_kernel_ptr = get_ptr(coulomb_kernel), + dipoles_kernel_ptr = get_ptr(dipoles_kernel), + elc_kernel_ptr = get_ptr(elc_kernel), + coulomb_u_kernel_ptr = get_ptr(coulomb_u_kernel), + &nonbonded_ias = *nonbonded_ias, + &thermostat = *thermostat, &bonded_ias = *bonded_ias, #ifdef COLLISION_DETECTION - &collision_detection = *collision_detection, -#endif - &box_geo = *box_geo](Particle &p1, Particle &p2, Distance const &d) { - auto const &ia_params = - nonbonded_ias.get_ia_param(p1.type(), p2.type()); - add_non_bonded_pair_force(p1, p2, d.vec21, sqrt(d.dist2), d.dist2, - ia_params, thermostat, box_geo, bonded_ias, - coulomb_kernel_ptr, dipoles_kernel_ptr, - elc_kernel_ptr, coulomb_u_kernel_ptr); + &collision_detection = *collision_detection, +#endif + &box_geo = *box_geo](Particle &p1, Particle &p2, + Distance const &d) { + auto const &ia_params = nonbonded_ias.get_ia_param(p1.type(), p2.type()); + add_non_bonded_pair_force( + p1, p2, d.vec21, sqrt(d.dist2), d.dist2, p1.q() * p2.q(), ia_params, + thermostat, box_geo, bonded_ias, coulomb_kernel_ptr, dipoles_kernel_ptr, + elc_kernel_ptr, coulomb_u_kernel_ptr); #ifdef COLLISION_DETECTION - if (not collision_detection.is_off()) { - collision_detection.detect_collision(p1, p2, d.dist2); - } + if (not collision_detection.is_off()) { + collision_detection.detect_collision(p1, p2, d.dist2); + } #endif - }, - *cell_structure, maximal_cutoff(), bonded_ias->maximal_cutoff(), - VerletCriterion<>{*this, cell_structure->get_verlet_skin(), - get_interaction_range(), coulomb_cutoff, dipole_cutoff, - collision_detection_cutoff}); + }; + short_range_loop(bond_kernel, pair_kernel, *cell_structure, maximal_cutoff(), + bonded_ias->maximal_cutoff(), verlet_criterion); + +#endif // SHARED_MEMORY_PARALLELISM constraints->add_forces(particles, get_sim_time()); oif_global->calculate_forces(); @@ -298,6 +408,9 @@ void calc_long_range_forces(const ParticleRange &particles) { } #ifdef NPT +void npt_add_virial_force_contribution(const Utils::Vector3d &virial) { + ::System::get_system().npt_add_virial_contribution(virial); +} void npt_add_virial_force_contribution(const Utils::Vector3d &force, const Utils::Vector3d &d) { ::System::get_system().npt_add_virial_contribution(force, d); diff --git a/src/core/forces.hpp b/src/core/forces.hpp index 3ba42ea3019..46e1d0ee27e 100644 --- a/src/core/forces.hpp +++ b/src/core/forces.hpp @@ -43,6 +43,7 @@ void calc_long_range_forces(ParticleRange const &particles); #ifdef NPT /** Update the NpT virial */ +void npt_add_virial_force_contribution(Utils::Vector3d const &virial); void npt_add_virial_force_contribution(Utils::Vector3d const &force, Utils::Vector3d const &d); void npt_add_virial_diagonalSum_contribution(double diagonal_sum); diff --git a/src/core/forces_cabana.hpp b/src/core/forces_cabana.hpp new file mode 100644 index 00000000000..1b3d24a4158 --- /dev/null +++ b/src/core/forces_cabana.hpp @@ -0,0 +1,186 @@ +/* + * Copyright (C) 2025 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#include "config/config.hpp" + +#ifdef SHARED_MEMORY_PARALLELISM + +#include "aosoa_pack.hpp" +#include "forces_inline.hpp" + +#include + +#include + +#if defined(__GNUG__) or defined(__clang__) +#define ESPRESSO_ATTR_ALWAYS_INLINE [[gnu::always_inline]] +#else +#define ESPRESSO_ATTR_ALWAYS_INLINE +#endif + +struct ForcesKernel { + BondedInteractionsMap const &bonded_ias; + InteractionsNonBonded const &nonbonded_ias; + Coulomb::ShortRangeForceKernel::kernel_type const *const coulomb_kernel; +#if defined(LONG_RANGE_KERNELS) + Dipoles::ShortRangeForceKernel::kernel_type const *const dipoles_kernel; + Coulomb::ShortRangeForceCorrectionsKernel::kernel_type const *elc_kernel; + Coulomb::ShortRangeEnergyKernel::kernel_type const *const coulomb_u_kernel; + Thermostat::Thermostat const &thermostat; +#endif + BoxGeometry const &box_geo; +#if defined(LONG_RANGE_KERNELS) or defined(EXCLUSIONS) + std::vector const &unique_particles; +#endif + CellStructure::ForceType const &local_force; +#ifdef ROTATION + CellStructure::ForceType const &local_torque; +#endif +#ifdef NPT + CellStructure::VirialType const &local_virial; +#endif + CellStructure::AoSoA_pack const &aosoa; + + ForcesKernel( + BondedInteractionsMap const &bonded_ias_, + InteractionsNonBonded const &nonbonded_ias_, + Coulomb::ShortRangeForceKernel::kernel_type const *coulomb_kernel_, +#if defined(LONG_RANGE_KERNELS) + Dipoles::ShortRangeForceKernel::kernel_type const *dipoles_kernel_, + Coulomb::ShortRangeForceCorrectionsKernel::kernel_type const *elc_kernel_, + Coulomb::ShortRangeEnergyKernel::kernel_type const *coulomb_u_kernel_, + Thermostat::Thermostat const &thermostat_, +#endif + BoxGeometry const &box_geo_, +#if defined(LONG_RANGE_KERNELS) or defined(EXCLUSIONS) + std::vector const &unique_particles_, +#endif + CellStructure::ForceType const &local_force_, +#ifdef ROTATION + CellStructure::ForceType const &local_torque_, +#endif +#ifdef NPT + CellStructure::VirialType const &local_virial_, +#endif + CellStructure::AoSoA_pack const &aosoa_) + : bonded_ias(bonded_ias_), nonbonded_ias(nonbonded_ias_), + coulomb_kernel(coulomb_kernel_), +#if defined(LONG_RANGE_KERNELS) + dipoles_kernel(dipoles_kernel_), elc_kernel(elc_kernel_), + coulomb_u_kernel(coulomb_u_kernel_), thermostat(thermostat_), +#endif + box_geo(box_geo_), +#if defined(LONG_RANGE_KERNELS) or defined(EXCLUSIONS) + unique_particles(unique_particles_), +#endif + local_force(local_force_), +#ifdef ROTATION + local_torque(local_torque_), +#endif +#ifdef NPT + local_virial(local_virial_), +#endif + aosoa(aosoa_) { + } + + ESPRESSO_ATTR_ALWAYS_INLINE KOKKOS_INLINE_FUNCTION void + operator()(int i, int j) const { + + auto const thread_id = omp_get_thread_num(); + + auto const &ia_params = + nonbonded_ias.get_ia_param(aosoa.type(i), aosoa.type(j)); + + ParticleForce pf{}; +#ifdef NPT + Utils::Vector3d virial{}; +#endif + auto const d = box_geo.get_mi_vector( + aosoa.position(i, 0), aosoa.position(i, 1), aosoa.position(i, 2), + aosoa.position(j, 0), aosoa.position(j, 1), aosoa.position(j, 2)); + auto const dist = d.norm(); + +#if defined(LONG_RANGE_KERNELS) or defined(EXCLUSIONS) + auto &p1 = *unique_particles.at(i); + auto &p2 = *unique_particles.at(j); +#endif + +#ifdef EXCLUSIONS + auto const do_nonbonded_flag = do_nonbonded(p1, p2); +#else +#if defined(LONG_RANGE_KERNELS) + auto constexpr do_nonbonded_flag = true; +#endif // LONG_RANGE_KERNELS +#endif + + if (dist < ia_params.max_cut) { +#ifdef EXCLUSIONS + if (do_nonbonded_flag) { +#endif + pf += calc_central_radial_force(ia_params, d, dist); +#ifdef EXCLUSIONS + } +#endif + } + +#if defined(LONG_RANGE_KERNELS) +#ifdef ELECTROSTATICS + auto const q1q2 = aosoa.charge(i) * aosoa.charge(j); +#else + auto constexpr q1q2 = 0.; +#endif + add_non_bonded_pair_force_with_p(p1, p2, pf, +#ifdef NPT + virial, +#endif // NPT + d, dist, dist * dist, q1q2, ia_params, + do_nonbonded_flag, thermostat, box_geo, + bonded_ias, coulomb_kernel, dipoles_kernel, + elc_kernel, coulomb_u_kernel); +#endif // LONG_RANGE_KERNELS + + local_force(i, thread_id, 0) += pf.f[0]; + local_force(i, thread_id, 1) += pf.f[1]; + local_force(i, thread_id, 2) += pf.f[2]; +#ifdef ROTATION + local_torque(i, thread_id, 0) += pf.torque[0]; + local_torque(i, thread_id, 1) += pf.torque[1]; + local_torque(i, thread_id, 2) += pf.torque[2]; +#endif + + auto const opf = calc_opposing_force(pf, d); + local_force(j, thread_id, 0) += opf.f[0]; + local_force(j, thread_id, 1) += opf.f[1]; + local_force(j, thread_id, 2) += opf.f[2]; +#ifdef ROTATION + local_torque(j, thread_id, 0) += opf.torque[0]; + local_torque(j, thread_id, 1) += opf.torque[1]; + local_torque(j, thread_id, 2) += opf.torque[2]; +#endif +#ifdef NPT + local_virial(thread_id, 0) += virial[0]; + local_virial(thread_id, 1) += virial[1]; + local_virial(thread_id, 2) += virial[2]; +#endif + } +}; + +#endif // SHARED_MEMORY_PARALLELISM diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index 9dcd726c511..381a1934347 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -29,6 +29,10 @@ #include "forces.hpp" +#if defined(ELECTROSTATICS) or defined(DIPOLES) or defined(DPD) or defined(NPT) +#define LONG_RANGE_KERNELS +#endif + #include "BoxGeometry.hpp" #include "actor/visitors.hpp" #include "bond_breakage/bond_breakage.hpp" @@ -169,25 +173,50 @@ inline ParticleForce calc_opposing_force(ParticleForce const &pf, return out; } -/** Calculate non-bonded forces between a pair of particles and update their - * forces and torques. - * @param[in,out] p1 particle 1. - * @param[in,out] p2 particle 2. - * @param[in] d vector between @p p1 and @p p2. - * @param[in] dist distance between @p p1 and @p p2. - * @param[in] dist2 distance squared between @p p1 and @p p2. - * @param[in] ia_params non-bonded interaction kernels. - * @param[in] thermostat thermostat. - * @param[in] box_geo box geometry. - * @param[in] bonded_ias bonded interaction kernels. - * @param[in] coulomb_kernel Coulomb force kernel. - * @param[in] dipoles_kernel Dipolar force kernel. - * @param[in] elc_kernel ELC force correction kernel. - * @param[in] coulomb_u_kernel Coulomb energy kernel. +/** + * For the interaction which need NO particle information */ -inline void add_non_bonded_pair_force( - Particle &p1, Particle &p2, Utils::Vector3d const &d, double dist, - double dist2, IA_parameters const &ia_params, +inline void add_non_bonded_pair_without_p( + ParticleForce &pf, Utils::Vector3d const &d, double dist, double q1q2, + IA_parameters const &ia_params, [[maybe_unused]] bool do_nonbonded, + Coulomb::ShortRangeForceKernel::kernel_type const *coulomb_kernel) { + + /***********************************************/ + /* non-bonded pair potentials */ + /***********************************************/ + + if (dist < ia_params.max_cut) { +#ifdef EXCLUSIONS + if (do_nonbonded) { +#endif + pf += calc_central_radial_force(ia_params, d, dist); +#ifdef EXCLUSIONS + } +#endif + } + + /***********************************************/ + /* short-range electrostatics */ + /***********************************************/ + +#ifdef ELECTROSTATICS + // real-space electrostatic charge-charge interaction + if (q1q2 != 0. and coulomb_kernel != nullptr) { + pf.f += (*coulomb_kernel)(q1q2, d, dist); + } +#endif // ELECTROSTATICS +} + +/** + * @brief For interactions which need particle information. + */ +inline void add_non_bonded_pair_force_with_p( + Particle &p1, Particle &p2, ParticleForce &pf, +#if defined(NPT) and defined(SHARED_MEMORY_PARALLELISM) + Utils::Vector3d &virial, +#endif + Utils::Vector3d const &d, double dist, double dist2, double q1q2, + IA_parameters const &ia_params, [[maybe_unused]] bool do_nonbonded, Thermostat::Thermostat const &thermostat, BoxGeometry const &box_geo, [[maybe_unused]] BondedInteractionsMap const &bonded_ias, Coulomb::ShortRangeForceKernel::kernel_type const *coulomb_kernel, @@ -195,17 +224,14 @@ inline void add_non_bonded_pair_force( Coulomb::ShortRangeForceCorrectionsKernel::kernel_type const *elc_kernel, Coulomb::ShortRangeEnergyKernel::kernel_type const *coulomb_u_kernel) { - ParticleForce pf{}; - /***********************************************/ /* non-bonded pair potentials */ /***********************************************/ if (dist < ia_params.max_cut) { #ifdef EXCLUSIONS - if (do_nonbonded(p1, p2)) { + if (do_nonbonded) { #endif - pf += calc_central_radial_force(ia_params, d, dist); #ifdef THOLE pf.f += thole_pair_force(p1, p2, ia_params, d, dist, bonded_ias, coulomb_kernel); @@ -222,7 +248,11 @@ inline void add_non_bonded_pair_force( /* electrostatic is calculated by energy */ /*********************************************************************/ #ifdef NPT +#ifdef SHARED_MEMORY_PARALLELISM + virial += hadamard_product(pf.f, d); +#else npt_add_virial_force_contribution(pf.f, d); +#endif #endif /***********************************************/ @@ -231,13 +261,16 @@ inline void add_non_bonded_pair_force( #ifdef ELECTROSTATICS // real-space electrostatic charge-charge interaction - auto const q1q2 = p1.q() * p2.q(); if (q1q2 != 0. and coulomb_kernel != nullptr) { pf.f += (*coulomb_kernel)(q1q2, d, dist); #ifdef NPT +#ifdef SHARED_MEMORY_PARALLELISM + virial[0] += (*coulomb_u_kernel)(p1, p2, q1q2, d, dist); +#else npt_add_virial_diagonalSum_contribution( (*coulomb_u_kernel)(p1, p2, q1q2, d, dist)); -#endif +#endif // SHARED_MEMORY_PARALLELISM +#endif // NPT #ifdef P3M if (elc_kernel) (*elc_kernel)(p1, p2, q1q2); @@ -255,8 +288,7 @@ inline void add_non_bonded_pair_force( auto const force = dpd_pair_force(p1.pos(), p1.v(), p1.id(), p2.pos(), p2.v(), p2.id(), *thermostat.dpd, box_geo, ia_params, d, dist, dist2); - p1.force() += force; - p2.force() -= force; + pf += force; } #endif @@ -270,13 +302,80 @@ inline void add_non_bonded_pair_force( pf += (*dipoles_kernel)(p1, p2, d, dist, dist2); } #endif +} + +/** Calculate non-bonded forces between a pair of particles and update their + * forces and torques. + * @param[in,out] p1 particle 1. + * @param[in,out] p2 particle 2. + * @param[in] d vector between @p p1 and @p p2. + * @param[in] dist distance between @p p1 and @p p2. + * @param[in] dist2 distance squared between @p p1 and @p p2. + * @param[in] q1q2 charge x charge between @p p1 and @p p2. + * @param[in] ia_params non-bonded interaction kernels. + * @param[in] thermostat thermostat. + * @param[in] box_geo box geometry. + * @param[in] bonded_ias bonded interaction kernels. + * @param[in] coulomb_kernel Coulomb force kernel. + * @param[in] dipoles_kernel Dipolar force kernel. + * @param[in] elc_kernel ELC force correction kernel. + * @param[in] coulomb_u_kernel Coulomb energy kernel. + */ +inline auto add_non_bonded_pair_force( + Particle &p1, Particle &p2, Utils::Vector3d const &d, double dist, + double dist2, double q1q2, IA_parameters const &ia_params, + Thermostat::Thermostat const &thermostat, BoxGeometry const &box_geo, + [[maybe_unused]] BondedInteractionsMap const &bonded_ias, + Coulomb::ShortRangeForceKernel::kernel_type const *coulomb_kernel, + Dipoles::ShortRangeForceKernel::kernel_type const *dipoles_kernel, + Coulomb::ShortRangeForceCorrectionsKernel::kernel_type const *elc_kernel, + Coulomb::ShortRangeEnergyKernel::kernel_type const *coulomb_u_kernel) { + + ParticleForce pf{}; +#if defined(NPT) and defined(SHARED_MEMORY_PARALLELISM) + Utils::Vector3d virial{}; +#endif + +#ifdef EXCLUSIONS + auto const do_nonbonded_flag = do_nonbonded(p1, p2); +#else +#if defined(LONG_RANGE_KERNELS) + auto constexpr do_nonbonded_flag = true; +#endif // LONG_RANGE_KERNELS +#endif + + if (dist < ia_params.max_cut) { +#ifdef EXCLUSIONS + if (do_nonbonded_flag) { +#endif + pf += calc_central_radial_force(ia_params, d, dist); +#ifdef EXCLUSIONS + } +#endif + } + +#if defined(LONG_RANGE_KERNELS) + add_non_bonded_pair_force_with_p( + p1, p2, pf, +#if defined(NPT) and defined(SHARED_MEMORY_PARALLELISM) + virial, +#endif + d, dist, dist2, q1q2, ia_params, do_nonbonded_flag, thermostat, box_geo, + bonded_ias, coulomb_kernel, dipoles_kernel, elc_kernel, coulomb_u_kernel); +#endif /***********************************************/ /* add total non-bonded forces to particles */ /***********************************************/ +#if defined(NPT) and defined(SHARED_MEMORY_PARALLELISM) + return std::pair{pf, virial}; +#elif defined(SHARED_MEMORY_PARALLELISM) + return pf; +#else p1.force_and_torque() += pf; p2.force_and_torque() += calc_opposing_force(pf, d); +#endif } /** Compute the bonded interaction force between particle pairs. diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 181a0322830..8106d42196e 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -505,7 +505,15 @@ int System::System::integrate(int n_steps, int reuse_forces) { if (propagation.integ_switch != INTEG_METHOD_STEEPEST_DESCENT) { lb_active = lb.is_solver_set(); ek_active = ek.is_ready_for_propagation(); +#ifdef SHARED_MEMORY_PARALLELISM + cell_structure->set_max_prefactor(5); +#endif + } +#ifdef SHARED_MEMORY_PARALLELISM + else { + cell_structure->set_max_prefactor(8); } +#endif auto const calc_md_steps_per_tau = [this](double tau) { return static_cast(std::round(tau / time_step)); }; diff --git a/src/core/npt.cpp b/src/core/npt.cpp index 6a554671a87..dcdaf080a13 100644 --- a/src/core/npt.cpp +++ b/src/core/npt.cpp @@ -148,4 +148,12 @@ void System::System::npt_add_virial_contribution(Utils::Vector3d const &force, npt_inst_pressure->p_vir += hadamard_product(force, d); } } + +void System::System::npt_add_virial_contribution( + Utils::Vector3d const &virial) { + if ((propagation->integ_switch == INTEG_METHOD_NPT_ISO_AND) or + (propagation->integ_switch == INTEG_METHOD_NPT_ISO_MTK)) { + npt_inst_pressure->p_vir += virial; + } +} #endif // NPT diff --git a/src/core/short_range_cabana.hpp b/src/core/short_range_cabana.hpp new file mode 100644 index 00000000000..1f1b7593fe0 --- /dev/null +++ b/src/core/short_range_cabana.hpp @@ -0,0 +1,184 @@ +/* + * Copyright (C) 2025 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#include "config/config.hpp" + +#ifdef SHARED_MEMORY_PARALLELISM + +#include "cell_system/CellStructure.hpp" + +#include "aosoa_pack.hpp" +#include "custom_verlet_list.hpp" +#include "forces_cabana.hpp" + +#include +#include + +#include +#include +#include + +ESPRESSO_ATTR_ALWAYS_INLINE inline void +commit_particle(Particle const &p, int const index, + CellStructure::AoSoA_pack &aosoa) { + aosoa.id(index) = p.id(); +#ifdef ELECTROSTATICS + aosoa.charge(index) = p.q(); +#endif + aosoa.type(index) = p.type(); + auto const &pos = p.pos(); + aosoa.position(index, 0) = pos[0]; + aosoa.position(index, 1) = pos[1]; + aosoa.position(index, 2) = pos[2]; +} + +ESPRESSO_ATTR_ALWAYS_INLINE inline void construct_verlet_list( + std::span cells, BoxGeometry const &box_geo, + CellStructure::ListType &verlet_list, auto const &verlet_criterion, + Kokkos::View const &id_to_index, int const max_id) { + + auto const distance_function = detail::MinimalImageDistance{box_geo}; + + // implementation detail: max_id refers to the max local particle id, + // but ghost particles from other ranks may have larger particle ids; + // -1 is used as a sentinel value for particle ids from other threads + + auto intra_kernel = [&cells, &distance_function, &verlet_criterion, + &id_to_index, &verlet_list, max_id](const int i) { + auto &local_particles = cells[i]->particles(); + for (auto it = local_particles.begin(); it != local_particles.end(); ++it) { + auto const &p1 = *it; + if (p1.id() <= max_id) { + auto const ii = id_to_index(p1.id()); + if (ii >= 0) { + // pairs in this cell + for (auto jt = std::next(it); jt != local_particles.end(); ++jt) { + if ((*jt).id() <= max_id) { + if (verlet_criterion(p1, *jt, distance_function(p1, *jt))) { + auto const jj = id_to_index((*jt).id()); + if (jj >= 0) { + verlet_list.addNeighborLB(ii, jj); + } + } + } + } + } + } + } + }; + + auto inter_kernel = [&cells, &distance_function, &verlet_criterion, + &id_to_index, &verlet_list, max_id](const int i) { + auto &local_particles = cells[i]->particles(); + for (auto it = local_particles.begin(); it != local_particles.end(); ++it) { + auto const &p1 = *it; + if (p1.id() <= max_id) { + auto const ii = id_to_index(p1.id()); + if (ii >= 0) { + // pairs with neighboring cells + for (auto &neighbor : cells[i]->neighbors().red()) { + for (auto const &p2 : neighbor->particles()) { + if (p2.id() <= max_id) { + if (verlet_criterion(p1, p2, distance_function(p1, p2))) { + auto const jj = id_to_index(p2.id()); + if (jj >= 0) { + verlet_list.addNeighbor(ii, jj); + } + } + } + } + } + } + } + } + }; + + Kokkos::parallel_for("inter", cells.size(), intra_kernel); + Kokkos::fence(); + + Kokkos::parallel_for("intra", cells.size(), inter_kernel); + Kokkos::fence(); +} + +ESPRESSO_ATTR_ALWAYS_INLINE inline void +update_cabana_state(CellStructure &cell_structure, auto const &verlet_criterion, + double const pair_cutoff) { + using execution_space = Kokkos::DefaultExecutionSpace; + auto const num_threads = execution_space().concurrency(); + auto const rebuild = + cell_structure.prepare_verlet_list_cabana(num_threads, pair_cutoff); + auto const &unique_particles = cell_structure.get_unique_particles(); + auto const max_id = cell_structure.get_cached_max_local_particle_id(); + auto &aosoa = cell_structure.get_aosoa(); + + // =================================================== + // Fill particle storage + // =================================================== + Kokkos::View id_to_index( + Kokkos::ViewAllocateWithoutInitializing("id_to_index"), max_id + 1); + Kokkos::deep_copy(id_to_index, -1); + + using policy_type = Kokkos::RangePolicy; + Kokkos::parallel_for( + "AoSoA write", policy_type(0, unique_particles.size()), + [&unique_particles, &aosoa, &id_to_index](int const index) { + auto const &p = *unique_particles.at(index); + commit_particle(p, index, aosoa); + id_to_index(p.id()) = index; + }); + Kokkos::fence(); + + // =================================================== + // Get Verlet pairs and fill Verlet list + // =================================================== + if (rebuild) { + cell_structure.rebuild_verlet_list_cabana( + [&](std::span cells, BoxGeometry const &box, + CellStructure::ListType &verlet_list) { + construct_verlet_list(std::move(cells), box, verlet_list, + verlet_criterion, id_to_index, max_id); + }); + } +} + +void cabana_short_range(auto const &bond_kernel, auto const &forces_kernel, + CellStructure &cell_structure, double pair_cutoff, + double bond_cutoff) { + using execution_space = Kokkos::DefaultExecutionSpace; + assert(cell_structure.get_resort_particles() == Cells::RESORT_NONE); + + if (bond_cutoff >= 0.) { + cell_structure.bond_loop(bond_kernel); + } + + // Cabana short range loop + if (pair_cutoff > 0.) { + auto const &verlet_list = cell_structure.get_verlet_list_cabana(); + Kokkos::RangePolicy policy( + 0, cell_structure.get_unique_particles().size()); + Cabana::neighbor_parallel_for(policy, forces_kernel, verlet_list, + Cabana::FirstNeighborsTag(), + Cabana::SerialOpTag()); + Kokkos::fence(); + } +} + +#endif // SHARED_MEMORY_PARALLELISM diff --git a/src/core/system/System.cpp b/src/core/system/System.cpp index 56abf6c4d3a..c96819568a4 100644 --- a/src/core/system/System.cpp +++ b/src/core/system/System.cpp @@ -66,6 +66,9 @@ System::System(Private) { box_geo = std::make_shared(); local_geo = std::make_shared(); cell_structure = std::make_shared(*box_geo); +#ifdef SHARED_MEMORY_PARALLELISM + cell_structure->set_kokkos_handle(::kokkos_handle); +#endif propagation = std::make_shared(); bonded_ias = std::make_shared(); thermostat = std::make_shared(); diff --git a/src/core/system/System.hpp b/src/core/system/System.hpp index 9ac77f21608..be53a43e7b0 100644 --- a/src/core/system/System.hpp +++ b/src/core/system/System.hpp @@ -38,7 +38,7 @@ class BoxGeometry; class LocalBox; -struct CellStructure; +class CellStructure; class Propagation; class InteractionsNonBonded; class BondedInteractionsMap; @@ -165,6 +165,7 @@ class System : public std::enable_shared_from_this { void npt_add_virial_contribution(double energy); void npt_add_virial_contribution(Utils::Vector3d const &force, Utils::Vector3d const &d); + void npt_add_virial_contribution(Utils::Vector3d const &virial); #endif // NPT /** @brief Calculate all forces. */ diff --git a/src/script_interface/code_info/CodeInfo.cpp b/src/script_interface/code_info/CodeInfo.cpp index 868988b4d08..4f3416f1e8b 100644 --- a/src/script_interface/code_info/CodeInfo.cpp +++ b/src/script_interface/code_info/CodeInfo.cpp @@ -61,6 +61,13 @@ Variant CodeInfo::do_call_method(std::string const &name, #else // SCAFACOS return make_vector_of_variants(std::vector(0)); #endif // SCAFACOS + } + if (name == "has_fast_math") { +#if defined(__FAST_MATH__) + return true; +#else + return false; +#endif } return {}; } diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 3039d83629b..517a10dc1c5 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -95,11 +95,12 @@ function(python_test) ${MPIEXEC} ${ESPRESSO_MPIEXEC_PREFLAGS} ${MPIEXEC_NUMPROC_FLAG} ${TEST_NUM_PROC} ${MPIEXEC_PREFLAGS} ${ESPRESSO_MPIEXEC_TMPDIR} ${CMAKE_BINARY_DIR}/pypresso ${PYPRESSO_OPTIONS} - ${TEST_FILE_CONFIGURED} ${TEST_ARGUMENTS} ${MPIEXEC_POSTFLAGS}) + ${TEST_FILE_CONFIGURED} --verbose ${TEST_ARGUMENTS} + ${MPIEXEC_POSTFLAGS}) else() add_test(NAME ${TEST_NAME} COMMAND ${CMAKE_BINARY_DIR}/pypresso ${PYPRESSO_OPTIONS} - ${TEST_FILE_CONFIGURED} ${TEST_ARGUMENTS}) + ${TEST_FILE_CONFIGURED} --verbose ${TEST_ARGUMENTS}) endif() if(${TEST_GPU_SLOTS} GREATER 0 AND ESPRESSO_BUILD_WITH_CUDA) diff --git a/testsuite/python/caliper.py b/testsuite/python/caliper.py index 265e2f741e7..a7dc44784cb 100644 --- a/testsuite/python/caliper.py +++ b/testsuite/python/caliper.py @@ -25,21 +25,23 @@ import sys import os -EXPECTED_LABELS = """ +EXPECTED_LABELS = f""" integrate Initial Force Calculation calculate_forces copy_particles_to_GPU init_forces_and_thermost calc_long_range_forces - short_range_loop + {'parallel short range' if espressomd.has_features( + ["SHARED_MEMORY_PARALLELISM"]) else 'short_range_loop'} copy_forces_from_GPU Integration loop calculate_forces copy_particles_to_GPU init_forces_and_thermost calc_long_range_forces - short_range_loop + {'parallel short range' if espressomd.has_features( + ["SHARED_MEMORY_PARALLELISM"]) else 'short_range_loop'} copy_forces_from_GPU calc_energies short_range_loop diff --git a/testsuite/python/integrator_exceptions.py b/testsuite/python/integrator_exceptions.py index f8ec6a21b8c..bc3c573f6e1 100644 --- a/testsuite/python/integrator_exceptions.py +++ b/testsuite/python/integrator_exceptions.py @@ -21,6 +21,7 @@ import espressomd.lees_edwards import espressomd.shapes import espressomd.propagation +import os import numpy as np import unittest as ut import unittest_decorators as utx @@ -33,6 +34,7 @@ class Test(ut.TestCase): msg = r'while calling method integrate\(\): ERROR: ' def setUp(self): + self.system.box_l = [1., 1., 1.] self.system.part.add(pos=(0, 0, 0)) self.system.integrator.set_vv() self.system.periodicity = 3 * [True] @@ -179,6 +181,59 @@ def test_npt_integrator(self): self.system.lees_edwards.protocol = None self.system.integrator.run(0) + @ut.skipIf(espressomd.code_info._CodeInfo().call_method("has_fast_math"), + "cannot run with fast-math optimizations") + @ut.skipIf(os.environ.get("UBSAN_OPTIONS"), + "cannot run with UBSAN instrumentation") + @ut.skipIf(espressomd.has_features("FPE"), + "cannot run with FPE instrumentation") + @utx.skipIfMissingFeatures(["NPT", "WCA"]) + def test_npt_integrator_negative_volume(self): + """Test for NpT with bad parameters.""" + + import tests_common + data = np.genfromtxt(tests_common.data_path("npt_lj_system.data")) + ref_box_l = np.max(data[:, 0:3]) + + system = self.system + system.part.clear() + system.cell_system.skin = 0. + + for barostat in ["Andersen", "MTK"]: + system.box_l = 3 * [ref_box_l] + system.time_step = 0.01 + if barostat == "Andersen": + piston = 1e-4 + else: + piston = 4.0 + direction = [True] * 3 + ext_pressure = 100.0 # Too large external pressure + system.part.add(pos=data[:, 0:3], v=data[:, 3:6]) + system.integrator.set_vv() + system.thermostat.set_npt(kT=1.0, gamma0=0.1, gammav=1e-3, seed=42) + system.integrator.set_isotropic_npt(ext_pressure=ext_pressure, + piston=piston, + direction=direction, + barostat=barostat) + + if barostat == "Andersen": + exception_msg = "" + try: + system.integrator.run(10) + except Exception as err: + exception_msg = f"{exception_msg}\n{err}" + try: + system.part.clear() + except Exception as err: + exception_msg = f"{exception_msg}\n{err}" + self.assertIn("the volume to become negative", exception_msg) + if barostat == "MTK": + # Volume cannot be negative within NPT ensemble based on MTK equation + self.assertGreater(float(np.prod(system.box_l)), 0.) + + system.part.clear() + system.box_l = [1., 1., 1.] + @utx.skipIfMissingFeatures("STOKESIAN_DYNAMICS") def test_stokesian_integrator(self): self.system.cell_system.skin = 0.4 diff --git a/testsuite/python/integrator_npt_stats.py b/testsuite/python/integrator_npt_stats.py index ed352fe1a1c..e6f9c7281a7 100644 --- a/testsuite/python/integrator_npt_stats.py +++ b/testsuite/python/integrator_npt_stats.py @@ -44,6 +44,7 @@ def setUp(self): def tearDown(self): self.system.part.clear() + self.system.non_bonded_inter.reset() self.system.thermostat.turn_off() self.system.integrator.set_vv() @@ -116,44 +117,6 @@ def test_compressibility_and_pressure(self): self.assertAlmostEqual(avpV_sim, 100., delta=1.) self.assertAlmostEqual(avpV_inst, 100., delta=1.) - def test_negative_volume(self): - """Test for NpT with bad parameters.""" - - data = np.genfromtxt(tests_common.data_path("npt_lj_system.data")) - ref_box_l = np.max(data[:, 0:3]) - - system = self.system - system.box_l = 3 * [ref_box_l] - dt = 0.01 - system.time_step = dt - if self.barostat == "Andersen": - piston = 0.0001 - else: - piston = 4.0 - - direction = [True] * 3 - ext_pressure = 100.0 # Too large external pressure - system.box_l = 3 * [ref_box_l] - system.part.add(pos=data[:, 0:3], type=len(data) * [2]) - system.part.all().pos = data[:, 0:3] - system.part.all().v = data[:, 3:6] - self.system.integrator.set_vv() - - system.thermostat.set_npt(kT=1.0, gamma0=0.1, gammav=0.001, seed=42) - system.integrator.set_isotropic_npt(ext_pressure=ext_pressure, - piston=piston, - direction=direction, - barostat=self.barostat) - - if self.barostat == "Andersen": - with self.assertRaises(Exception): - system.integrator.run(10) - elif self.barostat == "MTK": - with self.assertRaises(Exception): - system.integrator.run(10) - # Volume cannot be negative within NPT ensemble based on MTK equation - self.assertTrue(float(np.prod(system.box_l)) > 0.) - @utx.skipIfMissingFeatures("NPT") class IntegratorNPT_Andersen(IntegratorNPT, ut.TestCase): diff --git a/testsuite/python/npt_thermostat.py b/testsuite/python/npt_thermostat.py index aa4fc7b0219..3a1733e7e18 100644 --- a/testsuite/python/npt_thermostat.py +++ b/testsuite/python/npt_thermostat.py @@ -34,11 +34,16 @@ class NPTThermostat: def setUp(self): np.random.seed(42) + self.system.box_l = [2., 2., 2.] + self.system.time_step = 0.01 def tearDown(self): + self.system.non_bonded_inter.reset() self.system.part.clear() self.system.thermostat.turn_off() self.system.integrator.set_vv() + if espressomd.has_features("ELECTROSTATICS"): + self.system.electrostatics.clear() def test_01__rng(self): """Test for RNG consistency.""" @@ -200,6 +205,47 @@ def test_integrator_exceptions(self): system.integrator.set_isotropic_npt(ext_pressure=1., piston=1., direction=[0, 0, 0], barostat=self.barostat) + @utx.skipIfMissingFeatures(["WCA", "P3M"]) + def test_pressure_with_p3m(self): + """Test for NpT with P3M.""" + + data = np.genfromtxt(tests_common.data_path("npt_lj_system.data")) + ref_box_l = np.max(data[:, 0:3]) + p_ext = 1.0 + + system = self.system + system.box_l = 3 * [ref_box_l] + system.time_step = 0.01 + system.non_bonded_inter[2, 2].wca.set_params(epsilon=1., sigma=1.) + system.part.add(pos=data[:, 0:3], v=data[:, 3:6], type=len(data) * [2], + q=np.sign(np.arange(100) - 50 + 0.5)) + system.integrator.set_vv() + system.electrostatics.solver = espressomd.electrostatics.P3M( + prefactor=2.0, accuracy=1e-2, mesh=3 * [18], cao=5, tune=True) + + if self.barostat == "Andersen": + system.thermostat.set_npt(kT=1.0, gamma0=0.2, gammav=0.01, seed=42) + system.integrator.set_isotropic_npt( + ext_pressure=p_ext, piston=0.0001) + else: + system.thermostat.set_npt( + kT=1.0, gamma0=0.5, gammav=0.001, seed=42) + system.integrator.set_isotropic_npt( + ext_pressure=p_ext, piston=4.0, barostat=self.barostat) + + steps = int(0.1 / system.time_step) + + for _ in range(100): + system.integrator.run(steps) + p_sim = system.analysis.pressure()['total'] + p_kin = system.analysis.pressure()['kinetic'] + # virial of electrostatic force from system.analysis + p_vir = p_sim - p_kin + # virial of electrostatic force from instantaneous_pressure + p_inst_vir = system.analysis.get_instantaneous_pressure_virial() + + np.testing.assert_allclose(p_vir, p_inst_vir, rtol=1e-2, atol=1e-7) + @utx.skipIfMissingFeatures("NPT") class NPTThermostat_Andersen(NPTThermostat, ut.TestCase):