diff --git a/doc/tutorials/charged_system/charged_system.ipynb b/doc/tutorials/charged_system/charged_system.ipynb index 1a5c03b1df9..0e0e3641e22 100644 --- a/doc/tutorials/charged_system/charged_system.ipynb +++ b/doc/tutorials/charged_system/charged_system.ipynb @@ -44,7 +44,7 @@ "import espressomd.zn\n", "\n", "\n", - "espressomd.assert_features(['ELECTROSTATICS', 'P3M', 'WCA'])\n", + "espressomd.assert_features(['ELECTROSTATICS', 'P3M', 'WCA', 'EXTERNAL_FORCES'])\n", "\n", "import tqdm\n", "import numpy as np\n", diff --git a/doc/tutorials/convert.py b/doc/tutorials/convert.py index c52a5bdb898..f6cd752244c 100644 --- a/doc/tutorials/convert.py +++ b/doc/tutorials/convert.py @@ -188,7 +188,7 @@ def execute_notebook(nb, src, cell_separator, notebook_filepath): import importlib_wrapper as iw notebook_dirname = os.path.dirname(notebook_filepath) # disable OpenGL GUI - src_no_gui = iw.mock_es_visualization(src) + src_no_gui = iw.mock_es_visualization(src, force_mock=True) # update notebook with new code set_code_cells(nb, src_no_gui.split(cell_separator)) # execute notebook diff --git a/doc/tutorials/electrodes/electrodes_part1.ipynb b/doc/tutorials/electrodes/electrodes_part1.ipynb index 734fe8d7fb3..ebd84d7cee6 100644 --- a/doc/tutorials/electrodes/electrodes_part1.ipynb +++ b/doc/tutorials/electrodes/electrodes_part1.ipynb @@ -192,7 +192,7 @@ "import espressomd.electrostatics\n", "import espressomd.electrostatic_extensions\n", "\n", - "espressomd.assert_features(['ELECTROSTATICS'])\n", + "espressomd.assert_features(['ELECTROSTATICS', 'EXTERNAL_FORCES'])\n", "plt.rcParams.update({'font.size': 18})" ] }, @@ -251,11 +251,7 @@ "ICC_EPSILON_WALLS = 1e5 # epsilon outside the slit. Very large to mimic metal\n", "ICC_CONVERGENCE = 1e-3 # ICC numeric/performance parameters\n", "ICC_RELAXATION = 0.95\n", - "ICC_MAX_ITERATIONS = 1000\n", - "\n", - "# Lennard-Jones parameters\n", - "LJ_SIGMA = 1.0\n", - "LJ_EPSILON = 1.0 \n", + "ICC_MAX_ITERATIONS = 100\n", "\n", "# Particle parameters\n", "TYPES = {\"Cation\": 0, \"Anion\": 1 ,\"Electrodes\": 2}\n", @@ -309,7 +305,7 @@ "outputs": [], "source": [ "# SOLUTION CELL\n", - "elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=ELC_GAP, maxPWerror=MAX_PW_ERROR)" + "elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=ELC_GAP, maxPWerror=MAX_PW_ERROR, check_neutrality=False)" ] }, { diff --git a/maintainer/parsing/importlib_wrapper.py b/maintainer/parsing/importlib_wrapper.py index 5bc879338c8..30b8c20123a 100644 --- a/maintainer/parsing/importlib_wrapper.py +++ b/maintainer/parsing/importlib_wrapper.py @@ -41,6 +41,7 @@ def configure_and_import(filepath, script_suffix="", move_to_script_dir=True, mock_visualizer=True, + force_visualizer_mock=True, **parameters): """ Copy a Python script to a new location and alter some lines of code: @@ -69,6 +70,8 @@ def configure_and_import(filepath, if ``True``, substitute the visualizer with a ``Mock`` class in case of ``ImportError`` (use ``False`` if an ``ImportError`` is relevant to your test) + force_visualizer_mock : :obj:`bool` + if ``True``, always substitute the visualizer with a ``Mock`` class move_to_script_dir : :obj:`bool` if ``True``, move to the script's directory (useful when the script needs to load files hardcoded as relative paths, or when files are @@ -101,7 +104,7 @@ def configure_and_import(filepath, code = disable_matplotlib_gui(code) # disable OpenGL GUI in case of ImportError using MagicMock() if mock_visualizer: - code = mock_es_visualization(code) + code = mock_es_visualization(code, force_visualizer_mock) # save changes to a new file output_filepath = filepath.parent / \ f"{filepath.stem}_{script_suffix}_processed.py" @@ -433,7 +436,7 @@ def visit_ImportFrom(self, node): node.lineno, node.module, child.name, child.asname) -def mock_es_visualization(code): +def mock_es_visualization(code, force_mock=False): """ Replace ``import espressomd.visualization`` by a ``MagicMock`` when the visualization module is unavailable, by catching the @@ -451,6 +454,12 @@ def mock_es_visualization(code): import unittest.mock import espressomd {1} = unittest.mock.MagicMock() +""".lstrip() + if force_mock: + r_es_vis_mock = r""" +import unittest.mock +import espressomd +{1} = unittest.mock.MagicMock() """.lstrip() visitor = GetEspressomdVisualizerImports() diff --git a/src/config/include/config/config.hpp b/src/config/include/config/config.hpp index a3a8d14ea9f..f26264a1db5 100644 --- a/src/config/include/config/config.hpp +++ b/src/config/include/config/config.hpp @@ -93,4 +93,16 @@ #define MAX_OBJECTS_IN_FLUID 10000 #endif +/** + * @brief Special cutoff value for an inactive interaction. + * Non-bonded potentials that have this cutoff are never evaluated. + */ +inline constexpr double INACTIVE_CUTOFF = -1.; + +/** + * @brief Special cutoff value for an inactive bond. + * Bonds that have this cutoff are never evaluated. + */ +inline constexpr double BONDED_INACTIVE_CUTOFF = -1.; + #endif diff --git a/src/core/analysis/statistics.cpp b/src/core/analysis/statistics.cpp index e12ed4c224c..9ee19f5106c 100644 --- a/src/core/analysis/statistics.cpp +++ b/src/core/analysis/statistics.cpp @@ -149,10 +149,10 @@ Utils::Vector3d calc_linear_momentum(System::System const &system, if (include_particles) { momentum = reduce_over_local_particles( *(system.cell_structure), - [](Particle const &p, Utils::Vector3d &res) { - res += p.mass() * p.v(); + [](Utils::Vector3d &acc, Particle const &p) { + acc += p.mass() * p.v(); }, - [](Utils::Vector3d &a, Utils::Vector3d const &b) { a = a + b; }); + [](Utils::Vector3d &acc, Utils::Vector3d const &v) { acc = acc + v; }); } if (include_lbfluid and system.lb.is_solver_set()) { momentum += system.lb.get_momentum() * system.lb.get_lattice_speed(); diff --git a/src/core/aosoa_pack.hpp b/src/core/aosoa_pack.hpp index 3215fb01f88..5c12b606502 100644 --- a/src/core/aosoa_pack.hpp +++ b/src/core/aosoa_pack.hpp @@ -19,6 +19,8 @@ #pragma once +#include + #ifdef SHARED_MEMORY_PARALLELISM #include "cell_system/CellStructure.hpp" diff --git a/src/core/bonded_interactions/bonded_interaction_data.hpp b/src/core/bonded_interactions/bonded_interaction_data.hpp index a2a969bef71..4a6504d6091 100644 --- a/src/core/bonded_interactions/bonded_interaction_data.hpp +++ b/src/core/bonded_interactions/bonded_interaction_data.hpp @@ -24,6 +24,8 @@ * For more information on how to add new interactions, see @ref bondedIA_new. */ +#include + #include "angle_common.hpp" #include "angle_cosine.hpp" #include "angle_cossquare.hpp" @@ -56,11 +58,6 @@ #include #include -/* Special cutoff value for a disabled bond. - * Bonds that have this cutoff are not visited during bond evaluation. - */ -static constexpr double BONDED_INACTIVE_CUTOFF = -1.; - /** Interaction type for unused bonded interaction slots */ struct NoneBond { static constexpr int num = 0; diff --git a/src/core/cell_system/CellStructure.cpp b/src/core/cell_system/CellStructure.cpp index 2e477079ba3..0ce965428be 100644 --- a/src/core/cell_system/CellStructure.cpp +++ b/src/core/cell_system/CellStructure.cpp @@ -29,8 +29,10 @@ #include "BoxGeometry.hpp" #include "LocalBox.hpp" #include "Particle.hpp" +#include "aosoa_pack.hpp" #include "cell_system/CellStructureType.hpp" #include "communication.hpp" +#include "custom_verlet_list.hpp" #include "lees_edwards/lees_edwards.hpp" #include "particle_enumeration.hpp" #include "particle_reduction.hpp" @@ -43,6 +45,13 @@ #include +#ifdef SHARED_MEMORY_PARALLELISM +#include +#include +#include +#include +#endif + #include #include #include @@ -57,45 +66,29 @@ #include #include +CellStructure::~CellStructure() { #ifdef SHARED_MEMORY_PARALLELISM -#include "aosoa_pack.hpp" -#include "custom_verlet_list.hpp" -#include -#include -#include + clear_local_properties(); + // Kokkos handle can only be freed after all Cabana containers have been freed + m_kokkos_handle.reset(); #endif +} -CellStructure::~CellStructure() { #ifdef SHARED_MEMORY_PARALLELISM - if (m_local_force) { - m_local_force.reset(); - } +void CellStructure::clear_local_properties() { + m_local_force.reset(); #ifdef ROTATION - if (m_local_torque) { - m_local_torque.reset(); - } + 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(); + m_local_virial.reset(); #endif + m_aosoa.reset(); + m_particle_storage.reset(); + m_verlet_list_cabana.reset(); + m_rebuild_verlet_list_cabana = true; } -#ifdef SHARED_MEMORY_PARALLELISM - void CellStructure::set_kokkos_handle( std::shared_ptr handle) { m_kokkos_handle = std::move(handle); @@ -123,9 +116,10 @@ static auto estimate_max_counts(int max_prefactor, double pair_cutoff, return max_counts; } -void CellStructure::rebuild_local_properties(std::size_t const num_threads, - double const pair_cutoff) { +void CellStructure::rebuild_local_properties(double const pair_cutoff) { assert(m_kokkos_handle); + using execution_space = Kokkos::DefaultExecutionSpace; + auto const num_threads = execution_space().concurrency(); auto const num_part = get_unique_particles().size(); m_local_force = std::make_unique("local_force", num_part, num_threads); @@ -145,13 +139,17 @@ void CellStructure::rebuild_local_properties(std::size_t const num_threads, m_verlet_list_cabana = std::make_unique(0ul, num_part, max_counts); } +void CellStructure::reset_local_force() { + Kokkos::deep_copy(get_local_force(), 0.); +} + void CellStructure::reset_local_properties() { - Kokkos::deep_copy(get_local_force(), 0); + Kokkos::deep_copy(get_local_force(), 0.); #ifdef ROTATION - Kokkos::deep_copy(get_local_torque(), 0); + Kokkos::deep_copy(get_local_torque(), 0.); #endif #ifdef NPT - Kokkos::deep_copy(get_local_virial(), 0); + Kokkos::deep_copy(get_local_virial(), 0.); #endif } @@ -490,12 +488,12 @@ bool CellStructure::check_resort_required( Utils::Vector3d const &additional_offset) const { auto const lim = Utils::sqr(m_verlet_skin / 2.) - additional_offset.norm2(); - Reduction::AddPartialResultKernel add_partial = [lim](Particle const &p, - bool &result) { - if ((p.pos() - p.pos_at_last_verlet_update()).norm2() > lim) { - result = true; - } - }; + Reduction::AddPartialResultKernel add_partial = + [lim](bool &result, Particle const &p) { + if ((p.pos() - p.pos_at_last_verlet_update()).norm2() > lim) { + result = true; + } + }; Reduction::ReductionOp reduce_op = [](bool &acc, bool const &val) { acc |= val; diff --git a/src/core/cell_system/CellStructure.hpp b/src/core/cell_system/CellStructure.hpp index 32b8eaa5ae9..9358d72ff7d 100644 --- a/src/core/cell_system/CellStructure.hpp +++ b/src/core/cell_system/CellStructure.hpp @@ -729,8 +729,9 @@ class CellStructure : public System::Leaf { 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 rebuild_local_properties(double pair_cutoff); void reset_local_properties(); + void reset_local_force(); auto &get_local_force() { return *m_local_force; } #ifdef ROTATION @@ -742,6 +743,7 @@ class CellStructure : public System::Leaf { 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; } + void clear_local_properties(); [[nodiscard]] auto is_verlet_list_cabana_rebuild_needed() const { return m_rebuild_verlet_list_cabana or (not use_verlet_list); @@ -749,17 +751,16 @@ class CellStructure : public System::Leaf { /** * @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) { + [[nodiscard]] auto prepare_verlet_list_cabana(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); + rebuild_local_properties(cutoff); } else { // If we do not rebuild we can use the saved map reset_local_properties(); diff --git a/src/core/cell_system/HybridDecomposition.cpp b/src/core/cell_system/HybridDecomposition.cpp index e41d7b0e2ad..c6934e08070 100644 --- a/src/core/cell_system/HybridDecomposition.cpp +++ b/src/core/cell_system/HybridDecomposition.cpp @@ -154,10 +154,7 @@ void HybridDecomposition::resort(bool global, m_regular_decomposition.resort(global, diff); m_n_square.resort(global, diff); - /* basically do CellStructure::ghost_count() */ ghost_communicator(exchange_ghosts_comm(), m_box, GHOSTTRANS_PARTNUM); - - /* basically do CellStructure::ghost_update(unsigned data_parts) */ ghost_communicator(exchange_ghosts_comm(), m_box, map_data_parts(m_get_global_ghost_flags())); } diff --git a/src/core/cell_system/particle_enumeration.hpp b/src/core/cell_system/particle_enumeration.hpp index fa9e2adc575..73208761753 100644 --- a/src/core/cell_system/particle_enumeration.hpp +++ b/src/core/cell_system/particle_enumeration.hpp @@ -19,29 +19,27 @@ #pragma once +#include + #include "Cell.hpp" #include "CellStructure.hpp" -#include "config/config.hpp" #ifdef SHARED_MEMORY_PARALLELISM #include #endif #include -#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 + * @tparam Kernel Callable with signature void(std::size_t, Particle&) + * @param cs The cell structure 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) { @@ -49,26 +47,21 @@ inline void enumerate_local_particles(CellStructure const &cs, 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); + std::vector cell_offsets(local_cells.size(), std::size_t{0}); + std::exclusive_scan(local_cells.begin(), local_cells.end(), + cell_offsets.begin(), std::size_t{0}, + [](auto acc, auto const &cell) { + return acc + cell->particles().size(); + }); - // 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)); + auto const n_part = cell_particles.size(); + for (std::size_t p_index{0}; p_index < n_part; ++p_index) { + auto global_index = base_offset + p_index; + kernel(global_index, *(cell_particles.begin() + p_index)); } }); return; diff --git a/src/core/custom_verlet_list.hpp b/src/core/custom_verlet_list.hpp index d68e50e7909..8a80ceac828 100644 --- a/src/core/custom_verlet_list.hpp +++ b/src/core/custom_verlet_list.hpp @@ -19,6 +19,8 @@ #pragma once +#include + #ifdef SHARED_MEMORY_PARALLELISM #include @@ -95,14 +97,14 @@ class CustomVerletList : public Cabana::VerletList(0, counts.size()), - [&](const int i) { - const int count = counts(i); - auto *ptr = &neighbors(i, 0); - std::sort(ptr, ptr + count); - }); + Kokkos::parallel_for("custom_verlet_list::sort_neighbors", + Kokkos::RangePolicy( + std::size_t{0}, counts.size()), + [&](std::size_t const i) { + auto const count = counts(i); + auto *ptr = &neighbors(i, 0); + std::sort(ptr, ptr + count); + }); Kokkos::fence(); } @@ -135,8 +137,9 @@ class CustomVerletList : public Cabana::VerletList max_reduce(max_counts); Kokkos::parallel_reduce( "custom_verlet_list::reduce_max", - Kokkos::RangePolicy(0, counts.size()), - [&](const int i, int &value) { + Kokkos::RangePolicy(std::size_t{0}, + counts.size()), + [&](std::size_t const i, int &value) { if (counts(i) > value) value = counts(i); }, diff --git a/src/core/electrostatics/coulomb.cpp b/src/core/electrostatics/coulomb.cpp index bbddc9dedd0..ef3536e002e 100644 --- a/src/core/electrostatics/coulomb.cpp +++ b/src/core/electrostatics/coulomb.cpp @@ -172,7 +172,7 @@ double Solver::cutoff() const { if (impl->solver) { return std::visit(ShortRangeCutoff(), *impl->solver); } - return -1.0; + return INACTIVE_CUTOFF; } struct EventOnObservableCalc { diff --git a/src/core/electrostatics/elc.hpp b/src/core/electrostatics/elc.hpp index 2a7fe034064..dd330e15e22 100644 --- a/src/core/electrostatics/elc.hpp +++ b/src/core/electrostatics/elc.hpp @@ -115,13 +115,41 @@ struct elc_data { /** The space that is finally left. */ double space_box; - /// pairwise contributions from the lowest and top layers - template - void dielectric_layers_contribution(CoulombP3M const &, - BoxGeometry const &box_geo, +#ifdef SHARED_MEMORY_PARALLELISM + static auto copy_aosoa_vector_elc(std::size_t i, auto &slice) { + return Utils::Vector3d{slice(i, 0), slice(i, 1), slice(i, 2)}; + } + + /// pairwise contributions from lower and upper layers + void dielectric_layers_contribution(BoxGeometry const &box_geo, + std::size_t p1, std::size_t p2, + auto &aosoa, double q1q2, + auto &&kernel) const { + if (aosoa.position(p1, 2) < space_layer) { + auto const q_eff = delta_mid_bot * q1q2; + auto pos2 = copy_aosoa_vector_elc(p2, aosoa.position); + auto pos1 = copy_aosoa_vector_elc(p1, aosoa.position); + pos1[2] *= -1.; + auto const d = box_geo.get_mi_vector(pos2, pos1); + kernel(q_eff, d); + } + if (aosoa.position(p1, 2) > (box_h - space_layer)) { + auto const q_eff = delta_mid_top * q1q2; + auto const z = 2. * box_h - aosoa.position(p1, 2); + auto pos2 = copy_aosoa_vector_elc(p2, aosoa.position); + auto pos1 = copy_aosoa_vector_elc(p1, aosoa.position); + pos1[2] = 2. * box_h - pos1[2]; + auto const d = box_geo.get_mi_vector(pos2, pos1); + kernel(q_eff, d); + } + } +#endif // SHARED_MEMORY_PARALLELISM + + /// pairwise contributions from lower and upper layers + void dielectric_layers_contribution(BoxGeometry const &box_geo, Utils::Vector3d const &pos1, Utils::Vector3d const &pos2, double q1q2, - Kernel &&kernel) const { + auto &&kernel) const { if (pos1[2] < space_layer) { auto const q_eff = delta_mid_bot * q1q2; auto const d = box_geo.get_mi_vector(pos2, {pos1[0], pos1[1], -pos1[2]}); @@ -142,7 +170,7 @@ struct elc_data { auto energy = 0.; for (auto const &p : particles) { dielectric_layers_contribution( - p3m, box_geo, p.pos(), p.pos(), Utils::sqr(p.q()), + box_geo, p.pos(), p.pos(), Utils::sqr(p.q()), [&](double q1q2, Utils::Vector3d const &d) { energy += p3m.pair_energy(q1q2, d.norm()); }); @@ -156,7 +184,7 @@ struct elc_data { ParticleRange const &particles) const { for (auto &p : particles) { dielectric_layers_contribution( - p3m, box_geo, p.pos(), p.pos(), Utils::sqr(p.q()), + box_geo, p.pos(), p.pos(), Utils::sqr(p.q()), [&](double q1q2, Utils::Vector3d const &d) { p.force() += p3m.pair_force(q1q2, d, d.norm()); }); @@ -249,6 +277,34 @@ struct ElectrostaticLayerCorrection return {}; } +#ifdef SHARED_MEMORY_PARALLELISM + /** @brief Calculate short-range pair energy correction. */ + double pair_energy_correction(std::size_t p1, std::size_t p2, auto &aosoa, + double q1q2) const { + double energy = 0.; + if (elc.dielectric_contrast_on) { + energy = std::visit( + [this, &aosoa, p1, p2, q1q2](auto &p3m_ptr) { + auto const &p3m = *p3m_ptr; + auto energy = 0.; + elc.dielectric_layers_contribution( + *m_box_geo, p1, p2, aosoa, q1q2, + [&](double q_eff, Utils::Vector3d const &d) { + energy += p3m.pair_energy(q_eff, d.norm()); + }); + elc.dielectric_layers_contribution( + *m_box_geo, p2, p1, aosoa, q1q2, + [&](double q_eff, Utils::Vector3d const &d) { + energy += p3m.pair_energy(q_eff, d.norm()); + }); + return energy / 2.; + }, + base_solver); + } + return energy; + } +#endif // SHARED_MEMORY_PARALLELISM + /** @brief Calculate short-range pair energy correction. */ double pair_energy_correction(Particle const &p1, Particle const &p2, double q1q2) const { @@ -261,12 +317,12 @@ struct ElectrostaticLayerCorrection auto const &p3m = *p3m_ptr; auto energy = 0.; elc.dielectric_layers_contribution( - p3m, *m_box_geo, pos1, pos2, q1q2, + *m_box_geo, pos1, pos2, q1q2, [&](double q_eff, Utils::Vector3d const &d) { energy += p3m.pair_energy(q_eff, d.norm()); }); elc.dielectric_layers_contribution( - p3m, *m_box_geo, pos2, pos1, q1q2, + *m_box_geo, pos2, pos1, q1q2, [&](double q_eff, Utils::Vector3d const &d) { energy += p3m.pair_energy(q_eff, d.norm()); }); @@ -287,12 +343,12 @@ struct ElectrostaticLayerCorrection auto const &pos2 = p2.pos(); auto const &p3m = *p3m_ptr; elc.dielectric_layers_contribution( - p3m, *m_box_geo, pos1, pos2, q1q2, + *m_box_geo, pos1, pos2, q1q2, [&](double q_eff, Utils::Vector3d const &d) { p1.force() += p3m.pair_force(q_eff, d, d.norm()); }); elc.dielectric_layers_contribution( - p3m, *m_box_geo, pos2, pos1, q1q2, + *m_box_geo, pos2, pos1, q1q2, [&](double q_eff, Utils::Vector3d const &d) { p2.force() += p3m.pair_force(q_eff, d, d.norm()); }); @@ -333,7 +389,7 @@ struct ElectrostaticLayerCorrection /// the energy calculation double calc_energy(ParticleRange const &particles) const; - template void visit_base_solver(Visitor &&visitor) const { + void visit_base_solver(auto &&visitor) const { std::visit(visitor, base_solver); } }; diff --git a/src/core/electrostatics/icc.cpp b/src/core/electrostatics/icc.cpp index fb30648c90b..0f63fd3d765 100644 --- a/src/core/electrostatics/icc.cpp +++ b/src/core/electrostatics/icc.cpp @@ -43,11 +43,17 @@ #include "electrostatics/p3m.hpp" #include "errorhandling.hpp" #include "integrators/Propagation.hpp" +#include "short_range_cabana.hpp" #include "system/System.hpp" #include #include +#ifdef SHARED_MEMORY_PARALLELISM +#include +#include +#endif + #include #include #include @@ -70,12 +76,12 @@ static void force_calc_icc( Coulomb::ShortRangeForceKernel::result_type const &coulomb_kernel, Coulomb::ShortRangeForceCorrectionsKernel::result_type const &elc_kernel) { // reset forces - for (auto &p : particles) { - p.force() = {}; - } - for (auto &p : ghost_particles) { - p.force() = {}; - } + auto const reset_kernel = [](Particle &p) { p.force_and_torque() = {}; }; + cell_structure.for_each_local_particle(reset_kernel); + cell_structure.for_each_ghost_particle(reset_kernel); +#ifdef SHARED_MEMORY_PARALLELISM + cell_structure.reset_local_force(); +#endif // calc ICC forces cell_structure.non_bonded_loop( @@ -96,10 +102,7 @@ static void force_calc_icc( }); } -void ICCStar::iteration(CellStructure &cell_structure, - ParticleRange const &particles, - ParticleRange const &ghost_particles) { - +void ICCStar::iteration() { try { sanity_check(); } catch (std::runtime_error const &err) { @@ -108,7 +111,10 @@ void ICCStar::iteration(CellStructure &cell_structure, } auto &system = get_system(); + auto &cell_structure = *system.cell_structure; auto const &coulomb = system.coulomb; + auto const particles = cell_structure.local_particles(); + auto const ghost_particles = cell_structure.ghost_particles(); auto const prefactor = std::visit( [](auto const &ptr) { return ptr->prefactor; }, *coulomb.impl->solver); auto const pref = 1. / (prefactor * 2. * std::numbers::pi); @@ -116,6 +122,12 @@ void ICCStar::iteration(CellStructure &cell_structure, auto const elc_kernel = coulomb.pair_force_elc_kernel(); icc_cfg.citeration = 0; +#ifdef SHARED_MEMORY_PARALLELISM + using execution_space = Kokkos::DefaultExecutionSpace; + auto const &unique_particles = cell_structure.get_unique_particles(); + auto const &local_force = cell_structure.get_local_force(); +#endif // SHARED_MEMORY_PARALLELISM + auto global_max_rel_diff = 0.; for (int j = 0; j < icc_cfg.max_iterations; j++) { @@ -126,6 +138,22 @@ void ICCStar::iteration(CellStructure &cell_structure, elc_kernel); system.coulomb.calc_long_range_force(particles); cell_structure.ghosts_reduce_forces(); +#ifdef SHARED_MEMORY_PARALLELISM + // force reduction + int num_threads = execution_space().concurrency(); + kokkos_parallel_range_for>( + "reduction", std::size_t{0}, unique_particles.size(), + [&local_force, &unique_particles, num_threads](std::size_t const i) { + Utils::Vector3d force{}; + 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); + } + unique_particles.at(i)->force() += force; + }); + Kokkos::fence(); +#endif // SHARED_MEMORY_PARALLELISM auto max_rel_diff = 0.; @@ -193,6 +221,10 @@ void ICCStar::iteration(CellStructure &cell_structure, /* Update charges on ghosts. */ cell_structure.ghosts_update(Cells::DATA_PART_PROPERTIES); +#ifdef SHARED_MEMORY_PARALLELISM + // refresh local properties + update_aosoa_charges(cell_structure); +#endif icc_cfg.citeration++; @@ -274,9 +306,8 @@ struct SanityChecksICC { void ICCStar::sanity_check() const { sanity_checks_active_solver(); #ifdef NPT - if (get_system().propagation->integ_switch == INTEG_METHOD_NPT_ISO_AND || - get_system().propagation->integ_switch == INTEG_METHOD_NPT_ISO_MTK) { - throw std::runtime_error("ICC does not work in the NPT ensemble"); + if (get_system().has_npt_enabled()) { + throw std::runtime_error("ICC does not work in the NpT ensemble"); } #endif } @@ -290,12 +321,17 @@ void ICCStar::sanity_checks_active_solver() const { } } +bool System::System::has_icc_enabled() const { + return coulomb.impl->extension and + std::holds_alternative>( + *coulomb.impl->extension); +} + void System::System::update_icc_particles() { if (coulomb.impl->extension) { if (auto icc = std::get_if>( get_ptr(coulomb.impl->extension))) { - (**icc).iteration(*cell_structure, cell_structure->local_particles(), - cell_structure->ghost_particles()); + (**icc).iteration(); } } } diff --git a/src/core/electrostatics/icc.hpp b/src/core/electrostatics/icc.hpp index 72bb48cd1c6..9831663929a 100644 --- a/src/core/electrostatics/icc.hpp +++ b/src/core/electrostatics/icc.hpp @@ -97,8 +97,7 @@ struct ICCStar : public System::Leaf { * The main iterative scheme, where the surface element charges are calculated * self-consistently. */ - void iteration(CellStructure &cell_structure, ParticleRange const &particles, - ParticleRange const &ghost_particles); + void iteration(); void on_activation() const; void sanity_checks_active_solver() const; diff --git a/src/core/electrostatics/p3m.cpp b/src/core/electrostatics/p3m.cpp index a0250247e3b..4c8109c149c 100644 --- a/src/core/electrostatics/p3m.cpp +++ b/src/core/electrostatics/p3m.cpp @@ -54,14 +54,17 @@ #include "ParticleRange.hpp" #include "PropagationMode.hpp" #include "actor/visitors.hpp" +#include "aosoa_pack.hpp" #include "cell_system/CellStructure.hpp" #include "cell_system/CellStructureType.hpp" +#include "cell_system/particle_enumeration.hpp" #include "communication.hpp" #include "errorhandling.hpp" #include "integrators/Propagation.hpp" #include "npt.hpp" #include "p3m/send_mesh.hpp" #include "particle_reduction.hpp" +#include "short_range_cabana.hpp" #include "system/GpuParticleData.hpp" #include "system/System.hpp" #include "tuning.hpp" @@ -72,10 +75,6 @@ #include #include -#ifdef SHARED_MEMORY_PARALLELISM -#include -#endif - #include #include #include @@ -83,6 +82,11 @@ #include #include +#ifdef SHARED_MEMORY_PARALLELISM +#include +#include +#endif + #include #include #include @@ -139,12 +143,11 @@ void CoulombP3MImpl::count_charged_particles() { double local_q = 0.0; double local_q2 = 0.0; }; - Reduction::AddPartialResultKernel kernel = [](Particle const &p, - Res &res) { + Reduction::AddPartialResultKernel kernel = [](Res &acc, auto const &p) { if (p.q() != 0.0) { - res.local_n++; - res.local_q2 += Utils::sqr(p.q()); - res.local_q += p.q(); + acc.local_n++; + acc.local_q2 += Utils::sqr(p.q()); + acc.local_q += p.q(); } }; @@ -367,22 +370,40 @@ template struct AssignCharge { this->operator()(p3m, q, weights); } - template - void operator()(auto &p3m, combined_ranges const &p_q_pos_range) { #ifdef SHARED_MEMORY_PARALLELISM - // multi-threading -> cache sizes must be equal to the number of particles - auto const include_neutral_particles = Kokkos::num_threads() > 1; -#else - auto constexpr include_neutral_particles = false; -#endif + void operator()(auto &p3m, auto &cell_structure) { + auto const &aosoa = cell_structure.get_aosoa(); + auto const n_part = cell_structure.count_local_particles(); + p3m.inter_weights.zfill(n_part); // allocate buffer for parallel write + kokkos_parallel_range_for( + "InterpolationWeights", std::size_t{0u}, n_part, [&](auto p_index) { + Utils::Vector3d const p_pos{aosoa.position(p_index, 0), + aosoa.position(p_index, 1), + aosoa.position(p_index, 2)}; + auto constexpr memory_order = + std::remove_reference::type::memory_order; + auto const weights = + p3m_calculate_interpolation_weights( + p_pos, p3m.params.ai, p3m.local_mesh); + p3m.inter_weights.store_at(p_index, weights); + }); + // serial part of charge assignment + for (std::size_t p_index{0}; p_index < n_part; ++p_index) { + auto const weights = p3m.inter_weights.template load(p_index); + this->operator()(p3m, aosoa.charge(p_index), weights); + } + } +#else // SHARED_MEMORY_PARALLELISM + void operator()(auto &p3m, auto const &p_q_pos_range) { for (auto zipped : p_q_pos_range) { auto const p_q = boost::get<0>(zipped); - auto const &p_pos = boost::get<1>(zipped); - if (include_neutral_particles or p_q != 0.0) { + if (p_q != 0.0) { + auto const &p_pos = boost::get<1>(zipped); this->operator()(p3m, p_q, p_pos, p3m.inter_weights); } } } +#endif // SHARED_MEMORY_PARALLELISM }; } // namespace @@ -391,11 +412,16 @@ void CoulombP3MImpl::charge_assign( ParticleRange const &particles) { prepare_fft_mesh(true); +#ifdef SHARED_MEMORY_PARALLELISM + Utils::integral_parameter( + p3m.params.cao, p3m, *get_system().cell_structure); +#else // SHARED_MEMORY_PARALLELISM auto p_q_range = ParticlePropertyRange::charge_range(particles); auto p_pos_range = ParticlePropertyRange::pos_range(particles); Utils::integral_parameter( p3m.params.cao, p3m, boost::combine(p_q_range, p_pos_range)); +#endif // SHARED_MEMORY_PARALLELISM } template @@ -410,53 +436,50 @@ void CoulombP3MImpl::assign_charge( } } +namespace { template struct AssignForces { - template - void operator()(auto &p3m, double force_prefac, - combined_ranges const &p_q_force_range) const { + void operator()(auto &p3m, auto force_prefac, +#ifdef SHARED_MEMORY_PARALLELISM + CellStructure &cell_structure +#else + auto const &p_q_force_range +#endif + ) const { assert(cao == p3m.inter_weights.cao()); - auto const kernel = [&p3m](double pref, auto &p_force, - std::size_t p_index) { - if (pref != 0.) { - auto const weights = p3m.inter_weights.template load(p_index); - - Utils::Vector3d force{}; - p3m_interpolate(p3m.local_mesh, weights, - [&force, &p3m](int ind, double w) { - force[0u] += w * double(p3m.rs_E_fields[0u][ind]); - force[1u] += w * double(p3m.rs_E_fields[1u][ind]); - force[2u] += w * double(p3m.rs_E_fields[2u][ind]); - }); + auto const kernel = [&p3m](auto pref, auto &p_force, std::size_t p_index) { + auto const weights = p3m.inter_weights.template load(p_index); - p_force -= pref * force; - } - }; + Utils::Vector3d force{}; + p3m_interpolate(p3m.local_mesh, weights, + [&force, &p3m](int ind, double w) { + force[0u] += w * double(p3m.rs_E_fields[0u][ind]); + force[1u] += w * double(p3m.rs_E_fields[1u][ind]); + force[2u] += w * double(p3m.rs_E_fields[2u][ind]); + }); #ifdef SHARED_MEMORY_PARALLELISM - if (Kokkos::num_threads() > 1) { - std::vector q_vals; - std::vector f_ptrs; - q_vals.reserve(p3m.inter_weights.size()); - f_ptrs.reserve(p3m.inter_weights.size()); - for (auto zipped : p_q_force_range) { - q_vals.emplace_back(boost::get<0>(zipped)); - f_ptrs.emplace_back(&boost::get<1>(zipped)); - } - Kokkos::RangePolicy<> policy(std::size_t{0u}, q_vals.size()); - auto const *q_vals_data = q_vals.data(); - auto const *f_ptrs_data = f_ptrs.data(); - Kokkos::parallel_for( - "AssignForces", policy, KOKKOS_LAMBDA(std::size_t p_index) { - auto p_q = q_vals_data[p_index]; - auto &p_force = *f_ptrs_data[p_index]; - kernel(p_q * force_prefac, p_force, p_index); - }); - return; - } + auto const thread_id = omp_get_thread_num(); + p_force(p_index, thread_id, 0) -= pref * force[0]; + p_force(p_index, thread_id, 1) -= pref * force[1]; + p_force(p_index, thread_id, 2) -= pref * force[2]; +#else + p_force -= pref * force; #endif + }; +#ifdef SHARED_MEMORY_PARALLELISM + auto const n_part = p3m.inter_weights.size(); + auto const &aosoa = cell_structure.get_aosoa(); + auto &local_force = cell_structure.get_local_force(); + kokkos_parallel_range_for( + "AssignForces", std::size_t{0u}, n_part, [&](std::size_t p_index) { + if (auto const pref = aosoa.charge(p_index) * force_prefac) { + kernel(pref, local_force, p_index); + } + }); +#else // SHARED_MEMORY_PARALLELISM /* charged particle counter */ std::size_t p_index{0ul}; @@ -468,12 +491,25 @@ template struct AssignForces { ++p_index; } } +#endif // SHARED_MEMORY_PARALLELISM } }; +} // namespace -template +#ifdef SHARED_MEMORY_PARALLELISM +static auto calc_dipole_moment(boost::mpi::communicator const &comm, + auto const &cs, auto const &box_geo) { + auto const local_dip = reduce_over_local_particles( + cs, + [&box_geo](Utils::Vector3d &acc, Particle const &p) { + acc += p.q() * box_geo.unfolded_position(p.pos(), p.image_box()); + }, + [](Utils::Vector3d &a, Utils::Vector3d const &b) { a = a + b; }); + return boost::mpi::all_reduce(comm, local_dip, std::plus<>()); +} +#else // SHARED_MEMORY_PARALLELISM static auto calc_dipole_moment(boost::mpi::communicator const &comm, - combined_ranges const &p_q_unfolded_pos_range) { + auto const &p_q_unfolded_pos_range) { auto const local_dip = boost::accumulate(p_q_unfolded_pos_range, Utils::Vector3d{}, [](Utils::Vector3d const &dip, auto const &q_pos) { @@ -483,6 +519,7 @@ static auto calc_dipole_moment(boost::mpi::communicator const &comm, }); return boost::mpi::all_reduce(comm, local_dip, std::plus<>()); } +#endif // SHARED_MEMORY_PARALLELISM template void CoulombP3MImpl::kernel_ks_charge_density() { @@ -636,10 +673,7 @@ double CoulombP3MImpl::long_range_kernel( auto const &system = get_system(); auto const &box_geo = *system.box_geo; #ifdef NPT - auto const npt_flag = - force_flag and - ((system.propagation->integ_switch == INTEG_METHOD_NPT_ISO_AND) or - (system.propagation->integ_switch == INTEG_METHOD_NPT_ISO_MTK)); + auto const npt_flag = force_flag and system.has_npt_enabled(); #else auto constexpr npt_flag = false; #endif @@ -654,18 +688,28 @@ double CoulombP3MImpl::long_range_kernel( kernel_ks_charge_density(); - // Calculate the dipole term +#ifdef SHARED_MEMORY_PARALLELISM + auto &cell_structure = *system.cell_structure; + auto const &local_force = cell_structure.get_local_force(); + auto const &aosoa = cell_structure.get_aosoa(); +#else auto p_q_range = ParticlePropertyRange::charge_range(particles); auto p_force_range = ParticlePropertyRange::force_range(particles); auto p_unfolded_pos_range = ParticlePropertyRange::unfolded_pos_range(particles, box_geo); +#endif // SHARED_MEMORY_PARALLELISM // The dipole moment is only needed if we don't have metallic boundaries - auto const box_dipole = - (p3m.params.epsilon != P3M_EPSILON_METALLIC) - ? std::make_optional(calc_dipole_moment( - comm_cart, boost::combine(p_q_range, p_unfolded_pos_range))) - : std::nullopt; + auto const box_dipole = (p3m.params.epsilon != P3M_EPSILON_METALLIC) + ? std::make_optional(calc_dipole_moment( +#ifdef SHARED_MEMORY_PARALLELISM + comm_cart, cell_structure, box_geo)) +#else + comm_cart, + boost::combine(p_q_range, + p_unfolded_pos_range))) +#endif + : std::nullopt; auto const volume = box_geo.volume(); auto const pref = 4. * std::numbers::pi / volume / (2. * p3m.params.epsilon + 1.); @@ -677,19 +721,36 @@ double CoulombP3MImpl::long_range_kernel( // assign particle forces auto const force_prefac = prefactor / volume; +#ifdef SHARED_MEMORY_PARALLELISM + auto &particle_data = cell_structure; +#else + auto const particle_data = boost::combine(p_q_range, p_force_range); +#endif Utils::integral_parameter( - p3m.params.cao, p3m, force_prefac, - boost::combine(p_q_range, p_force_range)); + p3m.params.cao, p3m, force_prefac, particle_data); // add dipole forces // Eq. (3.19) @cite deserno00b if (box_dipole) { auto const dm = prefactor * pref * box_dipole.value(); +#ifdef SHARED_MEMORY_PARALLELISM + auto const n_part = cell_structure.count_local_particles(); + kokkos_parallel_range_for( + "AssignForcesBoxDipole", std::size_t{0u}, n_part, + [&aosoa, &local_force, dm](auto p_index) { + auto const thread_id = omp_get_thread_num(); + auto const q = aosoa.charge(p_index); + local_force(p_index, thread_id, 0) -= q * dm[0]; + local_force(p_index, thread_id, 1) -= q * dm[1]; + local_force(p_index, thread_id, 2) -= q * dm[2]; + }); +#else // SHARED_MEMORY_PARALLELISM for (auto zipped : boost::combine(p_q_range, p_force_range)) { auto p_q = boost::get<0>(zipped); auto &p_force = boost::get<1>(zipped); p_force -= p_q * dm; } +#endif // SHARED_MEMORY_PARALLELISM } } @@ -1065,8 +1126,7 @@ void CoulombP3MImpl::add_long_range_forces_gpu( ParticleRange const &particles) { if constexpr (Architecture == Arch::GPU) { #ifdef NPT - if ((get_system().propagation->integ_switch == INTEG_METHOD_NPT_ISO_AND) or - (get_system().propagation->integ_switch == INTEG_METHOD_NPT_ISO_MTK)) { + if (get_system().has_npt_enabled()) { get_system().npt_add_virial_contribution(long_range_energy(particles)); } #else diff --git a/src/core/electrostatics/p3m.hpp b/src/core/electrostatics/p3m.hpp index edb4e6a3222..3bad1a5e695 100644 --- a/src/core/electrostatics/p3m.hpp +++ b/src/core/electrostatics/p3m.hpp @@ -145,7 +145,7 @@ struct CoulombP3M : public Coulomb::Actor { /** Calculate real-space contribution of p3m Coulomb pair forces. */ Utils::Vector3d pair_force(double q1q2, Utils::Vector3d const &d, double dist) const { - if ((q1q2 == 0.) || dist >= p3m_params.r_cut || dist <= 0.) { + if (q1q2 == 0. or dist >= p3m_params.r_cut or dist <= 0.) { return {}; } auto const alpha = p3m_params.alpha; @@ -166,7 +166,7 @@ struct CoulombP3M : public Coulomb::Actor { /** Calculate real-space contribution of Coulomb pair energy. */ // Eq. (3.6) @cite deserno00b double pair_energy(double q1q2, double dist) const { - if ((q1q2 == 0.) || dist >= p3m_params.r_cut || dist <= 0.) { + if (q1q2 == 0. or dist >= p3m_params.r_cut or dist <= 0.) { return {}; } auto const adist = p3m_params.alpha * dist; diff --git a/src/core/electrostatics/solver.hpp b/src/core/electrostatics/solver.hpp index 207a5f30e39..56a59712bc4 100644 --- a/src/core/electrostatics/solver.hpp +++ b/src/core/electrostatics/solver.hpp @@ -19,7 +19,7 @@ #pragma once -#include "config/config.hpp" +#include #include "actor/optional.hpp" #include "actor/traits.hpp" @@ -63,6 +63,7 @@ struct Solver { Solver(); #else // ELECTROSTATICS Solver() = default; + constexpr double cutoff() const { return INACTIVE_CUTOFF; } #endif // ELECTROSTATICS using ShortRangeForceKernel = diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 52791093212..957b37b92b3 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -24,6 +24,8 @@ * The corresponding header file is forces.hpp. */ +#include + #include "BoxGeometry.hpp" #include "Particle.hpp" #include "ParticleRange.hpp" @@ -45,6 +47,7 @@ #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "npt.hpp" #include "rotation.hpp" +#include "short_range_cabana.hpp" #include "short_range_loop.hpp" #include "system/System.hpp" #include "thermostat.hpp" @@ -59,7 +62,6 @@ #endif #ifdef SHARED_MEMORY_PARALLELISM -#include "short_range_cabana.hpp" #include #endif @@ -92,12 +94,12 @@ static ParticleForce external_force(Particle const &p) { } /** Combined force initialization and Langevin noise application */ -void init_forces_and_thermostat(const CellStructure &cell_structure, - System::System &system) { +void init_forces_and_thermostat(System::System const &system) { #ifdef CALIPER CALI_CXX_MARK_FUNCTION; #endif + auto &cell_structure = *system.cell_structure; auto const &propagation = *system.propagation; auto const &thermostat = *system.thermostat; auto const kT = thermostat.kT; @@ -111,11 +113,10 @@ void init_forces_and_thermostat(const CellStructure &cell_structure, // Single pass over all local particles cell_structure.for_each_local_particle([&](Particle &p) { - // Initialize force with external forces (original init_forces logic) + // Initialize force with external forces p.force_and_torque() = external_force(p); - // Apply Langevin noise if thermostat is active (original - // thermostat_force_init logic) + // Apply Langevin noise if thermostat is active if (langevin_active) { auto const &langevin = *thermostat.langevin; if (propagation.should_propagate_with(p, PropagationMode::TRANS_LANGEVIN)) @@ -127,6 +128,9 @@ void init_forces_and_thermostat(const CellStructure &cell_structure, #endif } }); +#ifdef SHARED_MEMORY_PARALLELISM + cell_structure.reset_local_force(); +#endif // Initialize ghost forces (unchanged) init_forces_ghosts(cell_structure); @@ -166,18 +170,12 @@ void System::System::calculate_forces() { #ifdef COLLISION_DETECTION collision_detection->clear_queue(); + auto const collision_detection_cutoff = collision_detection->cutoff(); +#else + auto const collision_detection_cutoff = INACTIVE_CUTOFF; #endif bond_breakage->clear_queue(); auto particles = cell_structure->local_particles(); -#ifdef ELECTROSTATICS - if (coulomb.impl->extension) { - if (auto icc = std::get_if>( - get_ptr(coulomb.impl->extension))) { - auto ghost_particles = cell_structure->ghost_particles(); - (**icc).iteration(*cell_structure, particles, ghost_particles); - } - } -#endif // ELECTROSTATICS #ifdef NPT if (propagation->used_propagations & PropagationMode::TRANS_LANGEVIN_NPT) { // reset virial part of instantaneous pressure @@ -185,56 +183,53 @@ void System::System::calculate_forces() { } #endif // Use combined function instead of two separate calls - init_forces_and_thermostat(*cell_structure, *this); - - calc_long_range_forces(particles); auto const elc_kernel = coulomb.pair_force_elc_kernel(); auto const coulomb_kernel = coulomb.pair_force_kernel(); auto const dipoles_kernel = dipoles.pair_force_kernel(); auto const coulomb_u_kernel = coulomb.pair_energy_kernel(); - -#ifdef ELECTROSTATICS - auto const coulomb_cutoff = coulomb.cutoff(); -#else - auto const coulomb_cutoff = INACTIVE_CUTOFF; -#endif - -#ifdef DIPOLES - auto const dipole_cutoff = dipoles.cutoff(); -#else - auto const dipole_cutoff = INACTIVE_CUTOFF; -#endif -#ifdef COLLISION_DETECTION - auto const collision_detection_cutoff = collision_detection->cutoff(); -#else - auto const collision_detection_cutoff = INACTIVE_CUTOFF; -#endif + auto *const virial = get_npt_virial(); // interaction kernel is defined auto bond_kernel = [coulomb_kernel_ptr = get_ptr(coulomb_kernel), &bonded_ias = *bonded_ias, - &bond_breakage = *bond_breakage, + &bond_breakage = *bond_breakage, virial, &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); + box_geo, virial, coulomb_kernel_ptr); }; VerletCriterion<> const verlet_criterion{*this, cell_structure->get_verlet_skin(), get_interaction_range(), - coulomb_cutoff, - dipole_cutoff, + coulomb.cutoff(), + dipoles.cutoff(), collision_detection_cutoff}; #ifdef SHARED_MEMORY_PARALLELISM #ifdef CALIPER - CALI_MARK_BEGIN("parallel short range"); + CALI_MARK_BEGIN("convert particles AoS to SoA"); #endif - using execution_space = Kokkos::DefaultExecutionSpace; update_cabana_state(*cell_structure, verlet_criterion, get_interaction_range()); +#ifdef CALIPER + CALI_MARK_END("convert particles AoS to SoA"); +#endif +#endif +#ifdef ELECTROSTATICS + if (coulomb.impl->extension) { + update_icc_particles(); + } +#endif // ELECTROSTATICS + init_forces_and_thermostat(*this); + calc_long_range_forces(particles); + +#ifdef SHARED_MEMORY_PARALLELISM +#ifdef CALIPER + CALI_MARK_BEGIN("parallel short range"); +#endif + using execution_space = Kokkos::DefaultExecutionSpace; auto const &unique_particles = cell_structure->get_unique_particles(); auto const &local_force = cell_structure->get_local_force(); #ifdef ROTATION @@ -247,20 +242,13 @@ void System::System::calculate_forces() { 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, + *thermostat, *box_geo, unique_particles, local_force, #ifdef ROTATION local_torque, #endif #ifdef NPT - local_virial, + virial, local_virial, #endif aosoa); @@ -268,13 +256,14 @@ void System::System::calculate_forces() { 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::RangePolicy policy(std::size_t{0}, + unique_particles.size()); Kokkos::parallel_for("reduction", policy, [&local_force, #ifdef ROTATION &local_torque, #endif - &unique_particles, num_threads](const int i) { + &unique_particles, num_threads](std::size_t const i) { Utils::Vector3d force{}; #ifdef ROTATION Utils::Vector3d torque{}; @@ -297,13 +286,13 @@ void System::System::calculate_forces() { 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); + if (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 @@ -329,6 +318,7 @@ void System::System::calculate_forces() { coulomb_u_kernel_ptr = get_ptr(coulomb_u_kernel), &nonbonded_ias = *nonbonded_ias, &thermostat = *thermostat, &bonded_ias = *bonded_ias, + virial, #ifdef COLLISION_DETECTION &collision_detection = *collision_detection, #endif @@ -337,8 +327,8 @@ void System::System::calculate_forces() { 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); + thermostat, box_geo, bonded_ias, virial, 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); @@ -350,6 +340,7 @@ void System::System::calculate_forces() { bonded_ias->maximal_cutoff(), verlet_criterion); #endif // SHARED_MEMORY_PARALLELISM + constraints->add_forces(particles, get_sim_time()); oif_global->calculate_forces(); @@ -406,16 +397,3 @@ void calc_long_range_forces(const ParticleRange &particles) { Dipoles::get_dipoles().calc_long_range_force(particles); #endif // DIPOLES } - -#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); -} -void npt_add_virial_diagonalSum_contribution(double diagonal_sum) { - ::System::get_system().npt_add_virial_contribution(diagonal_sum); -} -#endif diff --git a/src/core/forces.hpp b/src/core/forces.hpp index 46e1d0ee27e..be223854f61 100644 --- a/src/core/forces.hpp +++ b/src/core/forces.hpp @@ -27,24 +27,17 @@ * Implementation in forces.cpp. */ -#include "ParticleRange.hpp" +#include -#include +#include "ParticleRange.hpp" +#include "cell_system/CellStructure.hpp" +#include "system/System.hpp" /** Combined force initialization and Langevin noise application. */ -void init_forces_and_thermostat(const CellStructure &cell_structure, - class System::System &system); +void init_forces_and_thermostat(System::System const &system); /** Set forces of all ghosts to zero */ void init_forces_ghosts(const CellStructure &cell_structure); /** Calculate long range forces (P3M, ...). */ 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); -#endif diff --git a/src/core/forces_cabana.hpp b/src/core/forces_cabana.hpp index 1b3d24a4158..d35637311e1 100644 --- a/src/core/forces_cabana.hpp +++ b/src/core/forces_cabana.hpp @@ -26,8 +26,12 @@ #include "aosoa_pack.hpp" #include "forces_inline.hpp" +#include + #include +#include + #include #if defined(__GNUG__) or defined(__clang__) @@ -40,21 +44,18 @@ 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 + Utils::Vector3d *const global_virial; CellStructure::VirialType const &local_virial; #endif CellStructure::AoSoA_pack const &aosoa; @@ -63,46 +64,36 @@ struct 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) + Thermostat::Thermostat const &thermostat_, BoxGeometry const &box_geo_, std::vector const &unique_particles_, -#endif CellStructure::ForceType const &local_force_, #ifdef ROTATION CellStructure::ForceType const &local_torque_, #endif #ifdef NPT + Utils::Vector3d *const global_virial_, 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_), + coulomb_kernel(coulomb_kernel_), dipoles_kernel(dipoles_kernel_), + elc_kernel(elc_kernel_), coulomb_u_kernel(coulomb_u_kernel_), + thermostat(thermostat_), box_geo(box_geo_), + unique_particles(unique_particles_), local_force(local_force_), #ifdef ROTATION local_torque(local_torque_), #endif #ifdef NPT - local_virial(local_virial_), + global_virial(global_virial_), local_virial(local_virial_), #endif aosoa(aosoa_) { } ESPRESSO_ATTR_ALWAYS_INLINE KOKKOS_INLINE_FUNCTION void - operator()(int i, int j) const { + operator()(std::size_t i, std::size_t j) const { auto const thread_id = omp_get_thread_num(); @@ -110,25 +101,25 @@ struct ForcesKernel { nonbonded_ias.get_ia_param(aosoa.type(i), aosoa.type(j)); ParticleForce pf{}; + #ifdef NPT Utils::Vector3d virial{}; + auto *const virial_handle = global_virial ? &virial : nullptr; +#else + Utils::Vector3d *const virial_handle = nullptr; #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) { @@ -141,21 +132,32 @@ struct ForcesKernel { #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 + add_non_bonded_pair_force_with_p( + p1, p2, pf, d, dist, dist * dist, q1q2, ia_params, do_nonbonded_flag, + thermostat, box_geo, bonded_ias, virial_handle, coulomb_kernel, + dipoles_kernel, elc_kernel, coulomb_u_kernel); + +#ifdef ELECTROSTATICS + // real-space electrostatic charge-charge interaction + if (q1q2 != 0. and coulomb_kernel != nullptr) { + pf.f += (*coulomb_kernel)(q1q2, d, dist); + } +#endif // ELECTROSTATICS +#ifdef DIPOLES + // real-space magnetic dipole-dipole interaction + if (dipoles_kernel) { + auto const d1d2 = p1.dipm() * p2.dipm(); + if (d1d2 != 0.) { + pf += (*dipoles_kernel)(d1d2, p1.calc_dip(), p2.calc_dip(), d, dist, + dist * dist); + } + } +#endif // DIPOLES local_force(i, thread_id, 0) += pf.f[0]; local_force(i, thread_id, 1) += pf.f[1]; @@ -176,9 +178,11 @@ struct ForcesKernel { 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]; + if (virial_handle) { + local_virial(thread_id, 0) += virial[0]; + local_virial(thread_id, 1) += virial[1]; + local_virial(thread_id, 2) += virial[2]; + } #endif } }; diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index 381a1934347..5fec6a0ef28 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -25,14 +25,10 @@ * Force calculation. */ -#include "config/config.hpp" +#include #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" @@ -211,14 +207,12 @@ inline void add_non_bonded_pair_without_p( * @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, + Particle &p1, Particle &p2, ParticleForce &pf, Utils::Vector3d const &d, + double dist, double dist2, double q1q2, IA_parameters const &ia_params, + [[maybe_unused]] bool do_nonbonded_flag, Thermostat::Thermostat const &thermostat, BoxGeometry const &box_geo, [[maybe_unused]] BondedInteractionsMap const &bonded_ias, + [[maybe_unused]] Utils::Vector3d *const virial, Coulomb::ShortRangeForceKernel::kernel_type const *coulomb_kernel, Dipoles::ShortRangeForceKernel::kernel_type const *dipoles_kernel, Coulomb::ShortRangeForceCorrectionsKernel::kernel_type const *elc_kernel, @@ -230,7 +224,7 @@ inline void add_non_bonded_pair_force_with_p( if (dist < ia_params.max_cut) { #ifdef EXCLUSIONS - if (do_nonbonded) { + if (do_nonbonded_flag) { #endif #ifdef THOLE pf.f += thole_pair_force(p1, p2, ia_params, d, dist, bonded_ias, @@ -248,12 +242,10 @@ inline void add_non_bonded_pair_force_with_p( /* 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 + if (virial) { + *virial += hadamard_product(pf.f, d); + } +#endif // NPT /***********************************************/ /* short-range electrostatics */ @@ -262,19 +254,17 @@ inline void add_non_bonded_pair_force_with_p( #ifdef ELECTROSTATICS // real-space electrostatic charge-charge interaction if (q1q2 != 0. and coulomb_kernel != nullptr) { +#if not defined(SHARED_MEMORY_PARALLELISM) pf.f += (*coulomb_kernel)(q1q2, d, dist); +#endif // not SHARED_MEMORY_PARALLELISM #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 // SHARED_MEMORY_PARALLELISM + if (virial) { + (*virial)[0] += (*coulomb_u_kernel)(p1, p2, q1q2, d, dist); + } #endif // NPT -#ifdef P3M - if (elc_kernel) + if (elc_kernel) { (*elc_kernel)(p1, p2, q1q2); -#endif // P3M + } } #endif // ELECTROSTATICS @@ -297,10 +287,16 @@ inline void add_non_bonded_pair_force_with_p( /***********************************************/ #ifdef DIPOLES +#if not defined(SHARED_MEMORY_PARALLELISM) // real-space magnetic dipole-dipole if (dipoles_kernel) { - pf += (*dipoles_kernel)(p1, p2, d, dist, dist2); + auto const d1d2 = p1.dipm() * p2.dipm(); + if (d1d2 != 0.) { + pf += + (*dipoles_kernel)(d1d2, p1.calc_dip(), p2.calc_dip(), d, dist, dist2); + } } +#endif // not SHARED_MEMORY_PARALLELISM #endif } @@ -316,34 +312,32 @@ inline void add_non_bonded_pair_force_with_p( * @param[in] thermostat thermostat. * @param[in] box_geo box geometry. * @param[in] bonded_ias bonded interaction kernels. + * @param[out] virial NpT virial. * @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( +inline void 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, + [[maybe_unused]] Utils::Vector3d *const virial, 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 +#ifndef SHARED_MEMORY_PARALLELISM if (dist < ia_params.max_cut) { #ifdef EXCLUSIONS if (do_nonbonded_flag) { @@ -353,26 +347,18 @@ inline auto add_non_bonded_pair_force( } #endif } +#endif // not SHARED_MEMORY_PARALLELISM -#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 + p1, p2, pf, d, dist, dist2, q1q2, ia_params, do_nonbonded_flag, + thermostat, box_geo, bonded_ias, virial, coulomb_kernel, dipoles_kernel, + elc_kernel, coulomb_u_kernel); /***********************************************/ /* 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 +#ifndef SHARED_MEMORY_PARALLELISM p1.force_and_torque() += pf; p2.force_and_torque() += calc_opposing_force(pf, d); #endif @@ -424,8 +410,8 @@ inline std::optional calc_bond_pair_force( } inline bool add_bonded_two_body_force( - Bonded_IA_Parameters const &iaparams, Particle &p1, Particle &p2, - BoxGeometry const &box_geo, + Bonded_IA_Parameters const &iaparams, BoxGeometry const &box_geo, + Particle &p1, Particle &p2, [[maybe_unused]] Utils::Vector3d *const virial, Coulomb::ShortRangeForceKernel::kernel_type const *kernel) { auto const dx = box_geo.get_mi_vector(p1.pos(), p2.pos()); @@ -446,7 +432,9 @@ inline bool add_bonded_two_body_force( p2.force() -= result.value(); #ifdef NPT - npt_add_virial_force_contribution(result.value(), dx); + if (virial) { + *virial += hadamard_product(result.value(), dx); + } #endif return false; } @@ -554,6 +542,7 @@ add_bonded_force(Particle &p1, int bond_id, std::span partners, BondedInteractionsMap const &bonded_ia_params, BondBreakage::BondBreakage &bond_breakage, BoxGeometry const &box_geo, + [[maybe_unused]] Utils::Vector3d *const virial, Coulomb::ShortRangeForceKernel::kernel_type const *kernel) { // Consider for bond breakage @@ -579,8 +568,8 @@ add_bonded_force(Particle &p1, int bond_id, std::span partners, case 0: return false; case 1: - return add_bonded_two_body_force(iaparams, p1, *partners[0], box_geo, - kernel); + return add_bonded_two_body_force(iaparams, box_geo, p1, *partners[0], + virial, kernel); case 2: return add_bonded_three_body_force(iaparams, box_geo, p1, *partners[0], *partners[1]); diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 8106d42196e..abcab32d980 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -551,8 +551,7 @@ int System::System::integrate(int n_steps, int reuse_forces) { } #ifdef NPT - if ((propagation.integ_switch != INTEG_METHOD_NPT_ISO_AND) && - (propagation.integ_switch != INTEG_METHOD_NPT_ISO_MTK)) + if (not has_npt_enabled()) #endif { resort_particles_if_needed(*this); @@ -571,8 +570,7 @@ int System::System::integrate(int n_steps, int reuse_forces) { #ifdef VIRTUAL_SITES_RELATIVE if (has_vs_rel()) { #ifdef NPT - if ((propagation.integ_switch == INTEG_METHOD_NPT_ISO_AND) or - (propagation.integ_switch == INTEG_METHOD_NPT_ISO_MTK)) { + if (has_npt_enabled()) { cell_structure->update_ghosts_and_resort_particle( Cells::DATA_PART_PROPERTIES); } @@ -708,8 +706,7 @@ int System::System::integrate(int n_steps, int reuse_forces) { cell_structure->update_verlet_stats(n_steps, n_verlet_updates); #ifdef NPT - if ((propagation.integ_switch == INTEG_METHOD_NPT_ISO_AND) or - (propagation.integ_switch == INTEG_METHOD_NPT_ISO_MTK)) { + if (has_npt_enabled()) { synchronize_npt_state(); } #endif diff --git a/src/core/integrators/velocity_verlet_npt_Andersen.cpp b/src/core/integrators/velocity_verlet_npt_Andersen.cpp index bab5736df32..1de0f20e684 100644 --- a/src/core/integrators/velocity_verlet_npt_Andersen.cpp +++ b/src/core/integrators/velocity_verlet_npt_Andersen.cpp @@ -161,7 +161,7 @@ static void velocity_verlet_npt_propagate_AVOVA_And( } boost::mpi::broadcast(::comm_cart, scal, 0); - /* stochastic reserviors for velocities; + /* stochastic reservoirs for velocities; * @f$ p(t+0.5*dt) = p(t+0.5*dt) \exp(- \gamma_0 dt / m) * + \sqrt{k_B T (1 - \exp(-2 \gamma_0 dt)}N(0,1) @f$ * and 2nd propagate positions with dt/2 while rescaling positions velocities diff --git a/src/core/integrators/velocity_verlet_npt_MTK.cpp b/src/core/integrators/velocity_verlet_npt_MTK.cpp index 39b9d25ff19..e28ba3d23fc 100644 --- a/src/core/integrators/velocity_verlet_npt_MTK.cpp +++ b/src/core/integrators/velocity_verlet_npt_MTK.cpp @@ -135,7 +135,7 @@ static void velocity_verlet_npt_propagate_AVOVA_MTK( velocity_verlet_npt_propagate_pos_MTK(nptiso, particles, time_step); velocity_verlet_npt_propagate_pos(particles, time_step); - /* stochastic reserviors for conjugate momentum for particles + /* stochastic reservoirs for conjugate momentum for particles * @f$ p(t+0.5*dt) = p(t+0.5*dt) \exp(- \gamma_0 dt / m) * + \sqrt{k_B T (1 - \exp(-2 \gamma_0 dt)}N(0,1) @f$ */ @@ -153,7 +153,7 @@ static void velocity_verlet_npt_propagate_AVOVA_MTK( /* 1st propagatation Volume with dt/2 * @f$ V(t+0.5*dt) = \exp(0.5 * dt * 3 * p_{\epsilon} / W) * V(t) @f$, - * stochastic reserviors for conjugate momentum for V + * stochastic reservoirs for conjugate momentum for V * @f$ p_{\epsilon} = p_{epsilon}(t) \exp(- \gamma_V dt / W) * + \sqrt{k_B T (1 - \exp(-2 \gamma_V dt / W)}N(0,1) @f$, * 2nd propagatation Volume with dt/2 diff --git a/src/core/magnetostatics/dipoles.cpp b/src/core/magnetostatics/dipoles.cpp index f3d4c4a1417..7e1818c8640 100644 --- a/src/core/magnetostatics/dipoles.cpp +++ b/src/core/magnetostatics/dipoles.cpp @@ -95,7 +95,7 @@ double Solver::cutoff() const { } } #endif - return -1.; + return INACTIVE_CUTOFF; } void Solver::on_observable_calc() { diff --git a/src/core/magnetostatics/dipoles_inline.hpp b/src/core/magnetostatics/dipoles_inline.hpp index a3184dc17b9..7aba5bc3264 100644 --- a/src/core/magnetostatics/dipoles_inline.hpp +++ b/src/core/magnetostatics/dipoles_inline.hpp @@ -52,10 +52,11 @@ struct ShortRangeForceKernel { #ifdef DP3M result_type operator()(std::shared_ptr const &ptr) const { auto const &actor = *ptr; - return kernel_type{[&actor](Particle const &p1, Particle const &p2, + return kernel_type{[&actor](double d1d2, Utils::Vector3d const &dip1, + Utils::Vector3d const &dip2, Utils::Vector3d const &d, double dist, double dist2) { - return actor.pair_force(p1, p2, d, dist2, dist); + return actor.pair_force(d1d2, dip1, dip2, d, dist, dist2); }}; } #endif // DP3M @@ -84,7 +85,7 @@ struct ShortRangeEnergyKernel { return kernel_type{[&actor](Particle const &p1, Particle const &p2, Utils::Vector3d const &d, double dist, double dist2) { - return actor.pair_energy(p1, p2, d, dist2, dist); + return actor.pair_energy(p1, p2, d, dist, dist2); }}; } #endif // DP3M diff --git a/src/core/magnetostatics/dp3m.cpp b/src/core/magnetostatics/dp3m.cpp index 0a6a83d4a3d..a9bda43f80c 100644 --- a/src/core/magnetostatics/dp3m.cpp +++ b/src/core/magnetostatics/dp3m.cpp @@ -51,6 +51,7 @@ #include "errorhandling.hpp" #include "integrators/Propagation.hpp" #include "npt.hpp" +#include "short_range_cabana.hpp" #include "system/System.hpp" #include "tuning.hpp" @@ -62,6 +63,11 @@ #include #include +#ifdef SHARED_MEMORY_PARALLELISM +#include +#include +#endif + #include #include #include @@ -157,6 +163,38 @@ void DipolarP3MImpl::init_cpu_kernels() { namespace { template struct AssignDipole { +#ifdef SHARED_MEMORY_PARALLELISM + void operator()(auto &dp3m, auto &cell_structure) { + using value_type = + typename std::remove_reference_t::value_type; + auto constexpr memory_order = Utils::MemoryOrder::ROW_MAJOR; + auto const &aosoa = cell_structure.get_aosoa(); + auto const &unique_particles = cell_structure.get_unique_particles(); + auto const n_part = cell_structure.count_local_particles(); + dp3m.inter_weights.zfill(n_part); // allocate buffer for parallel write + kokkos_parallel_range_for( + "InterpolationWeights", std::size_t{0u}, n_part, [&](auto p_index) { + Utils::Vector3d const p_pos{aosoa.position(p_index, 0), + aosoa.position(p_index, 1), + aosoa.position(p_index, 2)}; + auto const weights = + p3m_calculate_interpolation_weights( + p_pos, dp3m.params.ai, dp3m.local_mesh); + dp3m.inter_weights.store_at(p_index, weights); + }); + // serial part of dipole assignment + for (std::size_t p_index{0}; p_index < n_part; ++p_index) { + auto const weights = dp3m.inter_weights.template load(p_index); + auto const dip = unique_particles.at(p_index)->calc_dip(); + p3m_interpolate( + dp3m.local_mesh, weights, [&dip, &dp3m](int ind, double w) { + dp3m.mesh.rs_fields[0u][ind] += value_type(w * dip[0u]); + dp3m.mesh.rs_fields[1u][ind] += value_type(w * dip[1u]); + dp3m.mesh.rs_fields[2u][ind] += value_type(w * dip[2u]); + }); + } + } +#else // SHARED_MEMORY_PARALLELISM void operator()(auto &dp3m, Utils::Vector3d const &real_pos, Utils::Vector3d const &dip) const { using value_type = @@ -173,6 +211,7 @@ template struct AssignDipole { dp3m.inter_weights.template store(weights); } +#endif // SHARED_MEMORY_PARALLELISM }; } // namespace @@ -182,66 +221,132 @@ void DipolarP3MImpl::dipole_assign( dp3m.inter_weights.reset(dp3m.params.cao); /* prepare local FFT mesh */ - for (auto &rs_mesh_field : dp3m.mesh.rs_fields) - for (int j = 0; j < dp3m.local_mesh.size; j++) - rs_mesh_field[j] = 0.; + for (auto &rs_mesh_field : dp3m.mesh.rs_fields) { + std::ranges::fill_n(rs_mesh_field.data(), dp3m.local_mesh.size, 0.); + } +#ifdef SHARED_MEMORY_PARALLELISM + Utils::integral_parameter( + dp3m.params.cao, dp3m, *get_system().cell_structure); +#else // SHARED_MEMORY_PARALLELISM for (auto const &p : particles) { if (p.dipm() != 0.) { Utils::integral_parameter(dp3m.params.cao, dp3m, p.pos(), p.calc_dip()); } } +#endif // SHARED_MEMORY_PARALLELISM } namespace { template struct AssignTorques { void operator()(auto &dp3m, double prefac, int d_rs, - ParticleRange const &particles) const { - +#ifdef SHARED_MEMORY_PARALLELISM + CellStructure &cell_structure +#else + ParticleRange const &particles +#endif + ) const { + + assert(cao == dp3m.inter_weights.cao()); + + auto const kernel = [d_rs, &dp3m](auto const &pref, auto &p_torque, + std::size_t p_index) { + auto const weights = dp3m.inter_weights.template load(p_index); + Utils::Vector3d E{}; + p3m_interpolate(dp3m.local_mesh, weights, + [&E, &dp3m, d_rs](int ind, double w) { + E[d_rs] += w * double(dp3m.mesh.rs_scalar[ind]); + }); + + auto const torque = vector_product(pref, E); +#ifdef SHARED_MEMORY_PARALLELISM + auto const thread_id = omp_get_thread_num(); + p_torque(p_index, thread_id, 0) -= torque[0]; + p_torque(p_index, thread_id, 1) -= torque[1]; + p_torque(p_index, thread_id, 2) -= torque[2]; +#else + p_torque -= torque; +#endif + }; + +#ifdef SHARED_MEMORY_PARALLELISM + auto const n_part = dp3m.inter_weights.size(); + auto const &unique_particles = cell_structure.get_unique_particles(); + auto &local_torque = cell_structure.get_local_torque(); + kokkos_parallel_range_for( + "AssignTorques", std::size_t{0u}, n_part, [&](std::size_t p_index) { + auto const &p = *unique_particles.at(p_index); + if (p.dipm() != 0.) { + kernel(p.calc_dip() * prefac, local_torque, p_index); + } + }); +#else // SHARED_MEMORY_PARALLELISM /* magnetic particle index */ auto p_index = std::size_t{0ul}; for (auto &p : particles) { if (p.dipm() != 0.) { - auto const weights = dp3m.inter_weights.template load(p_index); - - Utils::Vector3d E{}; - p3m_interpolate(dp3m.local_mesh, weights, - [&E, &dp3m, d_rs](int ind, double w) { - E[d_rs] += w * double(dp3m.mesh.rs_scalar[ind]); - }); - - p.torque() -= vector_product(p.calc_dip(), prefac * E); + kernel(p.calc_dip() * prefac, p.torque(), p_index); ++p_index; } } +#endif // SHARED_MEMORY_PARALLELISM } }; template struct AssignForces { void operator()(auto &dp3m, double prefac, int d_rs, - ParticleRange const &particles) const { +#ifdef SHARED_MEMORY_PARALLELISM + CellStructure &cell_structure +#else + ParticleRange const &particles +#endif + ) const { + + assert(cao == dp3m.inter_weights.cao()); + auto const kernel = [d_rs, &dp3m](auto const &pref, auto &p_force, + std::size_t p_index) { + auto const weights = dp3m.inter_weights.template load(p_index); + + Utils::Vector3d E{}; + p3m_interpolate(dp3m.local_mesh, weights, [&E, &dp3m](int ind, double w) { + E[0u] += w * double(dp3m.mesh.rs_fields[0u][ind]); + E[1u] += w * double(dp3m.mesh.rs_fields[1u][ind]); + E[2u] += w * double(dp3m.mesh.rs_fields[2u][ind]); + }); + +#ifdef SHARED_MEMORY_PARALLELISM + auto const thread_id = omp_get_thread_num(); + p_force(p_index, thread_id, d_rs) += pref * E; +#else + p_force[d_rs] += pref * E; +#endif + }; + +#ifdef SHARED_MEMORY_PARALLELISM + auto const n_part = dp3m.inter_weights.size(); + auto const &unique_particles = cell_structure.get_unique_particles(); + auto &local_force = cell_structure.get_local_force(); + kokkos_parallel_range_for( + "AssignForces", std::size_t{0u}, n_part, [&](std::size_t p_index) { + auto const &p = *unique_particles.at(p_index); + if (p.dipm() != 0.) { + kernel(p.calc_dip() * prefac, local_force, p_index); + } + }); +#else // SHARED_MEMORY_PARALLELISM /* magnetic particle index */ auto p_index = std::size_t{0ul}; for (auto &p : particles) { if (p.dipm() != 0.) { - auto const weights = dp3m.inter_weights.template load(p_index); - - Utils::Vector3d E{}; - p3m_interpolate(dp3m.local_mesh, weights, - [&E, &dp3m](int ind, double w) { - E[0u] += w * double(dp3m.mesh.rs_fields[0u][ind]); - E[1u] += w * double(dp3m.mesh.rs_fields[1u][ind]); - E[2u] += w * double(dp3m.mesh.rs_fields[2u][ind]); - }); - - p.force()[d_rs] += p.calc_dip() * prefac * E; + kernel(p.calc_dip() * prefac, p.force(), p_index); ++p_index; } } +#endif // SHARED_MEMORY_PARALLELISM } }; } // namespace @@ -255,10 +360,7 @@ double DipolarP3MImpl::long_range_kernel( auto const &box_geo = *system.box_geo; auto const dipole_prefac = prefactor / Utils::int_pow<3>(dp3m.params.mesh[0]); #ifdef NPT - auto const npt_flag = - force_flag and - ((system.propagation->integ_switch == INTEG_METHOD_NPT_ISO_AND) or - (system.propagation->integ_switch == INTEG_METHOD_NPT_ISO_MTK)); + auto const npt_flag = force_flag and system.has_npt_enabled(); #else auto constexpr npt_flag = false; #endif @@ -376,8 +478,14 @@ double DipolarP3MImpl::long_range_kernel( dp3m.fft_buffers->perform_scalar_halo_spread(); /* Assign force component from mesh to particle */ auto const d_rs = (d + dp3m.mesh.ks_pnum) % 3; +#ifdef SHARED_MEMORY_PARALLELISM + auto &particle_data = *system.cell_structure; +#else + auto &particle_data = particles; +#endif Utils::integral_parameter( - dp3m.params.cao, dp3m, dipole_prefac * wavenumber, d_rs, particles); + dp3m.params.cao, dp3m, dipole_prefac * wavenumber, d_rs, + particle_data); } /*************************** @@ -432,9 +540,14 @@ double DipolarP3MImpl::long_range_kernel( dp3m.fft_buffers->perform_vector_halo_spread(); /* Assign force component from mesh to particle */ auto const d_rs = (d + dp3m.mesh.ks_pnum) % 3; +#ifdef SHARED_MEMORY_PARALLELISM + auto &particle_data = *system.cell_structure; +#else + auto &particle_data = particles; +#endif Utils::integral_parameter( dp3m.params.cao, dp3m, dipole_prefac * Utils::sqr(wavenumber), d_rs, - particles); + particle_data); } } /* if (dp3m.sum_mu2 > 0) */ } /* if (force_flag) */ diff --git a/src/core/magnetostatics/dp3m.hpp b/src/core/magnetostatics/dp3m.hpp index b8935799277..919d8482215 100644 --- a/src/core/magnetostatics/dp3m.hpp +++ b/src/core/magnetostatics/dp3m.hpp @@ -134,15 +134,13 @@ struct DipolarP3M : public Dipoles::Actor { /** Calculate real-space contribution of p3m dipolar pair forces and torques. * If NPT is compiled in, update the NpT virial. */ - inline ParticleForce pair_force(Particle const &p1, Particle const &p2, - Utils::Vector3d const &d, double dist2, - double dist) const { - if ((p1.dipm() == 0.) || (p2.dipm() == 0.) || dist >= dp3m_params.r_cut || - dist <= 0.) + inline ParticleForce pair_force(double d1d2, Utils::Vector3d const &dip1, + Utils::Vector3d const &dip2, + Utils::Vector3d const &d, double dist, + double dist2) const { + if (d1d2 == 0. or dist >= dp3m_params.r_cut or dist <= 0.) return {}; - auto const dip1 = p1.calc_dip(); - auto const dip2 = p2.calc_dip(); auto const alpsq = dp3m_params.alpha * dp3m_params.alpha; auto const adist = dp3m_params.alpha * dist; #if USE_ERFC_APPROXIMATION @@ -181,9 +179,9 @@ struct DipolarP3M : public Dipoles::Actor { auto const torque = prefactor * (-mixmj * B_r + mixr * (mjr * C_r)); #ifdef NPT #if USE_ERFC_APPROXIMATION - auto const fac = prefactor * p1.dipm() * p2.dipm() * exp_adist2; + auto const fac = prefactor * d1d2 * exp_adist2; #else - auto const fac = prefactor * p1.dipm() * p2.dipm(); + auto const fac = prefactor * d1d2; #endif auto const energy = fac * (mimj * B_r - mir * mjr * C_r); npt_add_virial_contribution(energy); @@ -193,9 +191,9 @@ struct DipolarP3M : public Dipoles::Actor { /** Calculate real-space contribution of dipolar pair energy. */ inline double pair_energy(Particle const &p1, Particle const &p2, - Utils::Vector3d const &d, double dist2, - double dist) const { - if ((p1.dipm() == 0.) || (p2.dipm() == 0.) || dist >= dp3m_params.r_cut || + Utils::Vector3d const &d, double dist, + double dist2) const { + if (p1.dipm() == 0. or p2.dipm() == 0. or dist >= dp3m_params.r_cut or dist <= 0.) return {}; diff --git a/src/core/magnetostatics/solver.hpp b/src/core/magnetostatics/solver.hpp index 58856326d0f..f5111408921 100644 --- a/src/core/magnetostatics/solver.hpp +++ b/src/core/magnetostatics/solver.hpp @@ -19,7 +19,7 @@ #pragma once -#include "config/config.hpp" +#include #include "actor/optional.hpp" #include "actor/traits.hpp" @@ -64,11 +64,12 @@ struct Solver { Solver(); #else // DIPOLES Solver() = default; + constexpr double cutoff() const { return INACTIVE_CUTOFF; } #endif // DIPOLES - using ShortRangeForceKernel = - std::function; + using ShortRangeForceKernel = std::function; using ShortRangeEnergyKernel = std::function; diff --git a/src/core/nonbonded_interactions/VerletCriterion.hpp b/src/core/nonbonded_interactions/VerletCriterion.hpp index a324d673dbc..484f902a9fd 100644 --- a/src/core/nonbonded_interactions/VerletCriterion.hpp +++ b/src/core/nonbonded_interactions/VerletCriterion.hpp @@ -21,8 +21,9 @@ #pragma once +#include + #include "Particle.hpp" -#include "config/config.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "system/System.hpp" diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp index 05aa9784e9b..0f14a5b3ec0 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp @@ -25,8 +25,9 @@ * Various procedures concerning interactions between particles. */ +#include + #include "TabulatedPotential.hpp" -#include "config/config.hpp" #include "system/Leaf.hpp" #include @@ -38,11 +39,6 @@ #include #include -/** Cutoff for deactivated interactions. Must be negative, so that even - * particles on top of each other don't interact by chance. - */ -constexpr double INACTIVE_CUTOFF = -1.; - /** Lennard-Jones with shift */ struct LJ_Parameters { double eps = 0.0; diff --git a/src/core/npt.cpp b/src/core/npt.cpp index dcdaf080a13..c0bdce72fdd 100644 --- a/src/core/npt.cpp +++ b/src/core/npt.cpp @@ -135,25 +135,8 @@ void System::System::npt_ensemble_init(bool recalc_forces) { } void System::System::npt_add_virial_contribution(double energy) { - if ((propagation->integ_switch == INTEG_METHOD_NPT_ISO_AND) or - (propagation->integ_switch == INTEG_METHOD_NPT_ISO_MTK)) { + if (has_npt_enabled()) { npt_inst_pressure->p_vir[0] += energy; } } - -void System::System::npt_add_virial_contribution(Utils::Vector3d const &force, - Utils::Vector3d const &d) { - if ((propagation->integ_switch == INTEG_METHOD_NPT_ISO_AND) or - (propagation->integ_switch == INTEG_METHOD_NPT_ISO_MTK)) { - 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/p3m/CMakeLists.txt b/src/core/p3m/CMakeLists.txt index f8f04882e7e..d7e817abfc8 100644 --- a/src/core/p3m/CMakeLists.txt +++ b/src/core/p3m/CMakeLists.txt @@ -43,6 +43,7 @@ if(ESPRESSO_BUILD_WITH_FFTW) espresso::instrumentation espresso::utils espresso::config + Random123 Boost::mpi MPI::MPI_CXX Heffte::Heffte diff --git a/src/core/p3m/interpolation.hpp b/src/core/p3m/interpolation.hpp index 42edd8eca65..b07ec3cf91d 100644 --- a/src/core/p3m/interpolation.hpp +++ b/src/core/p3m/interpolation.hpp @@ -72,6 +72,18 @@ class p3m_interpolation_cache { */ auto cao() const { return m_cao; } + /** + * @brief Fill cache with zero-initialized data. + * Meant for the parallel weight interpolation calculation. + * @param size Number of elements. + */ + void zfill(std::size_t size) { + assert(ca_frac.empty()); + assert(ca_fmp.empty()); + ca_fmp.resize(size); + ca_frac.resize(size * static_cast(m_cao * 3)); + } + /** * @brief Push back weights for one point. * @@ -89,6 +101,28 @@ class p3m_interpolation_cache { std::ranges::copy(weights.w_z, it); } + /** + * @brief Insert weights for one point. + * + * @tparam cao Interpolation order has to match the order + * set at last call to @ref p3m_interpolation_cache::reset. + * @param weights Interpolation weights to store. + */ + template + void store_at(std::size_t p_index, const InterpolationWeights &weights) { + assert(cao == m_cao); + assert(p_index < size()); + + ca_fmp[p_index] = weights.ind; + auto it = ca_frac.begin(); + std::advance(it, p_index * static_cast(m_cao * 3)); + std::ranges::copy(weights.w_x, it); + std::advance(it, m_cao); + std::ranges::copy(weights.w_y, it); + std::advance(it, m_cao); + std::ranges::copy(weights.w_z, it); + } + /** * @brief Load entry from the cache. * @@ -98,7 +132,7 @@ class p3m_interpolation_cache { * @tparam cao Interpolation order has to match the order * set at last call to @ref p3m_interpolation_cache::reset. * @param p_index Index of the entry to load. - * @return i-it interpolation weights. + * @return Interpolation weights. */ template InterpolationWeights load(std::size_t p_index) const { assert(cao == m_cao); diff --git a/src/core/particle_reduction.hpp b/src/core/particle_reduction.hpp index 9756257c60a..8f9abcad9d1 100644 --- a/src/core/particle_reduction.hpp +++ b/src/core/particle_reduction.hpp @@ -37,7 +37,7 @@ namespace Reduction { /** @brief Kernel that adds the result from a single particle to a reduction */ template using AddPartialResultKernel = - std::function; + std::function; /** @brief Join two partial reduciton results */ template @@ -111,7 +111,7 @@ ResultType reduce_over_local_particles( auto reducer = Reduction::make_kokkos_reducer( [&cells, add_partial](std::size_t const c_index, ResultType &res) { for (auto const &p : cells[c_index]->particles()) { - add_partial(p, res); + add_partial(res, p); } }, reduce_op); @@ -123,7 +123,7 @@ ResultType reduce_over_local_particles( auto const &particles = cells.front()->particles(); auto reducer = Reduction::make_kokkos_reducer( [&particles, add_partial](std::size_t const p_index, ResultType &res) { - add_partial(std::as_const(*(particles.begin() + p_index)), res); + add_partial(res, std::as_const(*(particles.begin() + p_index))); }, reduce_op); Kokkos::parallel_reduce( // loop over particles @@ -131,7 +131,7 @@ ResultType reduce_over_local_particles( return result; #else // SHARED_MEMORY_PARALLELISM for (auto const &p : cs.local_particles()) { - add_partial(p, result); + add_partial(result, p); } return result; #endif diff --git a/src/core/short_range_cabana.hpp b/src/core/short_range_cabana.hpp index 1f1b7593fe0..6dc530670e0 100644 --- a/src/core/short_range_cabana.hpp +++ b/src/core/short_range_cabana.hpp @@ -19,7 +19,7 @@ #pragma once -#include "config/config.hpp" +#include #ifdef SHARED_MEMORY_PARALLELISM @@ -36,8 +36,22 @@ #include #include +template > ESPRESSO_ATTR_ALWAYS_INLINE inline void -commit_particle(Particle const &p, int const index, +kokkos_parallel_range_for(auto const &name, auto start, auto end, + auto const &kernel) { + if (Kokkos::num_threads() > 1) { + KokkosRangePolicy policy(start, end); + Kokkos::parallel_for(name, policy, kernel); + } else { + for (auto p_index = start; p_index < end; ++p_index) { + kernel(p_index); + } + } +} + +ESPRESSO_ATTR_ALWAYS_INLINE inline void +commit_particle(Particle const &p, auto const index, CellStructure::AoSoA_pack &aosoa) { aosoa.id(index) = p.id(); #ifdef ELECTROSTATICS @@ -122,10 +136,10 @@ 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); + using policy_type = Kokkos::RangePolicy; + auto const rebuild = cell_structure.prepare_verlet_list_cabana(pair_cutoff); auto const &unique_particles = cell_structure.get_unique_particles(); + auto const n_part = unique_particles.size(); auto const max_id = cell_structure.get_cached_max_local_particle_id(); auto &aosoa = cell_structure.get_aosoa(); @@ -136,9 +150,8 @@ update_cabana_state(CellStructure &cell_structure, auto const &verlet_criterion, 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()), + kokkos_parallel_range_for( + "AoSoA write", std::size_t{0}, n_part, [&unique_particles, &aosoa, &id_to_index](int const index) { auto const &p = *unique_particles.at(index); commit_particle(p, index, aosoa); @@ -159,6 +172,23 @@ update_cabana_state(CellStructure &cell_structure, auto const &verlet_criterion, } } +#ifdef ELECTROSTATICS +ESPRESSO_ATTR_ALWAYS_INLINE inline void +update_aosoa_charges(CellStructure &cell_structure) { + using execution_space = Kokkos::DefaultExecutionSpace; + using policy_type = Kokkos::RangePolicy; + auto const &unique_particles = cell_structure.get_unique_particles(); + auto const n_part = unique_particles.size(); + auto &aosoa = cell_structure.get_aosoa(); + + kokkos_parallel_range_for( + "AoSoA update charges", std::size_t{0}, n_part, + [&unique_particles, &aosoa](std::size_t const index) { + aosoa.charge(index) = unique_particles.at(index)->q(); + }); +} +#endif + void cabana_short_range(auto const &bond_kernel, auto const &forces_kernel, CellStructure &cell_structure, double pair_cutoff, double bond_cutoff) { @@ -173,7 +203,7 @@ void cabana_short_range(auto const &bond_kernel, auto const &forces_kernel, if (pair_cutoff > 0.) { auto const &verlet_list = cell_structure.get_verlet_list_cabana(); Kokkos::RangePolicy policy( - 0, cell_structure.get_unique_particles().size()); + std::size_t{0}, cell_structure.get_unique_particles().size()); Cabana::neighbor_parallel_for(policy, forces_kernel, verlet_list, Cabana::FirstNeighborsTag(), Cabana::SerialOpTag()); diff --git a/src/core/system/System.cpp b/src/core/system/System.cpp index c96819568a4..9f05f5e2b4d 100644 --- a/src/core/system/System.cpp +++ b/src/core/system/System.cpp @@ -35,8 +35,10 @@ #include "communication.hpp" #include "electrostatics/icc.hpp" #include "errorhandling.hpp" +#include "nonbonded_interactions/VerletCriterion.hpp" #include "npt.hpp" #include "particle_node.hpp" +#include "short_range_cabana.hpp" #include "thermostat.hpp" #include "virtual_sites/relative.hpp" @@ -343,6 +345,9 @@ void System::on_particle_change() { /* the particle information is no longer valid */ invalidate_fetch_cache(); +#ifdef SHARED_MEMORY_PARALLELISM + cell_structure->clear_local_properties(); +#endif } void System::on_particle_charge_change() { @@ -360,7 +365,12 @@ void System::update_dependent_particles() { #endif #ifdef ELECTROSTATICS - update_icc_particles(); + if (has_icc_enabled()) { +#ifdef SHARED_MEMORY_PARALLELISM + rebuild_aosoa(); +#endif + update_icc_particles(); + } #endif // Here we initialize volume conservation @@ -384,8 +394,31 @@ void System::on_observable_calc() { #endif clear_particle_node(); +#ifdef SHARED_MEMORY_PARALLELISM + rebuild_aosoa(); +#endif } +#ifdef SHARED_MEMORY_PARALLELISM +void System::rebuild_aosoa() { +#ifdef COLLISION_DETECTION + auto const collision_detection_cutoff = collision_detection->cutoff(); +#else + auto const collision_detection_cutoff = INACTIVE_CUTOFF; +#endif + + VerletCriterion<> const verlet_criterion{*this, + cell_structure->get_verlet_skin(), + get_interaction_range(), + coulomb.cutoff(), + dipoles.cutoff(), + collision_detection_cutoff}; + + update_cabana_state(*cell_structure, verlet_criterion, + get_interaction_range()); +} +#endif // SHARED_MEMORY_PARALLELISM + void System::on_lees_edwards_change() { lb.on_lees_edwards_change(); } void System::update_local_geo() { @@ -397,12 +430,8 @@ void System::update_local_geo() { double System::maximal_cutoff() const { auto max_cut = INACTIVE_CUTOFF; max_cut = std::max(max_cut, get_min_global_cut()); -#ifdef ELECTROSTATICS max_cut = std::max(max_cut, coulomb.cutoff()); -#endif -#ifdef DIPOLES max_cut = std::max(max_cut, dipoles.cutoff()); -#endif if (::communicator.size > 1) { // If there is just one node, the bonded cutoff can be omitted // because bond partners are always on the local node. @@ -465,6 +494,9 @@ void System::on_integration_start() { } invalidate_fetch_cache(); +#ifdef SHARED_MEMORY_PARALLELISM + cell_structure->clear_local_properties(); +#endif #ifdef ADDITIONAL_CHECKS if (!Utils::Mpi::all_compare(::comm_cart, cell_structure->use_verlet_list)) { @@ -520,4 +552,20 @@ unsigned System::get_global_ghost_flags() const { return data_parts; } +#ifdef NPT +bool System::System::has_npt_enabled() const { + return (propagation->integ_switch == INTEG_METHOD_NPT_ISO_AND) or + (propagation->integ_switch == INTEG_METHOD_NPT_ISO_MTK); +} +#endif + +Utils::Vector3d *System::System::get_npt_virial() const { +#ifdef NPT + if (has_npt_enabled()) { + return &npt_inst_pressure->p_vir; + } +#endif + return nullptr; +} + } // namespace System diff --git a/src/core/system/System.hpp b/src/core/system/System.hpp index be53a43e7b0..4797fde87eb 100644 --- a/src/core/system/System.hpp +++ b/src/core/system/System.hpp @@ -19,7 +19,7 @@ #pragma once -#include "config/config.hpp" +#include #include "GpuParticleData.hpp" #include "ResourceCleanup.hpp" @@ -131,6 +131,9 @@ class System : public std::enable_shared_from_this { /** @brief Rebuild cell lists. Use e.g. after a skin change. */ void rebuild_cell_structure(); +#ifdef SHARED_MEMORY_PARALLELISM + void rebuild_aosoa(); +#endif /** @brief Calculate the maximal cutoff of all interactions. */ double maximal_cutoff() const; @@ -163,10 +166,9 @@ class System : public std::enable_shared_from_this { /** @brief Reinitialize the NpT state. */ void npt_ensemble_init(bool recalc_forces); 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); + bool has_npt_enabled() const; #endif // NPT + Utils::Vector3d *get_npt_virial() const; /** @brief Calculate all forces. */ void calculate_forces(); @@ -352,6 +354,7 @@ class System : public std::enable_shared_from_this { void update_local_geo(); #ifdef ELECTROSTATICS void update_icc_particles(); + bool has_icc_enabled() const; #endif // ELECTROSTATICS private: diff --git a/src/core/unit_tests/CMakeLists.txt b/src/core/unit_tests/CMakeLists.txt index d24b10cd6aa..db1d6b621dc 100644 --- a/src/core/unit_tests/CMakeLists.txt +++ b/src/core/unit_tests/CMakeLists.txt @@ -24,14 +24,22 @@ include(espresso_unit_test) espresso_unit_test(SRC RuntimeError_test.cpp DEPENDS Boost::serialization) espresso_unit_test(SRC RuntimeErrorCollector_test.cpp DEPENDS espresso::core Boost::mpi MPI::MPI_CXX NUM_PROC 2) -espresso_unit_test( - NAME EspressoSystemStandAlone_parallel_test SRC - EspressoSystemStandAlone_test.cpp DEPENDS espresso::core Boost::mpi - MPI::MPI_CXX NUM_PROC 4) -espresso_unit_test( - NAME EspressoSystemStandAlone_serial_test SRC - EspressoSystemStandAlone_test.cpp DEPENDS espresso::core Boost::mpi - MPI::MPI_CXX NUM_PROC 1) +foreach(TEST_NUM_PROC 1 3 4) + foreach(TEST_NUM_THREADS 1 2) + if(${ESPRESSO_BUILD_WITH_SHARED_MEMORY_PARALLELISM} OR ${TEST_NUM_THREADS} + EQUAL 1) + if(${ESPRESSO_TEST_NP} GREATER_EQUAL ${TEST_NUM_PROC}) + set(TEST_NAME + EspressoSystemStandAlone_test_${TEST_NUM_PROC}_mpi_ranks_${TEST_NUM_THREADS}_omp_threads + ) + espresso_unit_test( + NAME ${TEST_NAME} SRC EspressoSystemStandAlone_test.cpp DEPENDS + espresso::core Boost::mpi MPI::MPI_CXX NUM_PROC ${TEST_NUM_PROC} + NUM_THREADS ${TEST_NUM_THREADS}) + endif() + endif() + endforeach() +endforeach() espresso_unit_test(SRC EspressoSystem_test.cpp DEPENDS espresso::core Boost::mpi NUM_PROC 2) espresso_unit_test(SRC ResourceCleanup_test.cpp DEPENDS espresso::core diff --git a/src/core/unit_tests/particle_reduction_test.cpp b/src/core/unit_tests/particle_reduction_test.cpp index 763eee7290f..7fa3511ebdf 100644 --- a/src/core/unit_tests/particle_reduction_test.cpp +++ b/src/core/unit_tests/particle_reduction_test.cpp @@ -99,8 +99,8 @@ BOOST_AUTO_TEST_CASE(test_reduce_over_local_particles) { auto const *const p = cell_structure.get_local_particle(0); assert(p); - auto const kernel = [](Particle const &p, Utils::Vector3d &res) { - res += p.mass() * p.v(); + auto const kernel = [](Utils::Vector3d &acc, Particle const &p) { + acc += p.mass() * p.v(); }; auto const ref = p->mass() * p->v(); auto const res = reduce_over_local_particles( diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 517a10dc1c5..46b19ae2c91 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -232,6 +232,7 @@ python_test(FILE engine_langevin.py MAX_NUM_PROC 4) python_test(FILE engine_lb.py MAX_NUM_PROC 2 GPU_SLOTS 1) python_test(FILE engine_lb.py MAX_NUM_PROC 1 GPU_SLOTS 1 SUFFIX 1_core) python_test(FILE icc.py MAX_NUM_PROC 4) +python_test(FILE icc_electrodes.py MAX_NUM_PROC 1) python_test(FILE icc_interface.py MAX_NUM_PROC 1 GPU_SLOTS 1) python_test(FILE mass-and-rinertia_per_particle.py MAX_NUM_PROC 2 LABELS long) python_test(FILE interactions_bond_angle.py MAX_NUM_PROC 4) diff --git a/testsuite/python/caliper.py b/testsuite/python/caliper.py index a7dc44784cb..30771e10a02 100644 --- a/testsuite/python/caliper.py +++ b/testsuite/python/caliper.py @@ -25,23 +25,25 @@ import sys import os +HAS_CABANA = espressomd.has_features(["SHARED_MEMORY_PARALLELISM"]) + EXPECTED_LABELS = f""" integrate Initial Force Calculation calculate_forces copy_particles_to_GPU - init_forces_and_thermost + {'convert particles AoS to SoA' if HAS_CABANA else ''} + init_forces_and_thermostat calc_long_range_forces - {'parallel short range' if espressomd.has_features( - ["SHARED_MEMORY_PARALLELISM"]) else 'short_range_loop'} + {'parallel short range' if HAS_CABANA else 'short_range_loop'} copy_forces_from_GPU Integration loop calculate_forces copy_particles_to_GPU - init_forces_and_thermost + {'convert particles AoS to SoA' if HAS_CABANA else ''} + init_forces_and_thermostat calc_long_range_forces - {'parallel short range' if espressomd.has_features( - ["SHARED_MEMORY_PARALLELISM"]) else 'short_range_loop'} + {'parallel short range' if HAS_CABANA else 'short_range_loop'} copy_forces_from_GPU calc_energies short_range_loop @@ -70,9 +72,9 @@ def test_runtime_report(self): header = "Path\tMin time/rank\tMax time/rank\tAvg time/rank\tTime %" self.assertEqual(lines[0].split(), header.split(), msg=f"Caliper summary should start with '{header}'") - labels = [line[:30].strip() for line in lines[1:]] - labels_ref = [x.strip() for x in EXPECTED_LABELS.strip().split("\n") - if "GPU" not in x.upper() or has_cuda] + labels = [line[:36].rstrip() for line in lines[1:]] + labels_ref = [x.rstrip() for x in EXPECTED_LABELS.strip().split("\n") + if x.rstrip() and ("GPU" not in x.upper() or has_cuda)] self.assertEqual(labels[:len(labels_ref)], labels_ref, msg=f"Caliper returned this summary:\n{stderr}") diff --git a/testsuite/python/icc_electrodes.py b/testsuite/python/icc_electrodes.py new file mode 100644 index 00000000000..f273abd2090 --- /dev/null +++ b/testsuite/python/icc_electrodes.py @@ -0,0 +1,110 @@ +# +# Copyright (C) 2023-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 . +# + +import unittest as ut +import unittest_decorators as utx +import espressomd +import espressomd.electrostatics +import espressomd.electrostatic_extensions +import numpy as np + + +@utx.skipIfMissingFeatures(["ELECTROSTATICS", "EXTERNAL_FORCES"]) +class Test(ut.TestCase): + + def test_electrodes(self): + """ + This is a simplified version of the electrodes tutorial part 1. + Two particles between the electrode plates are moving away from each + other. Their force is compared against the analytical solution. + """ + box_l_x = 80. + box_l_y = 80. + box_l_z = 5. + elc_gap = 6. * box_l_z + + system = espressomd.System(box_l=[box_l_x, box_l_y, box_l_z + elc_gap]) + system.time_step = 0.01 + system.cell_system.skin = 0.4 + + pos1 = [box_l_x * 0.25, box_l_y / 2., box_l_z / 2.] + pos2 = [box_l_x * 0.75, box_l_y / 2., box_l_z / 2.] + p1 = system.part.add(pos=pos1, q=+1., fix=3 * [True]) + p2 = system.part.add(pos=pos2, q=-1., fix=3 * [True]) + + xs = np.linspace(0., box_l_x, num=int(box_l_x), endpoint=False) + ys = np.linspace(0., box_l_y, num=int(box_l_y), endpoint=False) + surf_q = 1. / (len(xs) * len(ys)) + icc_particles_bot = [] + icc_particles_top = [] + lz = box_l_z + for x in xs: + for y in ys: + icc_particles_bot.append( + system.part.add(pos=[x, y, 0.], q=-surf_q, fix=3 * [True])) + icc_particles_top.append( + system.part.add(pos=[x, y, lz], q=+surf_q, fix=3 * [True])) + + n_icc_part = len(icc_particles_top) + len(icc_particles_bot) + + p3m = espressomd.electrostatics.P3M( + prefactor=2., accuracy=1e-7, mesh=[100, 100, 96], cao=7, check_neutrality=False) + elc = espressomd.electrostatics.ELC( + actor=p3m, gap_size=elc_gap, maxPWerror=1e-7, check_neutrality=False) + system.electrostatics.solver = elc + icc = espressomd.electrostatic_extensions.ICC( + first_id=2, n_icc=n_icc_part, convergence=1e-3, eps_out=1., + ext_field=[0, 0, 0], max_iterations=100, relaxation=0.95, + normals=np.tile([[0, 0, 1], [0, 0, -1]], [n_icc_part // 2, 1]), + areas=np.ones(n_icc_part), sigmas=np.zeros(n_icc_part), + epsilons=1e5 * np.ones(n_icc_part) + ) + system.electrostatics.extension = icc + + n_points = 6 + r = np.logspace(0., box_l_z / 5., n_points) + elc_forces_axial = np.empty((n_points, 2)) + p1.pos = [0., box_l_y / 2., box_l_z / 2.] + for i in range(n_points): + p2.pos = [r[i], box_l_y / 2., box_l_z / 2.] + system.integrator.run(0, recalc_forces=True) + elc_forces_axial[i, 0] = p1.f[0] + elc_forces_axial[i, 1] = p2.f[0] + # reset ICC charges + for part in icc_particles_bot: + part.q = +surf_q + for part in icc_particles_bot: + part.q = -surf_q + + def calc_sum(x): + limit = 1000 + accumulator = 0. + for n in range(-limit + 1, limit + 1): + accumulator += (-1)**n * x / (x**2 + n**2)**(3. / 2.) + return accumulator + + ref_force = calc_sum(r / box_l_z) * p3m.prefactor / box_l_z**2 + np.testing.assert_allclose(np.log(elc_forces_axial[:, 0]), + np.log(ref_force), rtol=0.05) + np.testing.assert_allclose(np.log(-elc_forces_axial[:, 1]), + np.log(ref_force), rtol=0.05) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/python/icc_interface.py b/testsuite/python/icc_interface.py index c452d5b0ba7..6d13ca05f2e 100644 --- a/testsuite/python/icc_interface.py +++ b/testsuite/python/icc_interface.py @@ -237,7 +237,7 @@ def run_exceptions_npt(self, barostat): self.system.thermostat.set_npt(kT=1., gamma0=2., gammav=0.004, seed=42) self.system.integrator.set_isotropic_npt( ext_pressure=2., piston=0.001, barostat=barostat) - with self.assertRaisesRegex(RuntimeError, "ICC does not work in the NPT ensemble"): + with self.assertRaisesRegex(RuntimeError, "ICC does not work in the NpT ensemble"): self.system.electrostatics.extension = icc self.assertIsNone(self.system.electrostatics.extension) self.system.integrator.run(0) diff --git a/testsuite/python/long_range_actors.py b/testsuite/python/long_range_actors.py index 6ced5a0af0e..74616e5634f 100644 --- a/testsuite/python/long_range_actors.py +++ b/testsuite/python/long_range_actors.py @@ -105,7 +105,7 @@ def test_electrostatics_registration(self): self.system.electrostatics.extension = icc def run_with_npt(barostat): - with self.assertRaisesRegex(Exception, "ERROR: ICC does not work in the NPT ensemble"): + with self.assertRaisesRegex(Exception, "ERROR: ICC does not work in the NpT ensemble"): self.system.thermostat.set_npt( kT=1., gamma0=2., gammav=0., seed=42) self.system.integrator.set_isotropic_npt( diff --git a/testsuite/test_importlib_wrapper.py b/testsuite/test_importlib_wrapper.py index 076df361280..43b68efaf43 100644 --- a/testsuite/test_importlib_wrapper.py +++ b/testsuite/test_importlib_wrapper.py @@ -195,6 +195,14 @@ def test_mock_es_visualization(self): """ self.assertEqual(iw.mock_es_visualization(statement), expected[1:]) + statement = "import espressomd.zn" + expected = f"""\ +import unittest.mock +import espressomd +espressomd.zn = unittest.mock.MagicMock() +""" + self.assertEqual(iw.mock_es_visualization(statement, True), expected) + # test exceptions self.assertRaises(ValueError, iw.mock_es_visualization, "from espressomd.visualization import *")