script_interface/system/System.cpp

Go to the documentation of this file.

/*
 * Copyright (C) 2013-2026 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 <http://www.gnu.org/licenses/>.
 */
 
#include "System.hpp"
 
#include <config/config.hpp>
 
#include "core/BoxGeometry.hpp"
#include "core/Particle.hpp"
#include "core/ParticleRange.hpp"
#include "core/cell_system/CellStructure.hpp"
#include "core/cell_system/CellStructureType.hpp"
#include "core/cells.hpp"
#include "core/communication.hpp"
#include "core/exclusions.hpp"
#include "core/nonbonded_interactions/nonbonded_interaction_data.hpp"
#include "core/npt.hpp"
#include "core/object-in-fluid/oif_global_forces.hpp"
#include "core/particle_node.hpp"
#include "core/propagation.hpp"
#include "core/rotation.hpp"
#include "core/system/System.hpp"
#include "core/system/System.impl.hpp"
 
#include "script_interface/ObjectState.hpp"
#include "script_interface/accumulators/AutoUpdateAccumulators.hpp"
#include "script_interface/analysis/Analysis.hpp"
#include "script_interface/bond_breakage/BreakageSpecs.hpp"
#include "script_interface/cell_system/CellSystem.hpp"
#include "script_interface/collision_detection/CollisionDetection.hpp"
#include "script_interface/constraints/Constraints.hpp"
#include "script_interface/ek/Container.hpp"
#include "script_interface/electrostatics/Container.hpp"
#include "script_interface/galilei/ComFixed.hpp"
#include "script_interface/galilei/Galilei.hpp"
#include "script_interface/integrators/IntegratorHandle.hpp"
#include "script_interface/interactions/BondedInteractions.hpp"
#include "script_interface/interactions/NonBondedInteractions.hpp"
#include "script_interface/lb/Container.hpp"
#include "script_interface/lees_edwards/LeesEdwards.hpp"
#include "script_interface/magnetostatics/Container.hpp"
#include "script_interface/particle_data/ParticleHandle.hpp"
#include "script_interface/particle_data/ParticleList.hpp"
#include "script_interface/thermostat/thermostat.hpp"
 
#include <utils/Vector.hpp>
#include <utils/demangle.hpp>
#include <utils/math/vec_rotate.hpp>
#include <utils/mpi/gather_buffer.hpp>
#include <utils/serialization/pack.hpp>
 
#include <boost/mpi/collectives.hpp>
 
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <functional>
#include <initializer_list>
#include <memory>
#include <ranges>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <vector>
 
namespace ScriptInterface {
namespace System {
 
static bool system_created = false;
 
#ifdef ESPRESSO_EXCLUSIONS
static void set_exclusions(Particles::ParticleHandle &p,
                           Variant const &exclusions) {
  p.call_method("set_exclusions", {{"p_ids", exclusions}});
}
#endif // ESPRESSO_EXCLUSIONS
 
static void set_bonds(Particles::ParticleHandle &p, Variant const &bonds) {
  auto const bond_list_flat = get_value<std::vector<std::vector<int>>>(bonds);
  for (auto const &bond_flat : bond_list_flat) {
    auto const bond_id = bond_flat[0];
    auto const part_id =
        std::vector<int>{bond_flat.begin() + 1, bond_flat.end()};
    p.call_method("add_bond",
                  {{"bond_id", bond_id}, {"part_id", std::move(part_id)}});
  }
}
 
/** @brief Container for leaves of the system class. */
struct System::Leaves {
  Leaves() = default;
  std::shared_ptr<CellSystem::CellSystem> cell_system;
  std::shared_ptr<Integrators::IntegratorHandle> integrator;
  std::shared_ptr<Interactions::BondedInteractions> bonded_interactions;
#ifdef ESPRESSO_COLLISION_DETECTION
  std::shared_ptr<CollisionDetection::CollisionDetection> collision_detection;
#endif
  std::shared_ptr<Thermostat::Thermostat> thermostat;
  std::shared_ptr<Analysis::Analysis> analysis;
  std::shared_ptr<Galilei::ComFixed> comfixed;
  std::shared_ptr<Galilei::Galilei> galilei;
  std::shared_ptr<BondBreakage::BreakageSpecs> bond_breakage;
  std::shared_ptr<LeesEdwards::LeesEdwards> lees_edwards;
  std::shared_ptr<Accumulators::AutoUpdateAccumulators>
      auto_update_accumulators;
  std::shared_ptr<Constraints::Constraints> constraints;
  std::shared_ptr<Interactions::NonBondedInteractions> non_bonded_inter;
#ifdef ESPRESSO_ELECTROSTATICS
  std::shared_ptr<Coulomb::Container> electrostatics;
#endif
#ifdef ESPRESSO_DIPOLES
  std::shared_ptr<Dipoles::Container> magnetostatics;
#endif
  std::shared_ptr<LB::Container> lb;
  std::shared_ptr<EK::Container> ek;
  std::shared_ptr<Particles::ParticleList> part;
 
