From 713eb08f451a9e70edff4915604267ecc4d069d3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 8 Dec 2023 01:21:52 +0100 Subject: [PATCH] Remove virtual flag --- samples/immersed_boundary/addSoft.py | 2 +- .../sampleImmersedBoundary.py | 2 +- src/core/Particle.hpp | 29 +++-------- src/core/ParticleRange.hpp | 13 ++--- src/core/PropagationModes.hpp | 6 +++ src/core/collision.cpp | 1 - src/core/forces.cpp | 22 +++++--- src/core/ghosts.cpp | 6 +-- src/core/integrate.cpp | 4 +- src/core/integrate.hpp | 3 -- src/core/integrators/brownian_inline.hpp | 15 ++---- src/core/integrators/steepest_descent.cpp | 23 ++++----- .../integrators/stokesian_dynamics_inline.hpp | 6 +-- .../integrators/velocity_verlet_inline.hpp | 13 +---- src/core/integrators/velocity_verlet_npt.cpp | 9 ---- src/core/lb/particle_coupling.cpp | 9 ++-- src/core/lb/particle_coupling.hpp | 25 ++++++--- src/core/thermostat.cpp | 11 ---- src/core/thermostat.hpp | 12 ++--- src/core/thermostats/brownian_inline.hpp | 13 +---- src/core/thermostats/langevin_inline.hpp | 9 +--- .../unit_tests/lb_particle_coupling_test.cpp | 18 +++---- src/core/virtual_sites.hpp | 2 +- src/core/virtual_sites/lb_tracers.cpp | 8 +-- src/python/espressomd/particle_data.py | 3 +- src/python/espressomd/thermostat.pxd | 2 - src/python/espressomd/thermostat.pyx | 35 ++----------- .../particle_data/ParticleHandle.cpp | 26 ++++------ testsuite/python/analyze_chains.py | 5 +- testsuite/python/analyze_mass_related.py | 18 ++++--- testsuite/python/brownian_dynamics.py | 21 +++----- testsuite/python/collision_detection.py | 16 +++--- testsuite/python/field_test.py | 5 +- testsuite/python/langevin_thermostat.py | 16 ++---- testsuite/python/lb_thermo_virtual.py | 51 ++++--------------- testsuite/python/lees_edwards.py | 2 +- testsuite/python/npt_thermostat.py | 7 +-- testsuite/python/observables.py | 9 ++-- testsuite/python/particle.py | 7 ++- testsuite/python/rotate_system.py | 8 +-- testsuite/python/test_checkpoint.py | 2 - testsuite/python/virtual_sites_relative.py | 2 +- .../python/virtual_sites_tracers_common.py | 10 ++-- 43 files changed, 176 insertions(+), 330 deletions(-) diff --git a/samples/immersed_boundary/addSoft.py b/samples/immersed_boundary/addSoft.py index e851bed8d45..c0baed64073 100644 --- a/samples/immersed_boundary/addSoft.py +++ b/samples/immersed_boundary/addSoft.py @@ -31,7 +31,7 @@ def AddSoft(system, comX, comY, comZ, k1, k2): X = float(line[0]) + comX Y = float(line[1]) + comY Z = float(line[2]) + comZ - system.part.add(id=i, pos=[X, Y, Z], virtual=True) + system.part.add(id=i, pos=[X, Y, Z]) # triangles import espressomd.interactions diff --git a/samples/immersed_boundary/sampleImmersedBoundary.py b/samples/immersed_boundary/sampleImmersedBoundary.py index c675a9f2eee..b8ac18771a1 100644 --- a/samples/immersed_boundary/sampleImmersedBoundary.py +++ b/samples/immersed_boundary/sampleImmersedBoundary.py @@ -77,7 +77,7 @@ agrid=1, density=1, kinematic_viscosity=1, tau=system.time_step, ext_force_density=[force, 0, 0]) system.lb = lbf -system.thermostat.set_lb(LB_fluid=lbf, gamma=1.0, act_on_virtual=False) +system.thermostat.set_lb(LB_fluid=lbf, gamma=1.0) # Setup boundaries wall_shapes = [None] * 2 diff --git a/src/core/Particle.hpp b/src/core/Particle.hpp index f2f5bb3f64c..8e275b0ecb6 100644 --- a/src/core/Particle.hpp +++ b/src/core/Particle.hpp @@ -16,8 +16,8 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef ESPRESSO_CORE_PARTICLE_HPP -#define ESPRESSO_CORE_PARTICLE_HPP + +#pragma once #include "config/config.hpp" @@ -79,13 +79,6 @@ struct ParticleProperties { /** determines which propagation schemes should be applied to the particle **/ int propagation = PropagationMode::TRANS_SYSTEM_DEFAULT; -#ifdef VIRTUAL_SITES - /** is particle virtual */ - bool is_virtual = false; -#else // VIRTUAL_SITES - static constexpr bool is_virtual = false; -#endif // VIRTUAL_SITES - #ifdef ROTATION /** Bitfield for the particle axes of rotation. * Values: @@ -234,12 +227,9 @@ struct ParticleProperties { #ifdef DIPOLE_FIELD_TRACKING ar &dip_fld; #endif -#ifdef VIRTUAL_SITES - ar &is_virtual; #ifdef VIRTUAL_SITES_RELATIVE ar &vs_relative; #endif -#endif // VIRTUAL_SITES #ifdef THERMOSTAT_PER_PARTICLE ar γ @@ -524,18 +514,15 @@ struct Particle { // NOLINT(bugprone-exception-escape) auto const &mu_E() const { return p.mu_E; } auto &mu_E() { return p.mu_E; } #endif -#ifdef VIRTUAL_SITES - auto &virtual_flag() { return p.is_virtual; } - auto const &virtual_flag() const { return p.is_virtual; } - auto is_virtual() const { return p.is_virtual; } - void set_virtual(bool const virt_flag) { p.is_virtual = virt_flag; } + auto is_virtual() const { + return (p.propagation & (PropagationMode::TRANS_VS_RELATIVE | + PropagationMode::ROT_VS_RELATIVE | + PropagationMode::TRANS_LB_TRACER)) != 0; + } #ifdef VIRTUAL_SITES_RELATIVE auto const &vs_relative() const { return p.vs_relative; } auto &vs_relative() { return p.vs_relative; } #endif // VIRTUAL_SITES_RELATIVE -#else - constexpr auto is_virtual() const { return p.is_virtual; } -#endif #ifdef THERMOSTAT_PER_PARTICLE auto const &gamma() const { return p.gamma; } auto &gamma() { return p.gamma; } @@ -633,5 +620,3 @@ BOOST_IS_BITWISE_SERIALIZABLE(ParticleRattle) #ifdef VIRTUAL_SITES_RELATIVE BOOST_IS_BITWISE_SERIALIZABLE(decltype(ParticleProperties::vs_relative)) #endif - -#endif diff --git a/src/core/ParticleRange.hpp b/src/core/ParticleRange.hpp index fa5af366b4d..b050c60cb7f 100644 --- a/src/core/ParticleRange.hpp +++ b/src/core/ParticleRange.hpp @@ -16,8 +16,8 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef ESPRESSO_SRC_CORE_PARTICLE_RANGE_HPP -#define ESPRESSO_SRC_CORE_PARTICLE_RANGE_HPP + +#pragma once #include "CellParticleIterator.hpp" #include "Particle.hpp" @@ -51,14 +51,11 @@ class ParticleRange : public boost::iterator_range { template ParticleRangeFiltered filter(Predicate pred) const { - return {boost::make_filter_iterator(PropagationPredicate(pred), - begin(), end()), - boost::make_filter_iterator(PropagationPredicate(pred), - end(), end())}; + auto const predicate = PropagationPredicate(pred); + return {boost::make_filter_iterator(predicate, begin(), end()), + boost::make_filter_iterator(predicate, end(), end())}; } private: base_type::difference_type mutable m_size = -1; }; - -#endif diff --git a/src/core/PropagationModes.hpp b/src/core/PropagationModes.hpp index b9dae892dad..c3270d84fb3 100644 --- a/src/core/PropagationModes.hpp +++ b/src/core/PropagationModes.hpp @@ -19,6 +19,9 @@ #pragma once +/** + * @brief Bitfield for propagation modes. + */ enum PropagationMode { NONE = 0, TRANS_SYSTEM_DEFAULT = 1 << 0, @@ -34,6 +37,9 @@ enum PropagationMode { ROT_BROWNIAN = 1 << 10 }; +/** + * @brief Check allowlist of valid propagation modes combinations. + */ inline bool is_valid_propagation_combination(int propagation) { if (propagation == NONE) return true; diff --git a/src/core/collision.cpp b/src/core/collision.cpp index 88941d72a6c..c67aea601c5 100644 --- a/src/core/collision.cpp +++ b/src/core/collision.cpp @@ -356,7 +356,6 @@ static void place_vs_and_relate_to_particle(CellStructure &cell_structure, auto p_vs = cell_structure.add_particle(std::move(new_part)); vs_relate_to(*p_vs, get_part(cell_structure, relate_to), box_geo, min_global_cut); - p_vs->set_virtual(true); p_vs->type() = collision_params.vs_particle_type; } diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 73e56d70200..2af1c03a9b2 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -96,14 +96,22 @@ inline ParticleForce thermostat_force(Particle const &p, double time_step, return {}; } + auto const prop = p.propagation(); + auto const translation = prop & (PropagationMode::TRANS_LANGEVIN | + PropagationMode::TRANS_SYSTEM_DEFAULT); #ifdef ROTATION - return {friction_thermo_langevin(langevin, p, time_step, kT), - p.can_rotate() ? convert_vector_body_to_space( - p, friction_thermo_langevin_rotation( - langevin, p, time_step, kT)) - : Utils::Vector3d{}}; + auto const rotation = prop & (PropagationMode::ROT_LANGEVIN | + PropagationMode::TRANS_SYSTEM_DEFAULT) and + p.can_rotate(); + return {translation ? friction_thermo_langevin(langevin, p, time_step, kT) + : Utils::Vector3d{}, + rotation ? convert_vector_body_to_space( + p, friction_thermo_langevin_rotation(langevin, p, + time_step, kT)) + : Utils::Vector3d{}}; #else - return friction_thermo_langevin(langevin, p, time_step, kT); + return translation ? friction_thermo_langevin(langevin, p, time_step, kT) + : Utils::Vector3d{}; #endif } @@ -254,7 +262,7 @@ void System::System::calculate_forces(double kT) { immersed_boundaries.volume_conservation(*cell_structure); if (lb.is_solver_set()) { - LB::couple_particles(thermo_virtual, particles, ghost_particles, time_step); + LB::couple_particles(particles, ghost_particles, time_step); } #ifdef CUDA diff --git a/src/core/ghosts.cpp b/src/core/ghosts.cpp index 3c439957c43..c513656ffb9 100644 --- a/src/core/ghosts.cpp +++ b/src/core/ghosts.cpp @@ -126,11 +126,7 @@ serialize_and_reduce(Archive &ar, Particle &p, unsigned int data_parts, BoxGeometry const &box_geo, Utils::Vector3d const *ghost_shift) { if (data_parts & GHOSTTRANS_PROPRTS) { - ar &p.id() & p.mol_id() & p.type(); -#ifdef VIRTUAL_SITES - ar &p.virtual_flag(); -#endif - ar &p.propagation(); + ar &p.id() & p.mol_id() & p.type() & p.propagation(); #ifdef ROTATION ar &p.rotation(); #ifdef ROTATIONAL_INERTIA diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 47e16aa0ebd..6cf8e85b234 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -332,7 +332,7 @@ static bool integrator_step_1(ParticleRange const &particles, double kT, } #endif // STOKESIAN_DYNAMICS - increment_sim_time(time_step); + sim_time += time_step; return false; } @@ -671,8 +671,6 @@ double get_time_step() { return System::get_system().get_time_step(); } double get_sim_time() { return sim_time; } -void increment_sim_time(double amount) { sim_time += amount; } - void set_time(double value) { ::sim_time = value; ::recalc_forces = true; diff --git a/src/core/integrate.hpp b/src/core/integrate.hpp index 9fc107d4346..8c6f5dcef62 100644 --- a/src/core/integrate.hpp +++ b/src/core/integrate.hpp @@ -85,9 +85,6 @@ double get_time_step(); /** Get simulation time */ double get_sim_time(); -/** Increase simulation time (only on head node) */ -void increment_sim_time(double amount); - /** @brief Set the simulation time. */ void set_time(double value); diff --git a/src/core/integrators/brownian_inline.hpp b/src/core/integrators/brownian_inline.hpp index 625d04f0a66..e6575c9a258 100644 --- a/src/core/integrators/brownian_inline.hpp +++ b/src/core/integrators/brownian_inline.hpp @@ -19,8 +19,7 @@ * along with this program. If not, see . */ -#ifndef INTEGRATORS_BROWNIAN_INLINE_HPP -#define INTEGRATORS_BROWNIAN_INLINE_HPP +#pragma once #include "config/config.hpp" @@ -33,13 +32,10 @@ inline void brownian_dynamics_propagator(BrownianThermostat const &brownian, Particle &p, double time_step, double kT) { - // Don't propagate translational degrees of freedom of vs - if (!p.is_virtual() or thermo_virtual) { - p.pos() += bd_drag(brownian.gamma, p, time_step); - p.v() = bd_drag_vel(brownian.gamma, p); - p.pos() += bd_random_walk(brownian, p, time_step, kT); - p.v() += bd_random_walk_vel(brownian, p); - } + p.pos() += bd_drag(brownian.gamma, p, time_step); + p.v() = bd_drag_vel(brownian.gamma, p); + p.pos() += bd_random_walk(brownian, p, time_step, kT); + p.v() += bd_random_walk_vel(brownian, p); } #ifdef ROTATION @@ -55,4 +51,3 @@ inline void brownian_dynamics_rotator(BrownianThermostat const &brownian, p.omega() += bd_random_walk_vel_rot(brownian, p); } #endif // ROTATION -#endif // INTEGRATORS_BROWNIAN_INLINE_HPP diff --git a/src/core/integrators/steepest_descent.cpp b/src/core/integrators/steepest_descent.cpp index 2370da4b02e..0e7b9b3540b 100644 --- a/src/core/integrators/steepest_descent.cpp +++ b/src/core/integrators/steepest_descent.cpp @@ -55,19 +55,16 @@ bool steepest_descent_step(const ParticleRange &particles) { for (unsigned int j = 0; j < 3; j++) { // Skip, if coordinate is fixed if (!p.is_fixed_along(j)) { - // Skip positional increments of virtual particles - if (!p.is_virtual()) { - // Square of force on particle - f += Utils::sqr(p.force()[j]); - - // Positional increment, crop to maximum allowed by user - auto const dp = - std::clamp(params.gamma * p.force()[j], -params.max_displacement, - params.max_displacement); - - // Move particle - p.pos()[j] += dp; - } + // Square of force on particle + f += Utils::sqr(p.force()[j]); + + // Positional increment, crop to maximum allowed by user + auto const dp = + std::clamp(params.gamma * p.force()[j], -params.max_displacement, + params.max_displacement); + + // Move particle + p.pos()[j] += dp; } } #ifdef ROTATION diff --git a/src/core/integrators/stokesian_dynamics_inline.hpp b/src/core/integrators/stokesian_dynamics_inline.hpp index 760e79594f7..6516705e8b3 100644 --- a/src/core/integrators/stokesian_dynamics_inline.hpp +++ b/src/core/integrators/stokesian_dynamics_inline.hpp @@ -16,12 +16,13 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef STOKESIAN_DYNAMICS_INLINE_HPP -#define STOKESIAN_DYNAMICS_INLINE_HPP + +#pragma once #include "config/config.hpp" #ifdef STOKESIAN_DYNAMICS + #include "ParticleRange.hpp" #include "communication.hpp" #include "integrate.hpp" @@ -62,4 +63,3 @@ inline void stokesian_dynamics_step_1(const ParticleRange &particles, } #endif // STOKESIAN_DYNAMICS -#endif diff --git a/src/core/integrators/velocity_verlet_inline.hpp b/src/core/integrators/velocity_verlet_inline.hpp index 26eba7648f0..6391c4fde2a 100644 --- a/src/core/integrators/velocity_verlet_inline.hpp +++ b/src/core/integrators/velocity_verlet_inline.hpp @@ -16,8 +16,8 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef INTEGRATORS_VELOCITY_VERLET_HPP -#define INTEGRATORS_VELOCITY_VERLET_HPP + +#pragma once #include "config/config.hpp" @@ -34,9 +34,6 @@ */ inline void velocity_verlet_propagate_vel_pos_par(Particle &p, double time_step) { - // Don't propagate translational degrees of freedom of vs - if (p.is_virtual()) - return; for (unsigned int j = 0; j < 3; j++) { if (!p.is_fixed_along(j)) { /* Propagate velocities: v(t+0.5*dt) = v(t) + 0.5 * dt * a(t) */ @@ -54,10 +51,6 @@ inline void velocity_verlet_propagate_vel_pos_par(Particle &p, */ inline void velocity_verlet_propagate_vel_final_par(Particle &p, double time_step) { - // Virtual sites are not propagated during integration - if (p.is_virtual()) - return; - for (unsigned int j = 0; j < 3; j++) { if (!p.is_fixed_along(j)) { /* Propagate velocity: v(t+dt) = v(t+0.5*dt) + 0.5*dt * a(t+dt) */ @@ -65,5 +58,3 @@ inline void velocity_verlet_propagate_vel_final_par(Particle &p, } } } - -#endif diff --git a/src/core/integrators/velocity_verlet_npt.cpp b/src/core/integrators/velocity_verlet_npt.cpp index d8cfdf6c8ac..36b2a20f0cd 100644 --- a/src/core/integrators/velocity_verlet_npt.cpp +++ b/src/core/integrators/velocity_verlet_npt.cpp @@ -51,9 +51,6 @@ velocity_verlet_npt_propagate_vel_final(ParticleRangeNPT const &particles, nptiso.p_vel = {}; for (auto &p : particles) { - // Virtual sites are not propagated during integration - if (p.is_virtual()) - continue; auto const noise = friction_therm0_nptiso<2>(npt_iso, p.v(), p.id()); for (unsigned int j = 0; j < 3; j++) { if (!p.is_fixed_along(j)) { @@ -127,8 +124,6 @@ static void velocity_verlet_npt_propagate_pos(ParticleRangeNPT const &particles, /* propagate positions while rescaling positions and velocities */ for (auto &p : particles) { - if (p.is_virtual()) - continue; for (unsigned int j = 0; j < 3; j++) { if (!p.is_fixed_along(j)) { if (nptiso.geometry & ::nptgeom_dir[j]) { @@ -170,10 +165,6 @@ static void velocity_verlet_npt_propagate_vel(ParticleRangeNPT const &particles, nptiso.p_vel = {}; for (auto &p : particles) { - - // Don't propagate translational degrees of freedom of vs - if (p.is_virtual()) - continue; for (unsigned int j = 0; j < 3; j++) { if (!p.is_fixed_along(j)) { auto const noise = friction_therm0_nptiso<1>(npt_iso, p.v(), p.id()); diff --git a/src/core/lb/particle_coupling.cpp b/src/core/lb/particle_coupling.cpp index 416522eae59..10dcc8be663 100644 --- a/src/core/lb/particle_coupling.cpp +++ b/src/core/lb/particle_coupling.cpp @@ -241,9 +241,6 @@ Utils::Vector3d ParticleCoupling::get_noise_term(Particle const &p) const { } void ParticleCoupling::kernel(Particle &p) { - if (p.is_virtual() and not m_couple_virtual) - return; - auto const agrid = m_lb.get_agrid(); auto const &box_geo = *System::get_system().box_geo; @@ -303,7 +300,7 @@ static void lb_coupling_sanity_checks(Particle const &p) { } #endif -void couple_particles(bool couple_virtual, ParticleRange const &real_particles, +void couple_particles(ParticleRange const &real_particles, ParticleRange const &ghost_particles, double time_step) { #ifdef CALIPER CALI_CXX_MARK_FUNCTION; @@ -311,11 +308,11 @@ void couple_particles(bool couple_virtual, ParticleRange const &real_particles, if (lb_particle_coupling.couple_to_md) { auto &lb = System::get_system().lb; if (lb.is_solver_set()) { - ParticleCoupling coupling{lb, couple_virtual, time_step}; + ParticleCoupling coupling{lb, time_step}; CouplingBookkeeping bookkeeping{}; for (auto const &particle_range : {real_particles, ghost_particles}) { for (auto &p : particle_range) { - if (bookkeeping.should_be_coupled(p)) { + if (not LB::is_tracer(p) and bookkeeping.should_be_coupled(p)) { #if defined(THERMOSTAT_PER_PARTICLE) and defined(PARTICLE_ANISOTROPY) lb_coupling_sanity_checks(p); #endif diff --git a/src/core/lb/particle_coupling.hpp b/src/core/lb/particle_coupling.hpp index fce6d7ceb52..80549721fa8 100644 --- a/src/core/lb/particle_coupling.hpp +++ b/src/core/lb/particle_coupling.hpp @@ -129,21 +129,18 @@ extern LB::ParticleCouplingConfig lb_particle_coupling; namespace LB { /** @brief Calculate particle-lattice interactions. */ -void couple_particles(bool couple_virtual, ParticleRange const &real_particles, +void couple_particles(ParticleRange const &real_particles, ParticleRange const &ghost_particles, double time_step); class ParticleCoupling { LB::Solver &m_lb; - bool m_couple_virtual; bool m_thermalized; double m_time_step; double m_noise_pref_wo_gamma; public: - ParticleCoupling(LB::Solver &lb, bool couple_virtual, double time_step, - double kT) - : m_lb{lb}, m_couple_virtual{couple_virtual}, m_thermalized{kT != 0.}, - m_time_step{time_step} { + ParticleCoupling(LB::Solver &lb, double time_step, double kT) + : m_lb{lb}, m_thermalized{kT != 0.}, m_time_step{time_step} { assert(kT >= 0.); /* Eq. (16) @cite ahlrichs99a, without the gamma term. * The factor 12 comes from the fact that we use random numbers @@ -154,8 +151,8 @@ class ParticleCoupling { m_noise_pref_wo_gamma = std::sqrt(variance_inv * 2. * kT / time_step); } - ParticleCoupling(LB::Solver &lb, bool couple_virtual, double time_step) - : ParticleCoupling(lb, couple_virtual, time_step, + ParticleCoupling(LB::Solver &lb, double time_step) + : ParticleCoupling(lb, time_step, lb.get_kT() * Utils::sqr(lb.get_lattice_speed())) {} Utils::Vector3d get_noise_term(Particle const &p) const; @@ -188,6 +185,14 @@ class CouplingBookkeeping { public: /** @brief Determine if a given particle should be coupled. */ bool should_be_coupled(Particle const &p) { + auto const propagation = p.propagation(); + if (propagation == PropagationMode::NONE) { + return false; + } + if (propagation & (PropagationMode::TRANS_VS_RELATIVE | + PropagationMode::ROT_VS_RELATIVE)) { + return false; + } // real particles: always couple if (not p.is_ghost()) { return true; @@ -202,4 +207,8 @@ class CouplingBookkeeping { } }; +inline bool is_tracer(Particle const &p) { + return (p.propagation() & PropagationMode::TRANS_LB_TRACER) != 0; +} + } // namespace LB diff --git a/src/core/thermostat.cpp b/src/core/thermostat.cpp index 3b0ce275bad..cf8eab195c6 100644 --- a/src/core/thermostat.cpp +++ b/src/core/thermostat.cpp @@ -41,7 +41,6 @@ int thermo_switch = THERMO_OFF; double temperature = 0.0; -bool thermo_virtual = true; using Thermostat::GammaType; @@ -185,16 +184,6 @@ void mpi_set_langevin_gamma_rot(GammaType const &gamma) { mpi_call_all(mpi_set_langevin_gamma_rot_local, gamma); } -void mpi_set_thermo_virtual_local(bool thermo_virtual) { - ::thermo_virtual = thermo_virtual; -} - -REGISTER_CALLBACK(mpi_set_thermo_virtual_local) - -void mpi_set_thermo_virtual(bool thermo_virtual) { - mpi_call_all(mpi_set_thermo_virtual_local, thermo_virtual); -} - void mpi_set_temperature_local(double temperature) { ::temperature = temperature; try { diff --git a/src/core/thermostat.hpp b/src/core/thermostat.hpp index dfd0c507c36..88878664599 100644 --- a/src/core/thermostat.hpp +++ b/src/core/thermostat.hpp @@ -18,8 +18,9 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef CORE_THERMOSTAT_HPP -#define CORE_THERMOSTAT_HPP + +#pragma once + /** \file * Implementation in \ref thermostat.cpp. */ @@ -93,9 +94,6 @@ extern int thermo_switch; /** Temperature of the thermostat. */ extern double temperature; -/** True if the thermostat should act on virtual particles. */ -extern bool thermo_virtual; - #ifdef THERMOSTAT_PER_PARTICLE inline auto const & handle_particle_gamma(Thermostat::GammaType const &particle_gamma, @@ -409,8 +407,6 @@ void mpi_set_brownian_gamma_rot(Thermostat::GammaType const &gamma); void mpi_set_langevin_gamma(Thermostat::GammaType const &gamma); void mpi_set_langevin_gamma_rot(Thermostat::GammaType const &gamma); -void mpi_set_thermo_virtual(bool thermo_virtual); - void mpi_set_temperature(double temperature); void mpi_set_temperature_local(double temperature); @@ -421,5 +417,3 @@ void mpi_set_thermo_switch_local(int thermo_switch); void mpi_set_nptiso_gammas(double gamma0, double gammav); void mpi_set_nptiso_gammas_local(double gamma0, double gammav); #endif - -#endif diff --git a/src/core/thermostats/brownian_inline.hpp b/src/core/thermostats/brownian_inline.hpp index bc2ee7675f1..f010b03eeaa 100644 --- a/src/core/thermostats/brownian_inline.hpp +++ b/src/core/thermostats/brownian_inline.hpp @@ -19,8 +19,7 @@ * along with this program. If not, see . */ -#ifndef THERMOSTATS_BROWNIAN_INLINE_HPP -#define THERMOSTATS_BROWNIAN_INLINE_HPP +#pragma once #include "config/config.hpp" @@ -150,10 +149,6 @@ inline Utils::Vector3d bd_drag_vel(Thermostat::GammaType const &brownian_gamma, */ inline Utils::Vector3d bd_random_walk(BrownianThermostat const &brownian, Particle const &p, double dt, double kT) { - // skip the translation thermalizing for virtual sites unless enabled - if (p.is_virtual() and !thermo_virtual) - return {}; - Thermostat::GammaType sigma_pos = brownian.sigma_pos; #ifdef THERMOSTAT_PER_PARTICLE // override default if particle-specific gamma @@ -219,10 +214,6 @@ inline Utils::Vector3d bd_random_walk(BrownianThermostat const &brownian, */ inline Utils::Vector3d bd_random_walk_vel(BrownianThermostat const &brownian, Particle const &p) { - // skip the translation thermalizing for virtual sites unless enabled - if (p.is_virtual() and !thermo_virtual) - return {}; - auto const noise = Random::noise_gaussian( brownian.rng_counter(), brownian.rng_seed(), p.id()); Utils::Vector3d velocity = {}; @@ -384,5 +375,3 @@ bd_random_walk_vel_rot(BrownianThermostat const &brownian, Particle const &p) { return mask(p.rotation(), domega); } #endif // ROTATION - -#endif // THERMOSTATS_BROWNIAN_INLINE_HPP diff --git a/src/core/thermostats/langevin_inline.hpp b/src/core/thermostats/langevin_inline.hpp index 419d1314536..ec64500b51f 100644 --- a/src/core/thermostats/langevin_inline.hpp +++ b/src/core/thermostats/langevin_inline.hpp @@ -19,8 +19,7 @@ * along with this program. If not, see . */ -#ifndef THERMOSTATS_LANGEVIN_INLINE_HPP -#define THERMOSTATS_LANGEVIN_INLINE_HPP +#pragma once #include "config/config.hpp" @@ -44,11 +43,6 @@ inline Utils::Vector3d friction_thermo_langevin(LangevinThermostat const &langevin, Particle const &p, double time_step, double kT) { - // Early exit for virtual particles without thermostat - if (p.is_virtual() and !thermo_virtual) { - return {}; - } - // Determine prefactors for the friction and the noise term #ifdef THERMOSTAT_PER_PARTICLE auto const gamma = handle_particle_gamma(p.gamma(), langevin.gamma); @@ -99,4 +93,3 @@ friction_thermo_langevin_rotation(LangevinThermostat const &langevin, } #endif // ROTATION -#endif diff --git a/src/core/unit_tests/lb_particle_coupling_test.cpp b/src/core/unit_tests/lb_particle_coupling_test.cpp index 73ba0a71db0..50e512ceb55 100644 --- a/src/core/unit_tests/lb_particle_coupling_test.cpp +++ b/src/core/unit_tests/lb_particle_coupling_test.cpp @@ -183,7 +183,7 @@ BOOST_AUTO_TEST_CASE(rng) { lb_particle_coupling.gamma = 0.2; auto &lb = espresso::system->lb; - LB::ParticleCoupling coupling{lb, true, params.time_step, 1.}; + LB::ParticleCoupling coupling{lb, params.time_step, 1.}; BOOST_REQUIRE(lb_particle_coupling.rng_counter_coupling); BOOST_CHECK_EQUAL(lb_lbcoupling_get_rng_state(), 17); BOOST_CHECK(not lb_lbcoupling_is_seed_required()); @@ -209,7 +209,7 @@ BOOST_AUTO_TEST_CASE(rng) { BOOST_CHECK(step1_random1 != step2_random1); BOOST_CHECK(step1_random1 != step2_random2); - LB::ParticleCoupling coupling_unthermalized{lb, true, params.time_step, 0.}; + LB::ParticleCoupling coupling_unthermalized{lb, params.time_step, 0.}; auto const step3_norandom = coupling_unthermalized.get_noise_term(test_partcl_2); BOOST_CHECK((step3_norandom == Utils::Vector3d{0., 0., 0.})); @@ -218,7 +218,7 @@ BOOST_AUTO_TEST_CASE(rng) { BOOST_AUTO_TEST_CASE(drift_vel_offset) { Particle p{}; auto &lb = espresso::system->lb; - LB::ParticleCoupling coupling{lb, false, params.time_step}; + LB::ParticleCoupling coupling{lb, params.time_step}; BOOST_CHECK_EQUAL(coupling.lb_drift_velocity_offset(p).norm(), 0); Utils::Vector3d expected{}; #ifdef LB_ELECTROHYDRODYNAMICS @@ -261,7 +261,7 @@ BOOST_DATA_TEST_CASE(swimmer_force, bdata::make(kTs), kT) { // swimmer coupling { if (in_local_halo(p.pos())) { - LB::ParticleCoupling coupling{lb, true, params.time_step}; + LB::ParticleCoupling coupling{lb, params.time_step}; coupling.kernel(p); auto const interpolated = LB::get_force_to_be_applied(p.pos()); auto const expected = @@ -312,7 +312,7 @@ BOOST_DATA_TEST_CASE(particle_coupling, bdata::make(kTs), kT) { auto const gamma = 0.2; lb_lbcoupling_set_gamma(gamma); Particle p{}; - LB::ParticleCoupling coupling{lb, false, params.time_step}; + LB::ParticleCoupling coupling{lb, params.time_step}; auto expected = coupling.get_noise_term(p); #ifdef LB_ELECTROHYDRODYNAMICS p.mu_E() = Utils::Vector3d{-2., 1.5, 1.}; @@ -371,7 +371,7 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, set_particle_property(pid, &Particle::mu_E, Utils::Vector3d{-2., 1.5, 1.}); #endif - LB::ParticleCoupling coupling{lb, thermo_virtual, params.time_step}; + LB::ParticleCoupling coupling{lb, params.time_step}; auto const p_opt = copy_particle_to_head_node(comm, system, pid); auto expected = Utils::Vector3d{}; if (rank == 0) { @@ -443,8 +443,7 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, lb_lbcoupling_broadcast(); auto const particles = cell_structure.local_particles(); auto const ghost_particles = cell_structure.ghost_particles(); - LB::couple_particles(thermo_virtual, particles, ghost_particles, - params.time_step); + LB::couple_particles(particles, ghost_particles, params.time_step); auto const p_opt = copy_particle_to_head_node(comm, system, pid); if (rank == 0) { auto const &p = *p_opt; @@ -468,8 +467,7 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, } } // couple particle to LB - LB::couple_particles(thermo_virtual, particles, ghost_particles, - params.time_step); + LB::couple_particles(particles, ghost_particles, params.time_step); { auto const p_opt = copy_particle_to_head_node(comm, system, pid); if (rank == 0) { diff --git a/src/core/virtual_sites.hpp b/src/core/virtual_sites.hpp index 825316db5e2..f8ef445d510 100644 --- a/src/core/virtual_sites.hpp +++ b/src/core/virtual_sites.hpp @@ -50,7 +50,7 @@ inline void vs_relate_to(Particle &p_vs, Particle const &p_relate_to, // Set the particle id of the particle we want to relate to, the distance // and the relative orientation auto &vs_relative = p_vs.vs_relative(); - p_vs.propagation() = p_vs.propagation() | PropagationMode::TRANS_VS_RELATIVE; + p_vs.propagation() = PropagationMode::TRANS_VS_RELATIVE; vs_relative.to_particle_id = p_relate_to.id(); std::tie(vs_relative.rel_orientation, vs_relative.distance) = calculate_vs_relate_to_params(p_vs, p_relate_to, box_geo, min_global_cut); diff --git a/src/core/virtual_sites/lb_tracers.cpp b/src/core/virtual_sites/lb_tracers.cpp index 347fa396f97..936a6786e83 100644 --- a/src/core/virtual_sites/lb_tracers.cpp +++ b/src/core/virtual_sites/lb_tracers.cpp @@ -45,10 +45,6 @@ static bool lb_active_check(DeferredActiveLBChecks const &check) { return true; } -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_structure, double time_step) { DeferredActiveLBChecks check_lb_solver_set{}; @@ -68,7 +64,7 @@ void lb_tracers_add_particle_force_to_fluid(CellStructure &cell_structure, for (auto const &particle_range : {cell_structure.local_particles(), cell_structure.ghost_particles()}) { for (auto const &p : particle_range) { - if (!is_lb_tracer(p)) + if (!LB::is_tracer(p)) continue; if (!lb_active_check(check_lb_solver_set)) { return; @@ -94,7 +90,7 @@ void lb_tracers_propagate(CellStructure &cell_structure, double time_step) { // Advect particles for (auto &p : cell_structure.local_particles()) { - if (!is_lb_tracer(p)) + if (!LB::is_tracer(p)) continue; if (!lb_active_check(check_lb_solver_set)) { return; diff --git a/src/python/espressomd/particle_data.py b/src/python/espressomd/particle_data.py index 1161ead27ae..c49ec0e8201 100644 --- a/src/python/espressomd/particle_data.py +++ b/src/python/espressomd/particle_data.py @@ -416,12 +416,13 @@ def to_dict(self): """ pdict = self.get_params() - for k in ["director", "dip", "pos_folded", + for k in ["director", "dip", "pos_folded", "is_virtual", "image_box", "node", "lees_edwards_flag"]: if k in pdict: del pdict[k] if has_features("EXCLUSIONS"): pdict["exclusions"] = self.exclusions + pdict["propagation"] = self.propagation pdict["bonds"] = self.bonds return pdict diff --git a/src/python/espressomd/thermostat.pxd b/src/python/espressomd/thermostat.pxd index ccaceaacfc5..5830f75948a 100644 --- a/src/python/espressomd/thermostat.pxd +++ b/src/python/espressomd/thermostat.pxd @@ -25,7 +25,6 @@ from .utils cimport Vector3d cdef extern from "thermostat.hpp": double temperature int thermo_switch - cbool thermo_virtual int THERMO_OFF int THERMO_LANGEVIN int THERMO_LB @@ -103,7 +102,6 @@ cdef extern from "thermostat.hpp": void mpi_set_langevin_gamma(const double & gamma) void mpi_set_langevin_gamma_rot(const double & gamma) - void mpi_set_thermo_virtual(cbool thermo_virtual) void mpi_set_temperature(double temperature) void mpi_set_thermo_switch(int thermo_switch) diff --git a/src/python/espressomd/thermostat.pyx b/src/python/espressomd/thermostat.pyx index 80baef72f8b..2b912980fff 100644 --- a/src/python/espressomd/thermostat.pyx +++ b/src/python/espressomd/thermostat.pyx @@ -32,7 +32,7 @@ def AssertThermostatType(*allowed_thermostats): cdef class Thermostat: @AssertThermostatType(THERMO_LANGEVIN, THERMO_DPD) def set_langevin(self, kT=None, gamma=None, gamma_rotation=None, - act_on_virtual=False, seed=None): + seed=None): ... This will prefix an assertion that prevents setting up the Langevin @@ -104,7 +104,6 @@ cdef class Thermostat: if thmst["type"] == "LANGEVIN": self.set_langevin(kT=thmst["kT"], gamma=thmst["gamma"], gamma_rotation=thmst["gamma_rotation"], - act_on_virtual=thmst["act_on_virtual"], seed=thmst["seed"]) langevin_set_rng_counter(thmst["counter"]) if thmst["type"] == "LB": @@ -117,7 +116,6 @@ cdef class Thermostat: thmst["LB_fluid"].call_method("activate") self.set_lb( LB_fluid=thmst["LB_fluid"], - act_on_virtual=thmst["act_on_virtual"], gamma=thmst["gamma"], seed=thmst["rng_counter_fluid"]) thmst["LB_fluid"].call_method("deactivate") @@ -133,7 +131,6 @@ cdef class Thermostat: if thmst["type"] == "BROWNIAN": self.set_brownian(kT=thmst["kT"], gamma=thmst["gamma"], gamma_rotation=thmst["gamma_rotation"], - act_on_virtual=thmst["act_on_virtual"], seed=thmst["seed"]) brownian_set_rng_counter(thmst["counter"]) if thmst["type"] == "SD": @@ -155,7 +152,6 @@ cdef class Thermostat: lang_dict = {} lang_dict["type"] = "LANGEVIN" lang_dict["kT"] = temperature - lang_dict["act_on_virtual"] = thermo_virtual lang_dict["seed"] = langevin.rng_seed() lang_dict["counter"] = langevin.rng_counter() IF PARTICLE_ANISOTROPY: @@ -179,7 +175,6 @@ cdef class Thermostat: lang_dict = {} lang_dict["type"] = "BROWNIAN" lang_dict["kT"] = temperature - lang_dict["act_on_virtual"] = thermo_virtual lang_dict["seed"] = brownian.rng_seed() lang_dict["counter"] = brownian.rng_counter() IF PARTICLE_ANISOTROPY: @@ -204,7 +199,6 @@ cdef class Thermostat: lb_dict["LB_fluid"] = self._LB_fluid lb_dict["gamma"] = lb_lbcoupling_get_gamma() lb_dict["type"] = "LB" - lb_dict["act_on_virtual"] = thermo_virtual lb_dict["rng_counter_fluid"] = lb_lbcoupling_get_rng_state() thermo_list.append(lb_dict) if thermo_switch & THERMO_NPT_ISO: @@ -246,7 +240,6 @@ cdef class Thermostat: """ self._set_temperature(0.) - mpi_set_thermo_virtual(True) IF PARTICLE_ANISOTROPY: mpi_set_langevin_gamma(utils.make_Vector3d((0., 0., 0.))) mpi_set_brownian_gamma(utils.make_Vector3d((0., 0., 0.))) @@ -265,8 +258,7 @@ cdef class Thermostat: mpi_bcast_lb_particle_coupling() @AssertThermostatType(THERMO_LANGEVIN, THERMO_DPD) - def set_langevin(self, kT, gamma, gamma_rotation=None, - act_on_virtual=False, seed=None): + def set_langevin(self, kT, gamma, gamma_rotation=None, seed=None): """ Sets the Langevin thermostat. @@ -283,9 +275,6 @@ cdef class Thermostat: The same applies to ``gamma_rotation``, which requires the feature ``ROTATION`` to work properly. But also accepts three floats if ``PARTICLE_ANISOTROPY`` is also compiled in. - act_on_virtual : :obj:`bool`, optional - If ``True`` the thermostat will act on virtual sites, default is - ``False``. seed : :obj:`int` Initial counter value (or seed) of the philox RNG. Required on first activation of the Langevin thermostat. @@ -390,11 +379,8 @@ cdef class Thermostat: IF ROTATION: mpi_set_langevin_gamma_rot(gamma_rotation) - mpi_set_thermo_virtual(act_on_virtual) - @AssertThermostatType(THERMO_BROWNIAN) - def set_brownian(self, kT, gamma, gamma_rotation=None, - act_on_virtual=False, seed=None): + def set_brownian(self, kT, gamma, gamma_rotation=None, seed=None): """Sets the Brownian thermostat. Parameters @@ -410,9 +396,6 @@ cdef class Thermostat: The same applies to ``gamma_rotation``, which requires the feature ``ROTATION`` to work properly. But also accepts three floats if ``PARTICLE_ANISOTROPY`` is also compiled in. - act_on_virtual : :obj:`bool`, optional - If ``True`` the thermostat will act on virtual sites, default is - ``False``. seed : :obj:`int` Initial counter value (or seed) of the philox RNG. Required on first activation of the Brownian thermostat. @@ -518,15 +501,8 @@ cdef class Thermostat: IF ROTATION: mpi_set_brownian_gamma_rot(gamma_rotation) - mpi_set_thermo_virtual(act_on_virtual) - @AssertThermostatType(THERMO_LB, THERMO_DPD) - def set_lb( - self, - seed=None, - act_on_virtual=True, - LB_fluid=None, - gamma=0.0): + def set_lb(self, seed=None, LB_fluid=None, gamma=0.0): """ Sets the LB thermostat. @@ -538,8 +514,6 @@ cdef class Thermostat: seed : :obj:`int` Seed for the random number generator, required if kT > 0. Must be positive. - act_on_virtual : :obj:`bool`, optional - If ``True`` the thermostat will act on virtual sites (default). gamma : :obj:`float` Frictional coupling constant for the MD particle coupling. @@ -568,7 +542,6 @@ cdef class Thermostat: global thermo_switch mpi_set_thermo_switch(thermo_switch | THERMO_LB) - mpi_set_thermo_virtual(act_on_virtual) lb_lbcoupling_set_gamma(gamma) mpi_bcast_lb_particle_coupling() diff --git a/src/script_interface/particle_data/ParticleHandle.cpp b/src/script_interface/particle_data/ParticleHandle.cpp index 331fe5258aa..edab23a5d67 100644 --- a/src/script_interface/particle_data/ParticleHandle.cpp +++ b/src/script_interface/particle_data/ParticleHandle.cpp @@ -34,6 +34,7 @@ #include "core/rotation.hpp" #include "core/system/System.hpp" #include "core/virtual_sites.hpp" + #include #include @@ -249,18 +250,7 @@ ParticleHandle::ParticleHandle() { }, #endif // ELECTROSTATICS [this]() { return get_particle_data(m_pid).q(); }}, - {"virtual", -#ifdef VIRTUAL_SITES - [this](Variant const &value) { - set_particle_property(&Particle::virtual_flag, value); - }, -#else // VIRTUAL_SITES - [](Variant const &value) { - if (get_value(value)) { - throw std::runtime_error("Feature VIRTUAL_SITES not compiled in"); - } - }, -#endif // VIRTUAL_SITES + {"is_virtual", AutoParameter::read_only, [this]() { return get_particle_data(m_pid).is_virtual(); }}, #ifdef ROTATION {"director", @@ -459,6 +449,11 @@ ParticleHandle::ParticleHandle() { } set_particle_property( [&vs_relative](Particle &p) { p.vs_relative() = vs_relative; }); + if (vs_relative.to_particle_id != -1) { + set_particle_property([](Particle &p) { + p.propagation() = PropagationMode::TRANS_VS_RELATIVE; + }); + } }, [this]() { auto const vs_rel = get_particle_data(m_pid).vs_relative(); @@ -628,10 +623,8 @@ Variant ParticleHandle::do_call_method(std::string const &name, auto const [quat, dist] = calculate_vs_relate_to_params( p_current, p_relate_to, *system.box_geo, system.get_min_global_cut(), override_cutoff_check); - set_parameter("propagation", PropagationMode::TRANS_VS_RELATIVE); set_parameter("vs_relative", Variant{std::vector{ {other_pid, dist, quat2vector(quat)}}}); - set_parameter("virtual", true); #endif // VIRTUAL_SITES_RELATIVE #ifdef EXCLUSIONS } else if (name == "has_exclusion") { @@ -781,8 +774,9 @@ void ParticleHandle::do_construct(VariantMap const ¶ms) { context()->parallel_try_catch([&]() { // set particle properties (filter out read-only and deferred properties) std::set const skip = { - "pos_folded", "pos", "quat", "director", "id", "lees_edwards_flag", - "exclusions", "dip", "node", "image_box", "bonds", "__cpt_sentinel", + "lees_edwards_flag", "pos_folded", "pos", "quat", "director", "id", + "is_virtual", "exclusions", "dip", "node", "image_box", "bonds", + "__cpt_sentinel", }; #ifdef ROTATION // multiple parameters can potentially set the quaternion, but only one diff --git a/testsuite/python/analyze_chains.py b/testsuite/python/analyze_chains.py index 0fe3a7a34a0..0669485e81b 100644 --- a/testsuite/python/analyze_chains.py +++ b/testsuite/python/analyze_chains.py @@ -212,9 +212,10 @@ def test_exceptions(self): with self.assertRaisesRegex(ValueError, "needs at least 1 chain"): method(chain_start=0, number_of_chains=-1, chain_length=1) self.assertIsNone(analysis.call_method("unknown")) - if espressomd.has_features("VIRTUAL_SITES"): + if espressomd.has_features("VIRTUAL_SITES_RELATIVE"): with self.assertRaisesRegex(RuntimeError, "Center of mass is not well-defined"): - self.system.part.by_id(0).virtual = True + p = self.system.part.by_id(0) + p.vs_relative = (1, 0.01, (1., 0., 0., 0.)) analysis.calc_rg(chain_start=0, number_of_chains=num_poly, chain_length=num_mono) with self.assertRaisesRegex(RuntimeError, "Parameter 'analysis' is read-only"): diff --git a/testsuite/python/analyze_mass_related.py b/testsuite/python/analyze_mass_related.py index 6b114ad9332..d527667cabf 100644 --- a/testsuite/python/analyze_mass_related.py +++ b/testsuite/python/analyze_mass_related.py @@ -37,6 +37,7 @@ class AnalyzeMassRelated(ut.TestCase): @classmethod def setUpClass(cls): + from espressomd.propagation import Propagation cls.system.cell_system.skin = 0.4 cls.system.time_step = 0.01 cls.system.thermostat.turn_off() @@ -46,7 +47,8 @@ def setUpClass(cls): pos=np.random.random((10, 3)) * cls.box_l, type=[1] * 10) if espressomd.has_features("VIRTUAL_SITES"): cls.system.part.add( - pos=np.random.random((10, 3)) * cls.box_l, type=[0] * 10, virtual=[True] * 10) + pos=np.random.random((10, 3)) * cls.box_l, type=[0] * 10, + propagation=[Propagation.TRANS_LB_TRACER] * 10) all_partcls = cls.system.part.all() all_partcls.v = np.random.random( (len(all_partcls), 3)) + 0.1 @@ -76,7 +78,7 @@ def i_tensor(self, ids): def test_itensor(self): # Particles of type 0 I0 = self.i_tensor(self.system.part.select( - lambda p: (not p.virtual) and p.type == 0).id) + lambda p: (not p.is_virtual) and p.type == 0).id) np.testing.assert_allclose( I0, self.system.analysis.moment_of_inertia_matrix(p_type=0), atol=1E-9) # Particles of type 1 @@ -86,7 +88,7 @@ def test_itensor(self): def test_center_of_mass(self): no_virtual_type_0 = self.system.part.select( - lambda p: (not p.virtual) and p.type == 0) + lambda p: (not p.is_virtual) and p.type == 0) com = np.zeros(3) for p in no_virtual_type_0: com += p.pos * p.mass @@ -97,7 +99,7 @@ def test_center_of_mass(self): self.system.analysis.center_of_mass(p_type=0)) def test_galilei_transform(self): - no_virtual = self.system.part.select(virtual=False) + no_virtual = self.system.part.select(is_virtual=False) # Center of mass np.testing.assert_allclose( @@ -110,7 +112,7 @@ def test_galilei_transform(self): def test_angularmomentum(self): no_virtual_type_0 = self.system.part.select( - lambda p: (not p.virtual) and p.type == 0) + lambda p: (not p.is_virtual) and p.type == 0) am = np.zeros(3) for p in no_virtual_type_0: am += p.mass * np.cross(p.pos, p.v) @@ -120,14 +122,14 @@ def test_angularmomentum(self): self.system.analysis.angular_momentum(p_type=0)) def test_kinetic_energy(self): - no_virtual = self.system.part.select(virtual=False) + no_virtual = self.system.part.select(is_virtual=False) E_kin = 0.5 * \ np.sum(no_virtual.mass * np.sum(no_virtual.v**2, axis=1), axis=0) np.testing.assert_allclose( E_kin, self.system.analysis.energy()["kinetic"]) def test_kinetic_pressure(self): - no_virtual = self.system.part.select(virtual=False) + no_virtual = self.system.part.select(is_virtual=False) P_kin = np.sum( no_virtual.mass * np.sum(no_virtual.v**2, axis=1), axis=0) / (3 * self.system.volume()) @@ -142,7 +144,7 @@ def test_kinetic_pressure(self): expected_pressure_tensor, analyze_pressure_tensor) def test_gyration_radius(self): - if self.system.part.select(virtual=True): + if self.system.part.select(is_virtual=True): with self.assertRaisesRegex(RuntimeError, "not well-defined"): self.system.analysis.calc_rg(chain_start=0, number_of_chains=1, chain_length=len(self.system.part)) diff --git a/testsuite/python/brownian_dynamics.py b/testsuite/python/brownian_dynamics.py index 6c87186eff5..acc6e2d9837 100644 --- a/testsuite/python/brownian_dynamics.py +++ b/testsuite/python/brownian_dynamics.py @@ -127,16 +127,17 @@ def test_01__rng_per_particle(self): """Test for RNG consistency.""" self.check_rng(True) - @utx.skipIfMissingFeatures("VIRTUAL_SITES") def test_07__virtual(self): + from espressomd.propagation import Propagation system = self.system system.time_step = 0.01 - virtual = system.part.add(pos=[0, 0, 0], virtual=True, v=[1, 0, 0]) - physical = system.part.add(pos=[0, 0, 0], virtual=False, v=[1, 0, 0]) + virtual = system.part.add(pos=[0, 0, 0], v=[1, 0, 0], + propagation=Propagation.NONE) + physical = system.part.add(pos=[0, 0, 0], v=[1, 0, 0]) system.thermostat.set_brownian( - kT=0, gamma=1, gamma_rotation=1., act_on_virtual=False, seed=41) + kT=0, gamma=1, gamma_rotation=1., seed=41) system.integrator.set_brownian_dynamics() system.integrator.run(1) @@ -144,14 +145,6 @@ def test_07__virtual(self): np.testing.assert_almost_equal(np.copy(virtual.v), [1, 0, 0]) np.testing.assert_almost_equal(np.copy(physical.v), [0, 0, 0]) - virtual.pos = physical.pos = [0, 0, 0] - system.thermostat.set_brownian( - kT=0, gamma=1, gamma_rotation=1., act_on_virtual=True, seed=41) - system.integrator.run(1) - - np.testing.assert_almost_equal(np.copy(virtual.v), [0, 0, 0]) - np.testing.assert_almost_equal(np.copy(physical.v), [0, 0, 0]) - @utx.skipIfMissingFeatures(["ROTATION", "EXTERNAL_FORCES"]) def test_fix_rotation(self): system = self.system @@ -164,7 +157,7 @@ def test_fix_rotation(self): # torque only part.ext_torque = [0, 0, 1.3] system.thermostat.set_brownian( - kT=0, gamma=1, gamma_rotation=1.5, act_on_virtual=False, seed=41) + kT=0, gamma=1, gamma_rotation=1.5, seed=41) system.integrator.set_brownian_dynamics() system.integrator.run(3) np.testing.assert_allclose( @@ -174,7 +167,7 @@ def test_fix_rotation(self): # noise only part.ext_torque = 3 * [0.] system.thermostat.set_brownian( - kT=1, gamma=1, gamma_rotation=1.5, act_on_virtual=False, seed=41) + kT=1, gamma=1, gamma_rotation=1.5, seed=41) system.integrator.run(3) self.assertGreater(np.linalg.norm(part.omega_lab), 0.) diff --git a/testsuite/python/collision_detection.py b/testsuite/python/collision_detection.py index 141e65a0af7..8c9f9c744f5 100644 --- a/testsuite/python/collision_detection.py +++ b/testsuite/python/collision_detection.py @@ -160,7 +160,7 @@ def verify_state_after_bind_at_poc_pair(self, expected_np): # get partner p2 = system.part.by_id(p.bonds[0][1]) # Is that really a vs - self.assertTrue(p2.virtual) + self.assertTrue(p2.is_virtual) # Get base particles base_p1 = system.part.by_id(p.vs_relative[0]) base_p2 = system.part.by_id(p2.vs_relative[0]) @@ -187,7 +187,7 @@ def verify_state_after_bind_at_poc_triplet(self, expected_np): # From the two vs we find the two non-virtual particles for p in system.part: # Skip non-virtual - if not p.virtual: + if not p.is_virtual: continue # Skip vs that doesn't have a bond if p.bonds == (): @@ -204,7 +204,7 @@ def verify_state_after_bind_at_poc_triplet(self, expected_np): # Bond type self.assertEqual(p.bonds[0][0], self.bond_angle_vs) # Is that really a vs - self.assertTrue(p.virtual) + self.assertTrue(p.is_virtual) # Get base particles base_p1 = system.part.by_id(p.bonds[0][1]) base_p2 = system.part.by_id(p.bonds[0][2]) @@ -217,7 +217,7 @@ def verify_state_after_bind_at_poc_triplet(self, expected_np): # Check particle that did not take part in collision. self.assertEqual(len(parts_not_accounted_for), 1) p = system.part.by_id(parts_not_accounted_for[0]) - self.assertFalse(p.virtual) + self.assertFalse(p.is_virtual) self.assertEqual(p.bonds, ()) parts_not_accounted_for.remove(p.id) self.assertEqual(parts_not_accounted_for, []) @@ -243,8 +243,8 @@ def verify_bind_at_poc(self, p1, p2, vs1, vs2): self.assertTrue(vs1.bonds == bond_vs1 or vs2.bonds == bond_vs2) # Vs properties - self.assertTrue(vs1.virtual) - self.assertTrue(vs2.virtual) + self.assertTrue(vs1.is_virtual) + self.assertTrue(vs2.is_virtual) # vs_relative properties seen = [] @@ -368,8 +368,8 @@ def test_bind_at_point_of_collision_random(self): system.integrator.run(5000) # Analysis - virtual_sites = system.part.select(virtual=True) - non_virtual = system.part.select(virtual=False) + virtual_sites = system.part.select(is_virtual=True) + non_virtual = system.part.select(is_virtual=False) # Check bonds on non-virtual particles bonds = [] diff --git a/testsuite/python/field_test.py b/testsuite/python/field_test.py index 051fe326a58..c13b829f939 100644 --- a/testsuite/python/field_test.py +++ b/testsuite/python/field_test.py @@ -46,6 +46,7 @@ def tearDown(self): self.system.part.clear() def test_gravity(self): + from espressomd.propagation import Propagation g_const = np.array([1, 2, 3]) gravity = espressomd.constraints.Gravity(g=g_const) @@ -64,8 +65,8 @@ def test_gravity(self): self.assertAlmostEqual(self.system.analysis.energy()['total'], 0.) # Virtual sites don't feel gravity - if espressomd.has_features("VIRTUAL_SITES"): - p.virtual = True + if espressomd.has_features("VIRTUAL_SITES_INERTIALESS_TRACERS"): + p.propagation = Propagation.TRANS_LB_TRACER self.system.integrator.run(0) np.testing.assert_allclose(np.copy(p.f), 0) diff --git a/testsuite/python/langevin_thermostat.py b/testsuite/python/langevin_thermostat.py index 1d701e0eef6..86c9fa50abb 100644 --- a/testsuite/python/langevin_thermostat.py +++ b/testsuite/python/langevin_thermostat.py @@ -206,29 +206,23 @@ def test_03__friction_rot(self): np.testing.assert_allclose( np.copy(p.omega_body), ref_omega_body, atol=5E-4) - @utx.skipIfMissingFeatures("VIRTUAL_SITES") + @utx.skipIfMissingFeatures("VIRTUAL_SITES_RELATIVE") def test_07__virtual(self): system = self.system system.time_step = 0.01 - virtual = system.part.add(pos=[0, 0, 0], virtual=True, v=[1, 0, 0]) - physical = system.part.add(pos=[0, 0, 0], virtual=False, v=[1, 0, 0]) + virtual = system.part.add(pos=[0, 0, 0], v=[1, 0, 0]) + physical = system.part.add(pos=[0, 0, 0], v=[1, 0, 0]) + virtual.vs_relative = (physical.id, 0.01, (1., 0., 0., 0.)) system.thermostat.set_langevin( - kT=0, gamma=1, gamma_rotation=1., act_on_virtual=False, seed=41) + kT=0, gamma=1, gamma_rotation=1., seed=41) system.integrator.run(0) np.testing.assert_almost_equal(np.copy(virtual.f), [0, 0, 0]) np.testing.assert_almost_equal(np.copy(physical.f), [-1, 0, 0]) - system.thermostat.set_langevin( - kT=0, gamma=1, gamma_rotation=1., act_on_virtual=True, seed=41) - system.integrator.run(0) - - np.testing.assert_almost_equal(np.copy(virtual.f), [-1, 0, 0]) - np.testing.assert_almost_equal(np.copy(physical.f), [-1, 0, 0]) - if __name__ == "__main__": ut.main() diff --git a/testsuite/python/lb_thermo_virtual.py b/testsuite/python/lb_thermo_virtual.py index ab9875f4df4..73f3242231a 100644 --- a/testsuite/python/lb_thermo_virtual.py +++ b/testsuite/python/lb_thermo_virtual.py @@ -24,17 +24,10 @@ import numpy as np -@utx.skipIfMissingFeatures(["VIRTUAL_SITES"]) -class LBBoundaryThermoVirtualTest(ut.TestCase): - - """Test slip velocity of boundaries. - - In this simple test, a wall with slip velocity is - added and the fluid is checked if it has the same velocity. - """ +@utx.skipIfMissingFeatures(["WALBERLA", "VIRTUAL_SITES_INERTIALESS_TRACERS"]) +class Test(ut.TestCase): system = espressomd.System(box_l=[10.0, 10.0, 10.0]) - system.time_step = 1.0 system.cell_system.skin = 0.1 @@ -42,50 +35,24 @@ def tearDown(self): self.system.lb = None self.system.part.clear() - def check_virtual(self, fluid_class): + def test_virtual(self): + from espressomd.propagation import Propagation system = self.system - lb_fluid = fluid_class( + lb_fluid = espressomd.lb.LBFluidWalberla( agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=1.0) system.lb = lb_fluid - virtual = system.part.add(pos=[0, 0, 0], virtual=True, v=[1, 0, 0]) - physical = system.part.add(pos=[0, 0, 0], virtual=False, v=[1, 0, 0]) + virtual = system.part.add(pos=[0, 0, 0], v=[1, 0, 0], + propagation=Propagation.TRANS_LB_TRACER) + physical = system.part.add(pos=[0, 0, 0], v=[1, 0, 0]) - system.thermostat.set_lb( - LB_fluid=lb_fluid, - act_on_virtual=False, - gamma=1.0) + system.thermostat.set_lb(LB_fluid=lb_fluid, gamma=1.0) system.integrator.run(1) np.testing.assert_almost_equal(np.copy(virtual.f), [0, 0, 0]) np.testing.assert_almost_equal(np.copy(physical.f), [-1, 0, 0]) - system.thermostat.set_lb(LB_fluid=lb_fluid, act_on_virtual=True) - - virtual.v = [1, 0, 0] - physical.v = [1, 0, 0] - - system.lb = None - lb_fluid = fluid_class( - agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=1.0) - system.lb = lb_fluid - system.thermostat.set_lb(LB_fluid=lb_fluid, gamma=1.0) - virtual.pos = physical.pos - virtual.v = [1, 0, 0] - physical.v = [1, 0, 0] - system.integrator.run(1) - - # The forces are not exactly -1 because the fluid is not at - # rest anymore because of the previous check. - np.testing.assert_almost_equal(np.copy(physical.f), np.copy(virtual.f)) - np.testing.assert_almost_equal(np.copy(physical.f), [-1, 0, 0]) - np.testing.assert_almost_equal(np.copy(virtual.f), [-1, 0, 0]) - - @utx.skipIfMissingFeatures(["WALBERLA"]) - def test_lb_walberla(self): - self.check_virtual(espressomd.lb.LBFluidWalberla) - if __name__ == "__main__": ut.main() diff --git a/testsuite/python/lees_edwards.py b/testsuite/python/lees_edwards.py index 80a63de5b4b..17aecfe5483 100644 --- a/testsuite/python/lees_edwards.py +++ b/testsuite/python/lees_edwards.py @@ -616,7 +616,7 @@ def test_le_colldet(self): # generated VS. The other components are inherited from the real # particles. box_l = np.copy(system.box_l) - p_vs = system.part.select(virtual=True) + p_vs = system.part.select(is_virtual=True) np.testing.assert_array_almost_equal( np.minimum(np.abs(p_vs.pos[:, 0] - col_part1.pos[0]), np.abs(p_vs.pos[:, 0] - col_part2.pos[0])), 0.) diff --git a/testsuite/python/npt_thermostat.py b/testsuite/python/npt_thermostat.py index 70bdd088c1e..375dd16b096 100644 --- a/testsuite/python/npt_thermostat.py +++ b/testsuite/python/npt_thermostat.py @@ -140,12 +140,13 @@ def test_02__direction(self): np.testing.assert_allclose(box_l_rel, box_l_rel_ref, atol=1e-10) self.assertGreater(np.max(box_l_rel), 2) - @utx.skipIfMissingFeatures("VIRTUAL_SITES") def test_07__virtual(self): + from espressomd.propagation import Propagation system = self.system - virtual = system.part.add(pos=[0, 0, 0], virtual=True, v=[1, 0, 0]) - physical = system.part.add(pos=[0, 0, 0], virtual=False, v=[1, 0, 0]) + virtual = system.part.add(pos=[0, 0, 0], v=[1, 0, 0], + propagation=Propagation.NONE) + physical = system.part.add(pos=[0, 0, 0], v=[1, 0, 0]) system.thermostat.set_npt(kT=1.0, gamma0=2.0, gammav=0.04, seed=42) system.integrator.set_isotropic_npt(ext_pressure=2.0, piston=0.01) diff --git a/testsuite/python/observables.py b/testsuite/python/observables.py index eef1a9e8b6d..09974046a1b 100644 --- a/testsuite/python/observables.py +++ b/testsuite/python/observables.py @@ -20,6 +20,7 @@ import espressomd import numpy as np import espressomd.observables +from espressomd.propagation import Propagation def calc_com_x(system, x, id_list): @@ -28,7 +29,7 @@ def calc_com_x(system, x, id_list): masses = partcls.mass # Filter out virtual particles by using mass=0 for them - virtual = partcls.virtual + virtual = partcls.is_virtual masses[np.nonzero(virtual)] = 0. return np.average(getattr(partcls, x), weights=masses, axis=0) @@ -68,7 +69,7 @@ class Observables(ut.TestCase): if espressomd.has_features("VIRTUAL_SITES"): p = system.part.by_id(partcls.id[8]) - p.virtual = True + p.propagation = Propagation.TRANS_LB_TRACER def generate_test_for_pid_observable( _obs_class, _pprop_name, _agg_type=None): @@ -96,7 +97,7 @@ def func(self): # Get data from particles if pprop_name == "f": for p_id in id_list: - if self.system.part.by_id(p_id).virtual: + if self.system.part.by_id(p_id).is_virtual: id_list.remove(p_id) part_data = getattr(self.system.part.by_ids(id_list), pprop_name) @@ -234,7 +235,7 @@ def test_com_force(self): replace=False)) particles = self.system.part.select( - lambda p: p.id in id_list and not p.virtual) + lambda p: p.id in id_list and not p.is_virtual) np.testing.assert_allclose( np.sum(particles.f, axis=0), diff --git a/testsuite/python/particle.py b/testsuite/python/particle.py index 3ca8ba623ed..22d60b0e015 100644 --- a/testsuite/python/particle.py +++ b/testsuite/python/particle.py @@ -186,15 +186,14 @@ def test_gamma_rot_single(self): "dip", np.array([0.5, -0.5, 3])) test_dipm = generateTestForScalarProperty("dipm", -9.7) - if espressomd.has_features(["VIRTUAL_SITES"]): - test_virtual = generateTestForScalarProperty("virtual", True) - @utx.skipIfMissingFeatures(["VIRTUAL_SITES_RELATIVE"]) def test_vs_relative(self): self.system.part.add(id=0, pos=(0, 0, 0)) p1 = self.system.part.add(id=1, pos=(0, 0, 0)) + self.assertFalse(p1.is_virtual) p1.vs_relative = (0, 5.0, (0.5, -0.5, -0.5, -0.5)) p1.vs_quat = [1, 2, 3, 4] + self.assertTrue(p1.is_virtual) np.testing.assert_array_equal(p1.vs_quat, [1, 2, 3, 4]) res = p1.vs_relative self.assertEqual(res[0], 0, f"vs_relative: {res}") @@ -395,7 +394,7 @@ def check(feature, prop, throwing_values, valid_value=None): check("MASS", "mass", [1.1, 0., -1.], 1.) check("ELECTROSTATICS", "q", [1., -1.], 0.) - check("VIRTUAL_SITES", "virtual", [True], False) + check("VIRTUAL_SITES", "is_virtual", [True], False) def test_parallel_property_setters(self): system = self.system diff --git a/testsuite/python/rotate_system.py b/testsuite/python/rotate_system.py index 66a6903a9da..ded7eda4e18 100644 --- a/testsuite/python/rotate_system.py +++ b/testsuite/python/rotate_system.py @@ -70,12 +70,14 @@ def test_no_mass(self): # Check that virtual sites do not influence the center of mass # calculation - if espressomd.has_features("VIRTUAL_SITES"): - p2 = system.part.add(pos=p1.pos, virtual=True) + if espressomd.has_features("VIRTUAL_SITES_RELATIVE"): + vs_dist = 0.01 + p2 = system.part.add(pos=p1.pos) + p2.vs_relative = (p1.id, vs_dist, (1., 0., 0., 0.)) system.rotate_system(phi=pi / 2., theta=pi / 2., alpha=-pi / 2.) np.testing.assert_allclose(np.copy(p0.pos), [6, 4, 4]) np.testing.assert_allclose(np.copy(p1.pos), [4, 6, 6]) - np.testing.assert_allclose(np.copy(p2.pos), [4, 6, 6]) + np.testing.assert_allclose(np.copy(p2.pos), [4, 6, 6 + vs_dist]) if __name__ == "__main__": diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index 434655212e1..27759b9b368 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -511,7 +511,6 @@ def test_thermostat_Langevin(self): self.assertEqual(thmst['kT'], 1.0) self.assertEqual(thmst['seed'], 42) self.assertEqual(thmst['counter'], 0) - self.assertFalse(thmst['act_on_virtual']) np.testing.assert_array_equal(thmst['gamma'], 3 * [2.0]) if espressomd.has_features('ROTATION'): np.testing.assert_array_equal(thmst['gamma_rotation'], 3 * [2.0]) @@ -524,7 +523,6 @@ def test_thermostat_Brownian(self): self.assertEqual(thmst['kT'], 1.0) self.assertEqual(thmst['seed'], 42) self.assertEqual(thmst['counter'], 0) - self.assertFalse(thmst['act_on_virtual']) np.testing.assert_array_equal(thmst['gamma'], 3 * [2.0]) if espressomd.has_features('ROTATION'): np.testing.assert_array_equal(thmst['gamma_rotation'], 3 * [2.0]) diff --git a/testsuite/python/virtual_sites_relative.py b/testsuite/python/virtual_sites_relative.py index f16f31c69ef..f091a5562e6 100644 --- a/testsuite/python/virtual_sites_relative.py +++ b/testsuite/python/virtual_sites_relative.py @@ -185,7 +185,7 @@ def test_pos_vel_forces(self): propagation=Propagation.TRANS_VS_RELATIVE) p.vs_auto_relate_to(p1) # Was the particle made virtual - self.assertTrue(p.virtual) + self.assertTrue(p.is_virtual) # Are vs relative to id and vs_r = p.vs_relative # id diff --git a/testsuite/python/virtual_sites_tracers_common.py b/testsuite/python/virtual_sites_tracers_common.py index fb104d6ee72..622a8ecabc3 100644 --- a/testsuite/python/virtual_sites_tracers_common.py +++ b/testsuite/python/virtual_sites_tracers_common.py @@ -49,10 +49,7 @@ def set_lb(self, ext_force_density=(0, 0, 0), dir_walls=2): kT=0.0, agrid=1., density=1., kinematic_viscosity=1.8, tau=self.system.time_step, ext_force_density=ext_force_density) self.system.lb = self.lbf - self.system.thermostat.set_lb( - LB_fluid=self.lbf, - act_on_virtual=False, - gamma=1) + self.system.thermostat.set_lb(LB_fluid=self.lbf, gamma=1) # Setup boundaries normal = [0, 0, 0] @@ -155,9 +152,8 @@ def test_advection(self): system.lb = None def test_zz_exceptions_without_lb(self): - """Check behaviour without lb. Ignore non-virtual particles, complain on - virtual ones. - + """ + Check behaviour without LB. Ignore real particles, complain on tracers. """ self.set_lb() system = self.system