diff --git a/src/core/EspressoSystemStandAlone.cpp b/src/core/EspressoSystemStandAlone.cpp index 5291191b49e..c4fbded6625 100644 --- a/src/core/EspressoSystemStandAlone.cpp +++ b/src/core/EspressoSystemStandAlone.cpp @@ -27,8 +27,6 @@ #include "integrate.hpp" #include "magnetostatics/dipoles.hpp" #include "system/System.hpp" -#include "virtual_sites.hpp" -#include "virtual_sites/VirtualSitesOff.hpp" #include @@ -46,9 +44,6 @@ EspressoSystemStandAlone::EspressoSystemStandAlone(int argc, char **argv) { Communication::init(mpi_env); // default-construct global state of the system -#ifdef VIRTUAL_SITES - set_virtual_sites(std::make_shared()); -#endif ::System::set_system(std::make_shared<::System::System>()); } diff --git a/src/core/event.cpp b/src/core/event.cpp index 9541c21235a..6232a3c8073 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -48,7 +48,7 @@ #include "particle_node.hpp" #include "system/System.hpp" #include "thermostat.hpp" -#include "virtual_sites.hpp" +#include "virtual_sites/relative.hpp" #include @@ -350,7 +350,7 @@ unsigned global_ghost_flags() { void update_dependent_particles() { #ifdef VIRTUAL_SITES - virtual_sites()->update(); + vs_relative_update_particles(cell_structure); cells_update_ghosts(global_ghost_flags()); #endif diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 351bc605ad5..0c6c26a7ed6 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -49,7 +49,7 @@ #include "system/System.hpp" #include "thermostat.hpp" #include "thermostats/langevin_inline.hpp" -#include "virtual_sites.hpp" +#include "virtual_sites/relative.hpp" #include @@ -255,7 +255,7 @@ void force_calc(CellStructure &cell_structure, double time_step, double kT) { // VIRTUAL_SITES distribute forces #ifdef VIRTUAL_SITES - virtual_sites()->back_transfer_forces_and_torques(); + vs_relative_back_transfer_forces_and_torques(cell_structure); #endif // Communication Step: ghost forces diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 0cc0bbc15b4..ca325c2b44c 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -55,7 +55,8 @@ #include "rotation.hpp" #include "signalhandling.hpp" #include "thermostat.hpp" -#include "virtual_sites.hpp" +#include "virtual_sites/lb_tracers.hpp" +#include "virtual_sites/relative.hpp" #include #include @@ -340,7 +341,7 @@ int integrate(int n_steps, int reuse_forces) { lb_lbcoupling_deactivate(); #ifdef VIRTUAL_SITES - virtual_sites()->update(); + vs_relative_update_particles(cell_structure); #endif // Communication step: distribute ghost positions @@ -418,7 +419,7 @@ int integrate(int n_steps, int reuse_forces) { #endif #ifdef VIRTUAL_SITES - virtual_sites()->update(); + vs_relative_update_particles(cell_structure); #endif if (cell_structure.get_resort_particles() >= Cells::RESORT_LOCAL) @@ -432,7 +433,7 @@ int integrate(int n_steps, int reuse_forces) { force_calc(cell_structure, time_step, temperature); #ifdef VIRTUAL_SITES - virtual_sites()->after_force_calc(time_step); + lb_tracers_add_particle_force_to_fluid(cell_structure, time_step); #endif integrator_step_2(particles, temperature); LeesEdwards::run_kernel(); @@ -490,7 +491,7 @@ int integrate(int n_steps, int reuse_forces) { } #ifdef VIRTUAL_SITES - virtual_sites()->after_lb_propagation(time_step); + lb_tracers_propagate(cell_structure, time_step); #endif #ifdef COLLISION_DETECTION @@ -523,7 +524,7 @@ int integrate(int n_steps, int reuse_forces) { #endif #ifdef VIRTUAL_SITES - virtual_sites()->update(); + vs_relative_update_particles(cell_structure); #endif // Verlet list statistics diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index d8b8ab8f734..2adbca7fc5d 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -38,7 +38,7 @@ #include "pressure_inline.hpp" #include "short_range_loop.hpp" #include "system/System.hpp" -#include "virtual_sites.hpp" +#include "virtual_sites/relative.hpp" #include "config/config.hpp" @@ -111,7 +111,8 @@ std::shared_ptr calculate_pressure() { #ifdef VIRTUAL_SITES if (!obs_pressure.virtual_sites.empty()) { - auto const vs_pressure = virtual_sites()->pressure_tensor(); + auto const vs_pressure = + vs_relative_pressure_tensor(cell_structure.local_particles()); boost::copy(flatten(vs_pressure), obs_pressure.virtual_sites.begin()); } #endif diff --git a/src/core/virtual_sites.cpp b/src/core/virtual_sites.cpp index cb973e5a7fb..c8c664405e4 100644 --- a/src/core/virtual_sites.cpp +++ b/src/core/virtual_sites.cpp @@ -40,23 +40,12 @@ #include #include -namespace { -std::shared_ptr m_virtual_sites; -} - -const std::shared_ptr &virtual_sites() { return m_virtual_sites; } - -void set_virtual_sites(std::shared_ptr const &v) { - m_virtual_sites = v; - recalc_forces = true; -} - #ifdef VIRTUAL_SITES_RELATIVE /** Calculate the rotation quaternion and distance between two particles */ std::tuple, double> -calculate_vs_relate_to_params(Particle const &p_vs, - Particle const &p_relate_to) { +calculate_vs_relate_to_params(Particle const &p_vs, Particle const &p_relate_to, + bool override_cutoff_check) { // get the distance between the particles auto d = ::box_geo.get_mi_vector(p_vs.pos(), p_relate_to.pos()); @@ -64,8 +53,7 @@ calculate_vs_relate_to_params(Particle const &p_vs, // than minimum global cutoff. If so, warn user. auto const dist = d.norm(); auto const min_global_cut = get_min_global_cut(); - if (dist > min_global_cut && n_nodes > 1 && - not virtual_sites()->get_override_cutoff_check()) { + if (dist > min_global_cut && n_nodes > 1 && not override_cutoff_check) { runtimeErrorMsg() << "Warning: The distance between virtual and non-virtual particle (" << dist << ") is larger than the minimum global cutoff (" diff --git a/src/core/virtual_sites.hpp b/src/core/virtual_sites.hpp index bf540b3beef..f0dc7d6e9ed 100644 --- a/src/core/virtual_sites.hpp +++ b/src/core/virtual_sites.hpp @@ -24,16 +24,6 @@ #ifdef VIRTUAL_SITES #include "Particle.hpp" -#include "virtual_sites/VirtualSites.hpp" - -#include - -/** @brief Get active virtual sites implementation */ -const std::shared_ptr &virtual_sites(); - -/** @brief Set active virtual sites implementation */ -void set_virtual_sites(std::shared_ptr const &v); - #ifdef VIRTUAL_SITES_RELATIVE #include @@ -42,7 +32,8 @@ void set_virtual_sites(std::shared_ptr const &v); std::tuple, double> calculate_vs_relate_to_params(Particle const &p_current, - Particle const &p_relate_to); + Particle const &p_relate_to, + bool override_cutoff_check = false); /** * @brief Setup a virtual site to track a real particle. diff --git a/src/core/virtual_sites/CMakeLists.txt b/src/core/virtual_sites/CMakeLists.txt index 13401606e58..a59c29840e2 100644 --- a/src/core/virtual_sites/CMakeLists.txt +++ b/src/core/virtual_sites/CMakeLists.txt @@ -17,7 +17,5 @@ # along with this program. If not, see . # -target_sources( - espresso_core - PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/VirtualSitesInertialessTracers.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/VirtualSitesRelative.cpp) +target_sources(espresso_core PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/lb_tracers.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/relative.cpp) diff --git a/src/core/virtual_sites/VirtualSites.hpp b/src/core/virtual_sites/VirtualSites.hpp deleted file mode 100644 index 0d172059d22..00000000000 --- a/src/core/virtual_sites/VirtualSites.hpp +++ /dev/null @@ -1,75 +0,0 @@ -/* - * Copyright (C) 2010-2022 The ESPResSo project - * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 - * Max-Planck-Institute for Polymer Research, Theory Group - * - * 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 . - */ -#ifndef VIRTUAL_SITES_VIRTUAL_SITES_HPP -#define VIRTUAL_SITES_VIRTUAL_SITES_HPP - -/** \file - * This file contains routine to handle virtual sites - * Virtual sites are like particles, but they will not be integrated. - * Step performed for virtual sites: - * - update virtual sites - * - calculate forces - * - distribute forces - * - move non-virtual particles - * - update virtual sites - */ - -#include "config/config.hpp" - -#ifdef VIRTUAL_SITES -#include -#include - -/** @brief Base class for virtual sites implementations */ -class VirtualSites { -public: - VirtualSites() = default; - virtual ~VirtualSites() = default; - - /** - * @brief Update positions and velocities of virtual sites. - */ - virtual void update() const {} - /** Back-transfer forces (and torques) to non-virtual particles. */ - virtual void back_transfer_forces_and_torques() const {} - /** @brief Called after force calculation (and before rattle/shake) */ - virtual void after_force_calc(double) {} - virtual void after_lb_propagation(double) {} - /** @brief Pressure contribution. */ - virtual Utils::Matrix pressure_tensor() const { return {}; } - /** @brief Enable/disable quaternion calculations for vs.*/ - void set_have_quaternion(const bool &have_quaternion) { - m_have_quaternion = have_quaternion; - } - bool have_quaternions() const { return m_have_quaternion; } - /** @brief Enable/disable override for the vs cutoff check */ - void set_override_cutoff_check(const bool &override_cutoff_check) { - m_override_cutoff_check = override_cutoff_check; - } - bool get_override_cutoff_check() const { return m_override_cutoff_check; } - -private: - bool m_have_quaternion = false; - bool m_override_cutoff_check = false; -}; - -#endif -#endif diff --git a/src/core/virtual_sites/VirtualSitesOff.hpp b/src/core/virtual_sites/VirtualSitesOff.hpp deleted file mode 100644 index 71bd1e9f71e..00000000000 --- a/src/core/virtual_sites/VirtualSitesOff.hpp +++ /dev/null @@ -1,30 +0,0 @@ -/* - * Copyright (C) 2010-2022 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 . - */ -#ifndef VIRTUAL_SITES_VIRTUAL_SITES_OFF_HPP -#define VIRTUAL_SITES_VIRTUAL_SITES_OFF_HPP - -#include "config/config.hpp" -#ifdef VIRTUAL_SITES -#include "VirtualSites.hpp" - -/** @brief Do-nothing virtual-sites implementation */ -class VirtualSitesOff : public VirtualSites {}; - -#endif -#endif diff --git a/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp b/src/core/virtual_sites/lb_tracers.cpp similarity index 79% rename from src/core/virtual_sites/VirtualSitesInertialessTracers.cpp rename to src/core/virtual_sites/lb_tracers.cpp index 9daad17461f..1b8235677d6 100644 --- a/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp +++ b/src/core/virtual_sites/lb_tracers.cpp @@ -20,8 +20,7 @@ #ifdef VIRTUAL_SITES_INERTIALESS_TRACERS -#include "VirtualSitesInertialessTracers.hpp" - +#include "PropagationModes.hpp" #include "cells.hpp" #include "errorhandling.hpp" #include "forces.hpp" @@ -40,12 +39,17 @@ static bool lb_active_check() { return true; } -void VirtualSitesInertialessTracers::after_force_calc(double time_step) { +static bool is_lb_tracer(const Particle &p) { + return (p.propagation() & PropagationMode::TRANS_LB_TRACER); +} + +void lb_tracers_add_particle_force_to_fluid(CellStructure &cell_struct, + double time_step) { auto const to_lb_units = (lattice_switch == ActiveLB::NONE) ? 0. : 1. / LB::get_agrid(); // Distribute summed-up forces from physical particles to ghosts - init_forces_ghosts(cell_structure.ghost_particles()); + init_forces_ghosts(cell_struct.ghost_particles()); cells_update_ghosts(Cells::DATA_PART_FORCE); // Set to store ghost particles (ids) that have already been coupled @@ -53,9 +57,9 @@ void VirtualSitesInertialessTracers::after_force_calc(double time_step) { // Apply particle forces to the LB fluid at particle positions // For physical particles, also set particle velocity = fluid velocity for (auto const &particle_range : - {cell_structure.local_particles(), cell_structure.ghost_particles()}) { + {cell_struct.local_particles(), cell_struct.ghost_particles()}) { for (auto const &p : particle_range) { - if (!p.is_virtual()) + if (!is_lb_tracer(p)) continue; if (!lb_active_check()) { return; @@ -69,16 +73,16 @@ void VirtualSitesInertialessTracers::after_force_calc(double time_step) { } // Clear ghost forces to avoid double counting later - init_forces_ghosts(cell_structure.ghost_particles()); + init_forces_ghosts(cell_struct.ghost_particles()); } -void VirtualSitesInertialessTracers::after_lb_propagation(double time_step) { +void lb_tracers_propagate(CellStructure &cell_struct, double time_step) { auto const to_md_units = (lattice_switch == ActiveLB::NONE) ? 0. : LB::get_lattice_speed(); // Advect particles - for (auto &p : cell_structure.local_particles()) { - if (!p.is_virtual()) + for (auto &p : cell_struct.local_particles()) { + if (!is_lb_tracer(p)) continue; if (!lb_active_check()) { return; @@ -91,7 +95,7 @@ void VirtualSitesInertialessTracers::after_lb_propagation(double time_step) { } // Verlet list update check if ((p.pos() - p.pos_at_last_verlet_update()).norm2() > skin * skin) { - cell_structure.set_resort_particles(Cells::RESORT_LOCAL); + cell_struct.set_resort_particles(Cells::RESORT_LOCAL); } } } diff --git a/src/core/virtual_sites/VirtualSitesInertialessTracers.hpp b/src/core/virtual_sites/lb_tracers.hpp similarity index 71% rename from src/core/virtual_sites/VirtualSitesInertialessTracers.hpp rename to src/core/virtual_sites/lb_tracers.hpp index 0fbc25bb1be..81b5a896f4e 100644 --- a/src/core/virtual_sites/VirtualSitesInertialessTracers.hpp +++ b/src/core/virtual_sites/lb_tracers.hpp @@ -23,15 +23,8 @@ #ifdef VIRTUAL_SITES_INERTIALESS_TRACERS -#include "VirtualSites.hpp" - -/** @brief Virtual sites which are advected with an lb fluid. Forces on them are - * instantaneously transferred to the fluid - */ -class VirtualSitesInertialessTracers : public VirtualSites { - void after_force_calc(double time_step) override; - void after_lb_propagation(double time_step) override; -}; - -#endif // VIRTUAL_SITES_INERTIALESS_TRACERS +void lb_tracers_add_particle_force_to_fluid(CellStructure &cell_struct, + double time_step); +void lb_tracers_propagate(CellStructure &cell_struct, double time_step); +#endif #endif diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/relative.cpp similarity index 69% rename from src/core/virtual_sites/VirtualSitesRelative.cpp rename to src/core/virtual_sites/relative.cpp index 7da96d1221b..c827830aac0 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/relative.cpp @@ -17,11 +17,11 @@ * along with this program. If not, see . */ -#include "VirtualSitesRelative.hpp" - +#include "config/config.hpp" #ifdef VIRTUAL_SITES_RELATIVE #include "Particle.hpp" +#include "PropagationModes.hpp" #include "cells.hpp" #include "errorhandling.hpp" #include "forces.hpp" @@ -78,9 +78,6 @@ static Utils::Vector3d velocity(Particle const &p_ref, Particle const &p_vs) { * @return Pointer to real particle. */ static Particle *get_reference_particle(Particle const &p) { - if (!p.is_virtual()) { - return nullptr; - } auto const &vs_rel = p.vs_relative(); if (vs_rel.to_particle_id == -1) { runtimeErrorMsg() << "Particle with id " << p.id() @@ -109,32 +106,52 @@ static auto constraint_stress(Particle const &p_ref, Particle const &p_vs) { return tensor_product(-p_vs.force(), connection_vector(p_ref, p_vs)); } -void VirtualSitesRelative::update() const { - cell_structure.ghosts_update(Cells::DATA_PART_POSITION | - Cells::DATA_PART_MOMENTUM); +static bool is_vs_relative_trans(const Particle &p) { + return p.propagation() & PropagationMode::TRANS_VS_RELATIVE; +}; +static bool is_vs_relative_rot(const Particle &p) { + return p.propagation() & PropagationMode::ROT_VS_RELATIVE; +}; + +static bool is_vs_relative(const Particle &p) { + return (is_vs_relative_trans(p) or is_vs_relative_rot(p)); +}; + +void vs_relative_update_particles(CellStructure &cell_struct) { + cell_struct.ghosts_update(Cells::DATA_PART_POSITION | + Cells::DATA_PART_MOMENTUM); - auto const particles = cell_structure.local_particles(); + auto const particles = cell_struct.local_particles(); for (auto &p : particles) { + if (!is_vs_relative(p)) { + continue; + } + auto const *p_ref_ptr = get_reference_particle(p); if (!p_ref_ptr) continue; auto const &p_ref = *p_ref_ptr; - p.image_box() = p_ref.image_box(); - p.pos() = p_ref.pos() + connection_vector(p_ref, p); - p.v() = velocity(p_ref, p); - - if (box_geo.type() == BoxType::LEES_EDWARDS) { - auto push = LeesEdwards::Push(box_geo); - push(p, -1); // includes a position fold - } else { - fold_position(p.pos(), p.image_box(), box_geo); + + // position update + if (is_vs_relative_trans(p)) { + p.image_box() = p_ref.image_box(); + p.pos() = p_ref.pos() + connection_vector(p_ref, p); + p.v() = velocity(p_ref, p); + + if (box_geo.type() == BoxType::LEES_EDWARDS) { + auto push = LeesEdwards::Push(box_geo); + push(p, -1); // includes a position fold + } else { + fold_position(p.pos(), p.image_box(), box_geo); + } } - if (have_quaternions()) + // Orientation update + if (is_vs_relative_rot(p)) { p.quat() = p_ref.quat() * p.vs_relative().quat; + } } - if (cell_structure.check_resort_required(particles, skin)) { cell_structure.set_resort_particles(Cells::RESORT_LOCAL); } @@ -142,33 +159,43 @@ void VirtualSitesRelative::update() const { // Distribute forces that have accumulated on virtual particles to the // associated real particles -void VirtualSitesRelative::back_transfer_forces_and_torques() const { - cell_structure.ghosts_reduce_forces(); +void vs_relative_back_transfer_forces_and_torques(CellStructure &cell_struct) { + cell_struct.ghosts_reduce_forces(); - init_forces_ghosts(cell_structure.ghost_particles()); + init_forces_ghosts(cell_struct.ghost_particles()); // Iterate over all the particles in the local cells - for (auto &p : cell_structure.local_particles()) { + for (auto &p : cell_struct.local_particles()) { + if (!is_vs_relative(p)) + continue; + auto *p_ref_ptr = get_reference_particle(p); if (!p_ref_ptr) continue; - // Add forces and torques auto &p_ref = *p_ref_ptr; - p_ref.force() += p.force(); - p_ref.torque() += - vector_product(connection_vector(p_ref, p), p.force()) + p.torque(); + if (is_vs_relative_trans(p)) { + p_ref.force() += p.force(); + p_ref.torque() += vector_product(connection_vector(p_ref, p), p.force()); + } + + if (is_vs_relative_rot(p)) { + p_ref.torque() += p.torque(); + } } } // Rigid body contribution to scalar pressure and pressure tensor -Utils::Matrix VirtualSitesRelative::pressure_tensor() const { +Utils::Matrix +vs_relative_pressure_tensor(const ParticleRange &particles) { Utils::Matrix pressure_tensor = {}; - for (auto const &p : cell_structure.local_particles()) { - auto const *p_ref_ptr = get_reference_particle(p); - if (p_ref_ptr) { - pressure_tensor += constraint_stress(*p_ref_ptr, p); + for (auto const &p : particles) { + if (is_vs_relative_trans(p)) { + auto const *p_ref_ptr = get_reference_particle(p); + if (p_ref_ptr) { + pressure_tensor += constraint_stress(*p_ref_ptr, p); + } } } diff --git a/src/core/virtual_sites/VirtualSitesRelative.hpp b/src/core/virtual_sites/relative.hpp similarity index 67% rename from src/core/virtual_sites/VirtualSitesRelative.hpp rename to src/core/virtual_sites/relative.hpp index e7050f9c348..006f3b7a5cb 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.hpp +++ b/src/core/virtual_sites/relative.hpp @@ -23,22 +23,14 @@ #include "config/config.hpp" #ifdef VIRTUAL_SITES_RELATIVE -#include "VirtualSites.hpp" - +#include "cells.hpp" #include #include -/** @brief Virtual sites implementation for rigid bodies */ -class VirtualSitesRelative : public VirtualSites { -public: - VirtualSitesRelative() = default; - /** @copydoc VirtualSites::update */ - void update() const override; - /** @copydoc VirtualSites::back_transfer_forces_and_torques */ - void back_transfer_forces_and_torques() const override; - /** @copydoc VirtualSites::pressure_tensor */ - Utils::Matrix pressure_tensor() const override; -}; +void vs_relative_update_particles(CellStructure &cell_struct); +void vs_relative_back_transfer_forces_and_torques(CellStructure &cell_struct); +Utils::Matrix +vs_relative_pressure_tensor(const ParticleRange &particles); #endif diff --git a/src/python/espressomd/particle_data.py b/src/python/espressomd/particle_data.py index 547fdae8cd5..9deec812876 100644 --- a/src/python/espressomd/particle_data.py +++ b/src/python/espressomd/particle_data.py @@ -26,6 +26,7 @@ from .utils import check_type_or_throw_except from .code_features import assert_features, has_features from .script_interface import script_interface_register, ScriptInterfaceHelper +from .propagation import Propagation import itertools @@ -552,6 +553,14 @@ def exclusions(self, p_ids): assert_features("EXCLUSIONS") self.call_method("set_exclusions", p_ids=p_ids) + @property + def propagation(self): + return Propagation(self.get_parameter("propagation")) + + @propagation.setter + def propagation(self, value): + self.set_parameter("propagation", int(value)) + def vs_auto_relate_to(self, rel_to): """ Setup this particle as virtual site relative to the particle @@ -929,6 +938,7 @@ def __setattr__(self, name, value): if name != "chunk_size" and name != "id_selection" and name not in particle_attributes: raise AttributeError( f"ParticleHandle does not have the attribute {name}.") + print("particle set attr", name, value) super().__setattr__(name, value) def to_dict(self): @@ -1088,6 +1098,8 @@ def add(self, *args, **kwargs): def _place_new_particle(self, p_dict): bonds = [] + if "propagation" in p_dict: p_dict["propagation"] = int( + p_dict["propagation"]) if "bonds" in p_dict: bonds = p_dict.pop("bonds") if nesting_level(bonds) == 1: diff --git a/src/python/espressomd/propagation.py b/src/python/espressomd/propagation.py index c1ffe6ab43f..2c736740399 100644 --- a/src/python/espressomd/propagation.py +++ b/src/python/espressomd/propagation.py @@ -5,11 +5,12 @@ class Propagation(IntFlag): NONE = 0 TRANS_SYSTEM_DEFAULT = 1 TRANS_LANGEVIN = 2 - TRANS_VS_RELATIVE = 4 - TRANS_LB_MOMENTUM_EXCHANGE = 8 - TRANS_LB_TRACER = 16 - TRANS_BROWNIAN = 32 - TRANS_STOKESIAN = 64 - ROT_LANGEVIN = 128 - ROT_VS_RELATIVE = 256 - ROT_BROWNIAN = 512 + TRANS_LANGEVIN_NPT = 4 + TRANS_VS_RELATIVE = 8 + TRANS_LB_MOMENTUM_EXCHANGE = 16 + TRANS_LB_TRACER = 32 + TRANS_BROWNIAN = 64 + TRANS_STOKESIAN = 128 + ROT_LANGEVIN = 256 + ROT_VS_RELATIVE = 512 + ROT_BROWNIAN = 1024 diff --git a/src/python/espressomd/script_interface.pyx b/src/python/espressomd/script_interface.pyx index d663a96adb4..37637a4dd58 100644 --- a/src/python/espressomd/script_interface.pyx +++ b/src/python/espressomd/script_interface.pyx @@ -447,20 +447,21 @@ class ScriptInterfaceHelper(PScriptInterface): return list(self.__dict__.keys()) + self._valid_parameters() def __getattr__(self, attr): - if attr in self._valid_parameters(): - return self.get_parameter(attr) - if attr in self.__dict__: return self.__dict__[attr] + if attr in self._valid_parameters(): + return self.get_parameter(attr) + raise AttributeError( f"Object '{self.__class__.__name__}' has no attribute '{attr}'") def __setattr__(self, attr, value): - if attr in self._valid_parameters(): - self.set_params(**{attr: value}) - else: + if attr in self.__class__.__dict__ or ( + not attr in self._valid_parameters()): super().__setattr__(attr, value) + else: + self.set_params(**{attr: value}) def __delattr__(self, attr): if attr in self._valid_parameters(): diff --git a/src/python/espressomd/system.py b/src/python/espressomd/system.py index a8bc698a247..fc777fad098 100644 --- a/src/python/espressomd/system.py +++ b/src/python/espressomd/system.py @@ -38,7 +38,6 @@ from . import lees_edwards from . import particle_data from . import thermostat -from . import virtual_sites from .code_features import has_features, assert_features from . import utils @@ -183,8 +182,6 @@ def __init__(self, **kwargs): setable_properties = ["box_l", "min_global_cut", "periodicity", "time", "time_step", "force_cap", "max_oif_objects"] - if has_features("VIRTUAL_SITES"): - setable_properties.append("_active_virtual_sites_handle") self.call_method("set_system_handle", obj=_System(**kwargs)) self.integrator = integrate.IntegratorHandle() @@ -221,9 +218,6 @@ def __init__(self, **kwargs): self.non_bonded_inter = interactions.NonBondedInteractions() self.part = particle_data.ParticleList() self.thermostat = thermostat.Thermostat() - if has_features("VIRTUAL_SITES"): - self._active_virtual_sites_handle = virtual_sites.ActiveVirtualSitesHandle( - implementation=virtual_sites.VirtualSitesOff()) # lock class self.call_method("lock_system_creation") @@ -250,8 +244,6 @@ def _restore_object(cls, so_callback, so_callback_args, state): def __getstate__(self): checkpointable_properties = ["integrator"] - if has_features("VIRTUAL_SITES"): - checkpointable_properties.append("_active_virtual_sites_handle") checkpointable_properties += [ "non_bonded_inter", "bonded_inter", "cell_system", "lees_edwards", "part", "analysis", "auto_update_accumulators", @@ -346,23 +338,6 @@ def max_cut_bonded(self): return self.cell_system.max_cut_bonded @property - def virtual_sites(self): - """ - Set the virtual site implementation. - - Requires feature ``VIRTUAL_SITES``. - - Type: :obj:`espressomd.virtual_sites.ActiveVirtualSitesHandle` - - """ - assert_features("VIRTUAL_SITES") - return self._active_virtual_sites_handle.implementation - - @virtual_sites.setter - def virtual_sites(self, value): - assert_features("VIRTUAL_SITES") - self._active_virtual_sites_handle.implementation = value - def change_volume_and_rescale_particles(self, d_new, dir="xyz"): """Change box size and rescale particle coordinates. diff --git a/testsuite/python/virtual_sites_relative.py b/testsuite/python/virtual_sites_relative.py index 6a575542d83..f4389d4e820 100644 --- a/testsuite/python/virtual_sites_relative.py +++ b/testsuite/python/virtual_sites_relative.py @@ -19,7 +19,7 @@ import unittest as ut import unittest_decorators as utx import espressomd -import espressomd.virtual_sites +from espressomd.propagation import Propagation import numpy as np import tests_common @@ -32,6 +32,7 @@ class VirtualSites(ut.TestCase): def setUp(self): self.system.box_l = [10.0, 10.0, 10.0] + self.system.cell_system.set_regular_decomposition( use_verlet_lists=True) @@ -40,7 +41,6 @@ def tearDown(self): self.system.thermostat.turn_off() self.system.integrator.set_vv() self.system.non_bonded_inter[0, 0].lennard_jones.deactivate() - self.system.virtual_sites = espressomd.virtual_sites.VirtualSitesOff() def multiply_quaternions(self, a, b): return np.array( @@ -58,8 +58,7 @@ def director_from_quaternion(self, quat): def verify_vs(self, vs, verify_velocity=True): """Verify virtual site position and velocity.""" - self.assertTrue(vs.virtual) - + self.assertTrue(vs.propagation & Propagation.TRANS_VS_RELATIVE) vs_r = vs.vs_relative # Get related particle @@ -81,52 +80,34 @@ def verify_vs(self, vs, verify_velocity=True): v_d - vs_r[1] * self.director_from_quaternion( self.multiply_quaternions(rel.quat, vs_r[2]))), 1E-6) - def test_aa_method_switching(self): - # Virtual sites should be disabled by default - self.assertIsInstance( - self.system.virtual_sites, - espressomd.virtual_sites.VirtualSitesOff) - self.assertFalse(self.system.virtual_sites.have_quaternion) - self.assertFalse(self.system.virtual_sites.override_cutoff_check) - - # Set properties - self.system.virtual_sites.have_quaternion = True - self.system.virtual_sites.override_cutoff_check = True - self.assertTrue(self.system.virtual_sites.have_quaternion) - self.assertTrue(self.system.virtual_sites.override_cutoff_check) - - # Switch implementation - self.system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative() - self.assertIsInstance( - self.system.virtual_sites, - espressomd.virtual_sites.VirtualSitesRelative) - self.assertFalse(self.system.virtual_sites.have_quaternion) - self.assertFalse(self.system.virtual_sites.override_cutoff_check) - def test_vs_quat(self): self.system.time_step = 0.01 self.system.min_global_cut = 0.23 # First check that quaternion of virtual particle is unchanged if - # have_quaternion is false. - self.system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative( - have_quaternion=False) - self.assertFalse(self.system.virtual_sites.have_quaternion) + # The rotation propagaoitn is not set to vs_relative p1 = self.system.part.add(pos=[1, 1, 1], rotation=3 * [True], omega_lab=[1, 1, 1]) - p2 = self.system.part.add(pos=[1, 1, 1], rotation=3 * [True]) + p2 = self.system.part.add( + pos=[ + 1, + 1, + 1], + rotation=3 * + [True], + propagation=Propagation.TRANS_VS_RELATIVE) p2.vs_auto_relate_to(p1) np.testing.assert_array_equal(np.copy(p2.quat), [1, 0, 0, 0]) self.system.integrator.run(1) np.testing.assert_array_equal(np.copy(p2.quat), [1, 0, 0, 0]) # Now check that quaternion of the virtual particle gets updated. - self.system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative( - have_quaternion=True) + p2.propagation = p2.propagation + Propagation.ROT_VS_RELATIVE self.system.integrator.run(1) self.assertRaises(AssertionError, np.testing.assert_array_equal, np.copy(p2.quat), [1, 0, 0, 0]) # co-aligned case p2.vs_quat = (1, 0, 0, 0) + p2.propagation = Propagation.TRANS_VS_RELATIVE + Propagation.ROT_VS_RELATIVE self.system.integrator.run(1) np.testing.assert_allclose( np.copy(p2.director), np.copy(p1.director), atol=1E-12) @@ -149,7 +130,6 @@ def test_vs_quat(self): def test_vs_exceptions(self): system = self.system - system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative() system.time_step = 0.01 system.cell_system.skin = 0.1 system.min_global_cut = 0.1 @@ -159,7 +139,7 @@ def test_vs_exceptions(self): p4 = system.part.add(pos=[1.0, 1.0, 1.0], rotation=3 * [True], id=4) # dangling virtual sites are not allowed with self.assertRaisesRegex(Exception, "Particle with id 4 is a dangling virtual site"): - p4.virtual = True + p4.propagation = Propagation.TRANS_VS_RELATIVE self.assertEqual(p4.vs_relative[0], -1) system.integrator.run(0, recalc_forces=True) p4.remove() @@ -174,19 +154,18 @@ def test_vs_exceptions(self): # relating to a deleted particle is not allowed with self.assertRaisesRegex(Exception, "No real particle with id 3 for virtual site with id 2"): p2.vs_auto_relate_to(p3) + p2.propagation = Propagation.TRANS_VS_RELATIVE p3.remove() system.integrator.run(0, recalc_forces=True) if system.cell_system.get_state()["n_nodes"] > 1: with self.assertRaisesRegex(Exception, r"The distance between virtual and non-virtual particle \([0-9\.]+\) is larger than the minimum global cutoff"): p2.vs_auto_relate_to(p1) # If overridden this check should not raise an exception - system.virtual_sites.override_cutoff_check = True p2.vs_auto_relate_to(p1) def test_pos_vel_forces(self): system = self.system system.cell_system.skin = 0.3 - system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative() system.time_step = 0.004 # Check setting of min_global_cut @@ -200,7 +179,10 @@ def test_pos_vel_forces(self): pos3 = (0.3, 0.5, 0.4) pos4 = (0.5, 0.5, 0.5) for pos in (pos2, pos3, pos4): - p = system.part.add(rotation=3 * [True], pos=pos) + p = system.part.add( + rotation=3 * [True], + pos=pos, + propagation=Propagation.TRANS_VS_RELATIVE) p.vs_auto_relate_to(p1) # Was the particle made virtual self.assertTrue(p.virtual) @@ -260,7 +242,6 @@ def run_test_lj(self): get lost or are outdated in the short range loop. """ system = self.system - system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative() # Parameters n = 40 phi = 0.6 @@ -291,12 +272,12 @@ def run_test_lj(self): id=3 * i + 1, pos=p3i.pos + p3i.director / 2., - type=0) + type=0, propagation=Propagation.TRANS_VS_RELATIVE) p3ip2 = system.part.add(rotation=3 * [True], id=3 * i + 2, pos=p3i.pos - p3i.director / 2., - type=0) + type=0, propagation=Propagation.TRANS_VS_RELATIVE) p3ip1.vs_auto_relate_to(p3i.id) self.verify_vs(p3ip1, verify_velocity=False) p3ip2.vs_auto_relate_to(p3i.id) @@ -370,7 +351,6 @@ def test_zz_pressure_tensor(self): system.cell_system.skin = 0.1 system.min_global_cut = 0.2 # Should not have a pressure - system.virtual_sites = espressomd.virtual_sites.VirtualSitesOff() pressure_tensor_vs = system.analysis.pressure_tensor()[ "virtual_sites", 0] p_vs = system.analysis.pressure()["virtual_sites", 0] @@ -378,10 +358,13 @@ def test_zz_pressure_tensor(self): np.testing.assert_allclose(p_vs, 0., atol=1e-10) # vs relative contrib - system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative() p0 = system.part.add(id=0, pos=(0.0, 0.0, 0.0)) - p1 = system.part.add(id=1, pos=(0.1, 0.1, 0.1), ext_force=(1, 2, 3)) - p2 = system.part.add(id=2, pos=(0.1, 0.0, 0.0), ext_force=(-1, 0, 0)) + p1 = system.part.add( + id=1, pos=( + 0.1, 0.1, 0.1), ext_force=( + 1, 2, 3), propagation=Propagation.TRANS_VS_RELATIVE) + p2 = system.part.add(id=2, pos=( + 0.1, 0.0, 0.0), ext_force=(-1, 0, 0), propagation=Propagation.TRANS_VS_RELATIVE) p1.vs_auto_relate_to(p0) p2.vs_auto_relate_to(p0) system.integrator.run(0, recalc_forces=True)