  ~Leaves() {
    // Clear containers whose elements call MPI callbacks upon destruction.
    // The containers lifetime is extended by the global context shared object
    // registry on the head node, which can cause MPI deadlocks if they still
    // contain elements.
    if (bonded_interactions) {
      bonded_interactions->clear();
    }
    if (constraints) {
      constraints->clear();
    }
  }
};
 
System::System() : m_instance{}, m_leaves{std::make_unique<Leaves>()} {
  auto const add_parameter =
      [this, ptr = m_leaves.get()](std::string key, auto Leaves::*member) {
        add_parameters({AutoParameter(
            key.c_str(),
            [this, ptr, member, key](Variant const &val) {
              auto &dst = ptr->*member;
              if (dst != nullptr) {
                throw WriteError(key);
              }
              dst = get_value<std::remove_reference_t<decltype(dst)>>(val);
              dst->bind_system(m_instance);
            },
            [ptr, member]() { return ptr->*member; })});
      };
 
  add_parameters({
      {"box_l",
       [this](Variant const &v) {
         context()->parallel_try_catch([&]() {
           auto const new_value = get_value<Utils::Vector3d>(v);
           if (not(new_value > Utils::Vector3d::broadcast(0.))) {
             throw std::domain_error("Attribute 'box_l' must be > 0");
           }
           m_instance->veto_boxl_change();
           m_instance->box_geo->set_length(new_value);
           m_instance->on_boxl_change();
         });
       },
       [this]() { return m_instance->box_geo->length(); }},
      {"periodicity",
       [this](Variant const &v) {
         auto const periodicity = get_value<Utils::Vector3b>(v);
         for (unsigned int i = 0u; i < 3u; ++i) {
           m_instance->box_geo->set_periodic(i, periodicity[i]);
         }
         context()->parallel_try_catch(
             [&]() { m_instance->on_periodicity_change(); });
       },
       [this]() {
         return Utils::Vector3b{m_instance->box_geo->periodic(0),
                                m_instance->box_geo->periodic(1),
                                m_instance->box_geo->periodic(2)};
       }},
      {"min_global_cut",
       [this](Variant const &v) {
         context()->parallel_try_catch([&]() {
           auto const new_value = get_value<double>(v);
           if (new_value < 0. and new_value != inactive_cutoff) {
             throw std::domain_error("Attribute 'min_global_cut' must be >= 0");
           }
           m_instance->set_min_global_cut(new_value);
         });
       },
       [this]() { return m_instance->get_min_global_cut(); }},
      {"max_oif_objects",
       [this](Variant const &v) {
         m_instance->oif_global->max_oif_objects = get_value<int>(v);
       },
       [this]() { return m_instance->oif_global->max_oif_objects; }},
 
  });
  // note: the order of leaves matters! e.g. bonds depend on thermostats,
  // and thus a thermostat object must be instantiated before the bonds
  add_parameter("cell_system", &Leaves::cell_system);
  add_parameter("integrator", &Leaves::integrator);
  add_parameter("thermostat", &Leaves::thermostat);
  add_parameter("analysis", &Leaves::analysis);
  add_parameter("comfixed", &Leaves::comfixed);
  add_parameter("galilei", &Leaves::galilei);
  add_parameter("bonded_inter", &Leaves::bonded_interactions);
#ifdef ESPRESSO_COLLISION_DETECTION
  add_parameter("collision_detection", &Leaves::collision_detection);
#endif
  add_parameter("bond_breakage", &Leaves::bond_breakage);
  add_parameter("lees_edwards", &Leaves::lees_edwards);
  add_parameter("auto_update_accumulators", &Leaves::auto_update_accumulators);
  add_parameter("constraints", &Leaves::constraints);
  add_parameter("non_bonded_inter", &Leaves::non_bonded_inter);
#ifdef ESPRESSO_ELECTROSTATICS
  add_parameter("electrostatics", &Leaves::electrostatics);
#endif
#ifdef ESPRESSO_DIPOLES
  add_parameter("magnetostatics", &Leaves::magnetostatics);
#endif
  add_parameter("lbcontainer", &Leaves::lb);
  add_parameters({
      {"ekcontainer",
       [this](Variant const &v) {
         if (is_none(v)) {
           m_leaves->ek->do_call_method("deactivate", {});
           if (not context()->is_head_node()) {
             return;
           }
           set_parameter("ekcontainer",
                         context()->make_shared("EK::Container", {}));
         } else {
           auto const detach_solver = [this]() {
             auto &solver = m_leaves->ek;
             if (solver) {
               solver->do_call_method("deactivate", {});
               solver->detach_system();
             }
             solver.reset();
           };
           auto const bind_solver = [this]() {
             auto &solver = m_leaves->ek;
             if (solver) {
               solver->bind_system(m_instance);
               solver->do_call_method("activate", {});
             }
           };
           auto &solver = m_leaves->ek;
           auto new_solver = get_value<std::shared_ptr<EK::Container>>(v);
           auto old_solver = solver;
           detach_solver();
           try {
             solver = new_solver;
             context()->parallel_try_catch([&]() { bind_solver(); });
           } catch (...) {
             detach_solver();
             solver = old_solver;
             bind_solver();
             throw;
           }
         }
       },
       [this]() { return m_leaves->ek; }},
  });
  add_parameter("part", &Leaves::part);
}
 
template <typename LeafType>
void System::do_set_default_parameter(std::string const &name) {
  assert(context()->is_head_node());
  auto const so_name = Utils::demangle<LeafType>().substr(17);
  set_parameter(name, Variant{context()->make_shared(so_name, {})});
}
 
void System::do_construct(VariantMap const &params) {
  /* When reloading the system state from a checkpoint file,
   * the order of global variables instantiation matters.
   * The @c box_l must be set before any other global variable.
   * All these globals re-initialize the cell system, and we
   * cannot use the default-constructed @c box_geo when e.g.
   * long-range interactions exist in the system, otherwise
   * runtime errors about the local geometry being smaller
   * than the interaction range would be raised.
   */
  context()->parallel_try_catch([&]() {
    if (not params.contains("box_l")) {
      throw std::domain_error("Required argument 'box_l' not provided.");
    }
    if (params.contains("_regular_constructor") and system_created) {
      throw std::runtime_error(
          "You can only have one instance of the system class at a time");
    }
  });
  m_instance = ::System::System::create();
  ::System::set_system(m_instance);
 
  // domain decomposition can only be set after box_l is set
  m_instance->set_cell_structure_topology(CellStructureType::NSQUARE);
  do_set_parameter("box_l", params.at("box_l"));
  m_instance->set_cell_structure_topology(CellStructureType::REGULAR);
 
  m_instance->lb.bind_system(m_instance);
  m_instance->ek.bind_system(m_instance);
 
  if (params.contains("_regular_constructor")) {
    std::set<std::string> const setable_properties = {
        "box_l",           "min_global_cut",
        "periodicity",     "time",
        "time_step",       "force_cap",
        "max_oif_objects", "_regular_constructor"};
    for (auto const &name : std::views::elements<0>(params)) {
      if (not setable_properties.contains(name)) {
        context()->parallel_try_catch([&name]() {
          throw std::domain_error(
              "Property '" + name +
              "' cannot be set via argument to System class");
        });
      }
    }
    for (std::string attr :
         {"min_global_cut", "periodicity", "max_oif_objects"}) {
      if (params.contains(attr)) {
        do_set_parameter(attr, params.at(attr));
      }
    }
    if (not context()->is_head_node()) {
      return;
    }
    auto integrator = std::dynamic_pointer_cast<Integrators::IntegratorHandle>(
        context()->make_shared("Integrators::IntegratorHandle", {}));
    set_parameter("integrator", integrator);
    for (std::string attr : {"time", "time_step", "force_cap"}) {
      if (params.contains(attr)) {
        integrator->set_parameter(attr, params.at(attr));
      }
    }
    // note: the order of leaves matters! e.g. bonds depend on thermostats,
    // and thus a thermostat object must be instantiated before the bonds
    do_set_default_parameter<CellSystem::CellSystem>("cell_system");
    do_set_default_parameter<Thermostat::Thermostat>("thermostat");
    do_set_default_parameter<Interactions::BondedInteractions>("bonded_inter");
#ifdef ESPRESSO_COLLISION_DETECTION
    do_set_default_parameter<CollisionDetection::CollisionDetection>(
        "collision_detection");
#endif
    do_set_default_parameter<Analysis::Analysis>("analysis");
    do_set_default_parameter<Galilei::ComFixed>("comfixed");
    do_set_default_parameter<Galilei::Galilei>("galilei");
    do_set_default_parameter<BondBreakage::BreakageSpecs>("bond_breakage");
    do_set_default_parameter<LeesEdwards::LeesEdwards>("lees_edwards");
    do_set_default_parameter<Accumulators::AutoUpdateAccumulators>(
        "auto_update_accumulators");
    do_set_default_parameter<Constraints::Constraints>("constraints");
    do_set_default_parameter<Interactions::NonBondedInteractions>(
        "non_bonded_inter");
#ifdef ESPRESSO_ELECTROSTATICS
    do_set_default_parameter<Coulomb::Container>("electrostatics");
#endif
#ifdef ESPRESSO_DIPOLES
    do_set_default_parameter<Dipoles::Container>("magnetostatics");
#endif
    do_set_default_parameter<LB::Container>("lbcontainer");
    do_set_default_parameter<EK::Container>("ekcontainer");
    do_set_default_parameter<Particles::ParticleList>("part");
  } else {
    for (auto const &key : get_parameter_insertion_order()) {
      if (key != "box_l" and params.contains(key)) {
        do_set_parameter(key, params.at(key));
      }
    }
  }
  if (not context()->is_head_node()) {
    return;
  }
  call_method("internal_attach_leaves", {});
}
 
static void rotate_system(CellStructure &cell_structure, double phi,
                          double theta, double alpha) {
  auto const particles = cell_structure.local_particles();
 
  // Calculate center of mass
  Utils::Vector3d local_com{};
  double local_mass = 0.0;
 
  for (auto const &p : particles) {
    if (not p.is_virtual()) {
      local_com += p.mass() * p.pos();
      local_mass += p.mass();
    }
  }
 
  auto const total_mass =
      boost::mpi::all_reduce(comm_cart, local_mass, std::plus<>());
  auto const com =
      boost::mpi::all_reduce(comm_cart, local_com, std::plus<>()) / total_mass;
 
  // Rotation axis in Cartesian coordinates
  Utils::Vector3d axis;
  axis[0] = std::sin(theta) * std::cos(phi);
  axis[1] = std::sin(theta) * std::sin(phi);
  axis[2] = std::cos(theta);
 
  // Rotate particle coordinates
  for (auto &p : particles) {
    // Move the center of mass of the system to the origin
    p.pos() = com + Utils::vec_rotate(axis, alpha, p.pos() - com);
#ifdef ESPRESSO_ROTATION
    local_rotate_particle(p, axis, alpha);
#endif
  }
 
  cell_structure.set_resort_particles(Cells::RESORT_GLOBAL);
  cell_structure.update_ghosts_and_resort_particle(Cells::DATA_PART_POSITION |
                                                   Cells::DATA_PART_PROPERTIES);
}
 
Variant System::do_call_method(std::string const &name,
                               VariantMap const &parameters) {
  if (name == "lock_system_creation") {
    system_created = true;
    return {};
  }
  if (name == "rescale_boxl") {
    auto const rescale_particles = [this](unsigned dir, double scale) {
      for (auto &p : m_instance->cell_structure->local_particles()) {
        p.pos()[dir] *= scale;
      }
    };
    auto const &box_geo = *m_instance->box_geo;
    auto const coord = get_value<int>(parameters, "coord");
    auto const length = get_value<double>(parameters, "length");
    assert(coord >= 0 and coord <= 3);
    context()->parallel_try_catch([&]() {
      if (length <= 0.) {
        throw std::domain_error("Parameter 'd_new' must be > 0");
      }
      m_instance->veto_boxl_change(true);
    });
    if (coord == 3) {
      auto const old_length_inv = box_geo.length_inv();
      auto const new_value = Utils::Vector3d::broadcast(length);
      auto const rescale = Utils::hadamard_product(new_value, old_length_inv);
      // when shrinking, rescale particles before the box resize
      for (auto axis = 0u; axis < 3u; ++axis) {
        if (rescale[axis] <= 1.) {
          rescale_particles(axis, rescale[axis]);
        }
      }
      m_instance->on_particle_change();
      m_instance->box_geo->set_length(new_value);
      m_instance->on_boxl_change();
      // when growing, rescale particles after the box resize
      for (auto axis = 0u; axis < 3u; ++axis) {
        if (rescale[axis] > 1.) {
          rescale_particles(axis, rescale[axis]);
        }
      }
      m_instance->on_particle_change();
    } else {
      auto const axis = static_cast<unsigned>(coord);
      auto const scale = length * box_geo.length_inv()[axis];
      auto new_value = box_geo.length();
      new_value[axis] = length;
      // when shrinking, rescale particles before the box resize
      if (scale <= 1.) {
        rescale_particles(axis, scale);
        m_instance->on_particle_change();
      }
      m_instance->box_geo->set_length(new_value);
      m_instance->on_boxl_change();
      // when growing, rescale particles after the box resize
      if (scale > 1.) {
        rescale_particles(axis, scale);
        m_instance->on_particle_change();
      }
    }
    return {};
  }
  if (name == "reaction_get_maximal_particle_id") {
    return ::get_maximal_particle_id();
  }
  if (name == "get_pids_of_type") {
    auto const type = get_value<int>(parameters, "ptype");
    std::vector<int> pids;
    for (auto const &p : get_system().cell_structure->local_particles()) {
      if (p.type() == type) {
        pids.emplace_back(p.id());
      }
    }
    Utils::Mpi::gather_buffer(pids, context()->get_comm());
    return pids;
  }
  if (name == "number_of_particles") {
    auto const type = get_value<int>(parameters, "type");
    int local_counter = 0;
    int global_counter = 0;
    for (auto const &p : get_system().cell_structure->local_particles()) {
      if (p.type() == type) {
        local_counter++;
      }
    }
    boost::mpi::reduce(::comm_cart, local_counter, global_counter,
                       std::plus<>(), 0);
    return global_counter;
  }
  if (name == "velocity_difference") {
    auto const pos1 = get_value<Utils::Vector3d>(parameters, "pos1");
    auto const pos2 = get_value<Utils::Vector3d>(parameters, "pos2");
    auto const v1 = get_value<Utils::Vector3d>(parameters, "v1");
    auto const v2 = get_value<Utils::Vector3d>(parameters, "v2");
    return m_instance->box_geo->velocity_difference(pos2, pos1, v2, v1);
  }
  if (name == "distance_vec") {
    auto const pos1 = get_value<Utils::Vector3d>(parameters, "pos1");
    auto const pos2 = get_value<Utils::Vector3d>(parameters, "pos2");
    return m_instance->box_geo->get_mi_vector(pos2, pos1);
  }
  if (name == "rotate_system") {
    rotate_system(*m_instance->cell_structure,
                  get_value<double>(parameters, "phi"),
                  get_value<double>(parameters, "theta"),
                  get_value<double>(parameters, "alpha"));
    m_instance->on_particle_change();
    m_instance->update_dependent_particles();
    return {};
  }
  if (name == "get_propagation_modes_enum") {
    return make_unordered_map_of_variants(propagation_flags_map());
  }
  if (name == "session_shutdown") {
    if (m_instance) {
      if (::System::is_same_system(m_instance.get())) {
        ::System::reset_system();
      }
      assert(m_instance.use_count() == 1l);
      m_leaves.reset();
      m_instance.reset();
    }
    return {};
  }
  if (name == "internal_attach_leaves") {
    m_leaves->part->attach(m_leaves->cell_system,
                           m_leaves->bonded_interactions);
#ifdef ESPRESSO_COLLISION_DETECTION
    m_leaves->collision_detection->attach(m_leaves->bonded_interactions);
#endif
    return {};
  }
  return {};
}
 
/**
 * @brief Serialize particles.
 * Particles need to be serialized here to reduce overhead,
 * and also to guarantee particles get instantiated after the cell structure
 * was instantiated (since they store a weak pointer to it).
 */
std::string System::get_internal_state() const {
  auto const p_ids = get_particle_ids();
  std::vector<std::string> object_states(p_ids.size());
 
  std::ranges::transform(p_ids, object_states.begin(), [this](auto const p_id) {
    auto p_obj = context()->make_shared(
        "Particles::ParticleHandle",
        {{"id", p_id}, {"__cell_structure", m_leaves->cell_system}});
    auto &p_handle = dynamic_cast<Particles::ParticleHandle &>(*p_obj);
    auto const packed_state = p_handle.serialize();
    // custom particle serialization
    auto state = Utils::unpack<ObjectState>(packed_state);
    state.name = "Particles::ParticleHandle";
    auto const bonds_view = p_handle.call_method("get_bonds_view", {});
    state.params.emplace_back(std::string{"bonds"}, pack(bonds_view));
#ifdef ESPRESSO_EXCLUSIONS
    auto const exclusions = p_handle.call_method("get_exclusions", {});
    state.params.emplace_back(std::string{"exclusions"}, pack(exclusions));
#endif // ESPRESSO_EXCLUSIONS
    state.params.emplace_back(std::string{"__cpt_sentinel"}, pack(None{}));
    return Utils::pack(state);
  });
 
  return Utils::pack(object_states);
}
 
void System::set_internal_state(std::string const &state) {
  auto const object_states = Utils::unpack<std::vector<std::string>>(state);
#ifdef ESPRESSO_EXCLUSIONS
  std::unordered_map<int, Variant> exclusions = {};
#endif // ESPRESSO_EXCLUSIONS
  std::unordered_map<int, Variant> bonds = {};
 
  for (auto const &packed_object : object_states) {
    auto state = Utils::unpack<ObjectState>(packed_object);
    VariantMap params = {};
    for (auto const &[name, packed_value] : state.params) {
      params[name] = unpack(packed_value, {});
    }
    auto const p_id = get_value<int>(params.at("id"));
    bonds[p_id] = params.extract("bonds").mapped();
#ifdef ESPRESSO_EXCLUSIONS
    exclusions[p_id] = params.extract("exclusions").mapped();
#endif // ESPRESSO_EXCLUSIONS
    params["__cell_structure"] = get_parameter("cell_system");
    context()->make_shared("Particles::ParticleHandle", params);
  }
 
  for (auto const p_id : get_particle_ids()) {
    auto p_obj = context()->make_shared(
        "Particles::ParticleHandle",
        {{"id", p_id}, {"__cell_structure", m_leaves->cell_system}});
    auto &p_handle = dynamic_cast<Particles::ParticleHandle &>(*p_obj);
    set_bonds(p_handle, bonds[p_id]);
#ifdef ESPRESSO_EXCLUSIONS
    set_exclusions(p_handle, exclusions[p_id]);
#endif // ESPRESSO_EXCLUSIONS
  }
}
 
} // namespace System
} // namespace ScriptInterface