diff --git a/src/core/Particle.hpp b/src/core/Particle.hpp index d3c6a2e0e7..46b9085e04 100644 --- a/src/core/Particle.hpp +++ b/src/core/Particle.hpp @@ -513,11 +513,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 is_virtual() const { return (p.propagation & (PropagationMode::TRANS_VS_RELATIVE | PropagationMode::ROT_VS_RELATIVE | PropagationMode::TRANS_LB_TRACER)) != 0; } +#else + constexpr auto is_virtual() const { return false; } +#endif // VIRTUAL_SITES #ifdef VIRTUAL_SITES_RELATIVE auto const &vs_relative() const { return p.vs_relative; } auto &vs_relative() { return p.vs_relative; } diff --git a/src/core/PropagationMode.hpp b/src/core/PropagationMode.hpp index 8afe90ff23..ef8e029e81 100644 --- a/src/core/PropagationMode.hpp +++ b/src/core/PropagationMode.hpp @@ -23,7 +23,7 @@ namespace PropagationMode { /** * @brief Flags to create bitmasks for propagation modes. */ -enum PropagationMode { +enum PropagationMode : int { NONE = 0, SYSTEM_DEFAULT = 1 << 0, TRANS_NEWTON = 1 << 1, diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index c5f1ab4e9f..4f01e75fb8 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -326,22 +326,14 @@ static bool should_propagate_with(Particle const &p, int mode) { } void System::System::thermostats_force_init(double kT) { - int bitmask = PropagationMode::TRANS_LANGEVIN | PropagationMode::ROT_LANGEVIN; - if (::used_propagations & bitmask) { - if (default_propagation & bitmask) { - bitmask |= PropagationMode::SYSTEM_DEFAULT; - } - for (auto &p : cell_structure->local_particles()) { - if (p.propagation() & bitmask) { - if (should_propagate_with(p, PropagationMode::TRANS_LANGEVIN)) - p.force() += friction_thermo_langevin(langevin, p, time_step, kT); + for (auto &p : cell_structure->local_particles()) { + if (should_propagate_with(p, PropagationMode::TRANS_LANGEVIN)) + p.force() += friction_thermo_langevin(langevin, p, time_step, kT); #ifdef ROTATION - if (should_propagate_with(p, PropagationMode::ROT_LANGEVIN)) - p.torque() += convert_vector_body_to_space( - p, friction_thermo_langevin_rotation(langevin, p, time_step, kT)); + if (should_propagate_with(p, PropagationMode::ROT_LANGEVIN)) + p.torque() += convert_vector_body_to_space( + p, friction_thermo_langevin_rotation(langevin, p, time_step, kT)); #endif - } - } } } @@ -356,8 +348,11 @@ static bool integrator_step_1(ParticleRange const &particles, double kT, return steepest_descent_step(particles); for (auto &p : particles) { - if (should_propagate_with(p, TRANS_VS_RELATIVE)) +#ifdef VIRTUAL_SITES + // virtual sites are updated later in the integration loop + if (p.is_virtual()) continue; +#endif if (should_propagate_with(p, TRANS_LB_MOMENTUM_EXCHANGE)) velocity_verlet_propagator_1(p, time_step); if (should_propagate_with(p, TRANS_NEWTON)) @@ -409,8 +404,11 @@ static void integrator_step_2(ParticleRange const &particles, double kT, if (integ_switch == INTEG_METHOD_STEEPEST_DESCENT) return; for (auto &p : particles) { - if (should_propagate_with(p, TRANS_VS_RELATIVE)) +#ifdef VIRTUAL_SITES + // virtual sites are updated later in the integration loop + if (p.is_virtual()) continue; +#endif if (should_propagate_with(p, TRANS_LB_MOMENTUM_EXCHANGE)) velocity_verlet_propagator_2(p, time_step); if (should_propagate_with(p, TRANS_NEWTON)) @@ -551,8 +549,8 @@ int System::integrate(int n_steps, int reuse_forces) { #ifdef NPT if (integ_switch == INTEG_METHOD_NPT_ISO) { - // TODO: workaround to force a ghost update - cell_structure->update_ghosts_and_resort_particle(0xffff); + cell_structure->update_ghosts_and_resort_particle( + Cells::DATA_PART_PROPERTIES); } #endif diff --git a/src/core/unit_tests/Particle_test.cpp b/src/core/unit_tests/Particle_test.cpp index 8eef353fc6..a03c2c2e51 100644 --- a/src/core/unit_tests/Particle_test.cpp +++ b/src/core/unit_tests/Particle_test.cpp @@ -24,6 +24,7 @@ #include #include "Particle.hpp" +#include "PropagationMode.hpp" #include "config/config.hpp" #include @@ -36,6 +37,7 @@ #include #include #include +#include #include #include @@ -290,4 +292,8 @@ BOOST_AUTO_TEST_CASE(particle_bitfields) { p.set_cannot_rotate_all_axes(); BOOST_CHECK(not p.can_rotate()); #endif + + static_assert( + std::is_same_v, + decltype(ParticleProperties::propagation)>); } diff --git a/src/python/espressomd/particle_data.py b/src/python/espressomd/particle_data.py index ac6ad60a08..bde3cb5218 100644 --- a/src/python/espressomd/particle_data.py +++ b/src/python/espressomd/particle_data.py @@ -372,7 +372,6 @@ class ParticleHandle(ScriptInterfaceHelper): is_virtual() Whether the particle is a virtual site. - """ _so_name = "Particles::ParticleHandle" diff --git a/src/script_interface/particle_data/ParticleHandle.cpp b/src/script_interface/particle_data/ParticleHandle.cpp index d5f4d1ffde..fab4cb4032 100644 --- a/src/script_interface/particle_data/ParticleHandle.cpp +++ b/src/script_interface/particle_data/ParticleHandle.cpp @@ -593,14 +593,10 @@ Variant ParticleHandle::do_call_method(std::string const &name, remove_particle(m_pid); }); } else if (name == "is_virtual") { -#ifdef VIRTUAL_SITES if (not context()->is_head_node()) { return {}; } return get_particle_data(m_pid).is_virtual(); -#else - return false; -#endif // VIRTUAL_SITES #ifdef VIRTUAL_SITES_RELATIVE } else if (name == "vs_relate_to") { if (not context()->is_head_node()) { diff --git a/src/script_interface/system/System.cpp b/src/script_interface/system/System.cpp index a92c06e667..2a2eaf3e05 100644 --- a/src/script_interface/system/System.cpp +++ b/src/script_interface/system/System.cpp @@ -59,7 +59,6 @@ #include #include #include -#include #include namespace ScriptInterface { @@ -218,8 +217,7 @@ static void rotate_system(CellStructure &cell_structure, double phi, } cell_structure.set_resort_particles(Cells::RESORT_GLOBAL); - // TODO: workaround to force a ghost update - cell_structure.update_ghosts_and_resort_particle(0xffff); + cell_structure.update_ghosts_and_resort_particle(Cells::DATA_PART_PROPERTIES); } /** Rescale all particle positions in direction @p dir by a factor @p scale. diff --git a/testsuite/python/integrator_brownian.py b/testsuite/python/integrator_brownian.py index a3362c0584..1c2620ea69 100644 --- a/testsuite/python/integrator_brownian.py +++ b/testsuite/python/integrator_brownian.py @@ -181,23 +181,27 @@ def test_propagation(self): for various combinations of propagation modes. """ Propagation = espressomd.propagation.Propagation - gamma_trans = 1.2 - gamma_rot = 1.5 - x0 = [0., 1., 2.] + aniso = espressomd.has_features("PARTICLE_ANISOTROPY") + gamma_trans = np.array([1.2, 1.2, 1.2]) if aniso else 1.2 + gamma_rot = np.array([1.5, 1.5, 1.5]) if aniso else 1.5 + x0 = np.array([0., 1., 2.]) v0 = np.array([1., 2., 3.]) o0 = np.array([2., 3., 4.]) ext_force = np.array([-1., +2., -4.]) ext_torque = np.array([+1., -3., +5.]) + def has_gamma(gamma): + return gamma[0] > 0. if np.shape(gamma) else gamma > 0. + def calc_trajectory(p): t = self.system.time gamma_t = gamma_trans gamma_r = gamma_rot if espressomd.has_features("THERMOSTAT_PER_PARTICLE") and \ p.propagation & Propagation.SYSTEM_DEFAULT: - if p.gamma[0] > 0.: + if has_gamma(p.gamma): gamma_t = p.gamma - if p.gamma_rot[0] > 0.: + if has_gamma(p.gamma_rot): gamma_r = p.gamma_rot if (p.propagation & (Propagation.SYSTEM_DEFAULT | Propagation.TRANS_BROWNIAN)): diff --git a/testsuite/python/integrator_langevin.py b/testsuite/python/integrator_langevin.py index c0c8f12db3..05785c738f 100644 --- a/testsuite/python/integrator_langevin.py +++ b/testsuite/python/integrator_langevin.py @@ -174,22 +174,18 @@ def test_03__friction_rot(self): """Test the rotational friction-only part of the thermostat.""" system = self.system + aniso = espressomd.has_features("PARTICLE_ANISOTROPY") o0 = np.array([5., 5., 5.]) - if espressomd.has_features("PARTICLE_ANISOTROPY"): - gamma_t = np.array([0.5, 2., 1.5]) - gamma_r = np.array([1.5, 0.7, 1.2]) - else: - gamma_t = 2. - gamma_r = 3. + gamma_t = np.array([0.5, 2., 1.5]) if aniso else 2. + gamma_r = np.array([1.5, 0.7, 1.2]) if aniso else 3. + rinertia = np.array([1., 1., 1.]) system.time = 0. system.time_step = 0.0001 p = system.part.add(pos=(0, 0, 0), omega_body=o0, rotation=3 * [True]) if espressomd.has_features("ROTATIONAL_INERTIA"): - rinertia = np.array((2, 2, 2)) + rinertia = np.array([2., 2., 2.]) p.rinertia = rinertia - else: - rinertia = np.array((1, 1, 1)) system.thermostat.set_langevin( kT=0, gamma=gamma_t, gamma_rotation=gamma_r, seed=41) @@ -317,23 +313,27 @@ def test_propagation(self): for various combinations of propagation modes. """ Propagation = espressomd.propagation.Propagation - gamma_trans = 1.2 - gamma_rot = 0.6 - x0 = [0., 1., 2.] + aniso = espressomd.has_features("PARTICLE_ANISOTROPY") + gamma_trans = np.array([1.2, 1.2, 1.2]) if aniso else 1.2 + gamma_rot = np.array([0.6, 0.6, 0.6]) if aniso else 0.6 + x0 = np.array([0., 1., 2.]) v0 = np.array([1., 2., 3.]) o0 = np.array([1., 2., 1.]) ext_force = np.array([-2., +2., -3.]) ext_torque = np.array([+0.25, -0.5, +1.]) + def has_gamma(gamma): + return gamma[0] > 0. if np.shape(gamma) else gamma > 0. + def calc_trajectory(p): t = self.system.time gamma_t = gamma_trans gamma_r = gamma_rot if espressomd.has_features("THERMOSTAT_PER_PARTICLE") and \ p.propagation & Propagation.SYSTEM_DEFAULT: - if p.gamma[0] > 0.: + if has_gamma(p.gamma): gamma_t = p.gamma - if p.gamma_rot[0] > 0.: + if has_gamma(p.gamma_rot): gamma_r = p.gamma_rot if (p.propagation & (Propagation.SYSTEM_DEFAULT | Propagation.TRANS_LANGEVIN)): diff --git a/testsuite/python/integrator_langevin_stats.py b/testsuite/python/integrator_langevin_stats.py index 6b4f5cf6be..882be1f709 100644 --- a/testsuite/python/integrator_langevin_stats.py +++ b/testsuite/python/integrator_langevin_stats.py @@ -117,6 +117,10 @@ def verify_diffusion(self, p, corr, kT, gamma): acf = c.result() tau = c.lag_times() + gamma = np.copy(gamma) + if not gamma.shape: + gamma = gamma * np.ones(3) + # Integrate with trapezoidal rule for i in range(3): I = np.trapz(acf[:, p.id, i], tau) @@ -130,21 +134,20 @@ def test_06__diffusion(self): kT = 1.37 dt = 0.1 system.time_step = dt + aniso = espressomd.has_features("PARTICLE_ANISOTROPY") # Translational gamma. We cannot test per-component, if rotation is on, # because body and space frames become different. gamma = 3.1 # Rotational gamma - gamma_rot_i = 4.7 - gamma_rot_a = [4.2, 1, 1.2] + gamma_rot = [4.2, 1, 1.2] if aniso else 4.7 # If we have langevin per particle: # Translation per_part_gamma = 1.63 # Rotational - per_part_gamma_rot_i = 2.6 - per_part_gamma_rot_a = [2.4, 3.8, 1.1] + per_part_gamma_rot = [2.4, 3.8, 1.1] if aniso else 2.6 # Particle with global thermostat params p_global = system.part.add(pos=(0, 0, 0)) @@ -155,25 +158,14 @@ def test_06__diffusion(self): if espressomd.has_features("THERMOSTAT_PER_PARTICLE"): p_gamma = system.part.add(pos=(0, 0, 0)) self.setup_diff_mass_rinertia(p_gamma) - if espressomd.has_features("PARTICLE_ANISOTROPY"): - p_gamma.gamma = 3 * [per_part_gamma] - if espressomd.has_features("ROTATION"): - p_gamma.gamma_rot = per_part_gamma_rot_a - else: - p_gamma.gamma = per_part_gamma - if espressomd.has_features("ROTATION"): - p_gamma.gamma_rot = per_part_gamma_rot_i + p_gamma.gamma = 3 * [per_part_gamma] if aniso else per_part_gamma + if espressomd.has_features("ROTATION"): + p_gamma.gamma_rot = per_part_gamma_rot # Thermostat setup if espressomd.has_features("ROTATION"): - if espressomd.has_features("PARTICLE_ANISOTROPY"): - # particle anisotropy and rotation - system.thermostat.set_langevin( - kT=kT, gamma=gamma, gamma_rotation=gamma_rot_a, seed=41) - else: - # Rotation without particle anisotropy - system.thermostat.set_langevin( - kT=kT, gamma=gamma, gamma_rotation=gamma_rot_i, seed=41) + system.thermostat.set_langevin( + kT=kT, gamma=gamma, gamma_rotation=gamma_rot, seed=41) else: # No rotation system.thermostat.set_langevin(kT=kT, gamma=gamma, seed=41) @@ -215,31 +207,15 @@ def test_06__diffusion(self): corr_omega.finalize() # Verify diffusion - # Translation - # Cast gammas to vector, to make checks independent of - # PARTICLE_ANISOTROPY - gamma = np.ones(3) * gamma - per_part_gamma = np.ones(3) * per_part_gamma self.verify_diffusion(p_global, corr_vel, kT, gamma) if espressomd.has_features("THERMOSTAT_PER_PARTICLE"): self.verify_diffusion(p_gamma, corr_vel, kT, per_part_gamma) - # Rotation if espressomd.has_features("ROTATION"): - # Decide on effective gamma rotation, since for rotation it is - # direction dependent - eff_gamma_rot = None - if espressomd.has_features("PARTICLE_ANISOTROPY"): - eff_gamma_rot = gamma_rot_a - eff_per_part_gamma_rot = per_part_gamma_rot_a - else: - eff_gamma_rot = gamma_rot_i * np.ones(3) - eff_per_part_gamma_rot = per_part_gamma_rot_i * np.ones(3) - - self.verify_diffusion(p_global, corr_omega, kT, eff_gamma_rot) + self.verify_diffusion(p_global, corr_omega, kT, gamma_rot) if espressomd.has_features("THERMOSTAT_PER_PARTICLE"): self.verify_diffusion( - p_gamma, corr_omega, kT, eff_per_part_gamma_rot) + p_gamma, corr_omega, kT, per_part_gamma_rot) def test_08__noise_correlation(self): """Checks that the Langevin noise is uncorrelated""" diff --git a/testsuite/python/integrator_npt.py b/testsuite/python/integrator_npt.py index d9fd4df49d..15ec84be5c 100644 --- a/testsuite/python/integrator_npt.py +++ b/testsuite/python/integrator_npt.py @@ -101,6 +101,7 @@ def calc_trajectory(p, x0): modes_trans = [ Propagation.NONE, Propagation.SYSTEM_DEFAULT, + Propagation.TRANS_NEWTON, Propagation.TRANS_LANGEVIN_NPT, ] for mode_trans in modes_trans: