From c2aa5437495c0c27b6fac99906f94e72b8db787a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 17 Dec 2020 16:28:42 +0100 Subject: [PATCH] core: Remove per-particle temperature feature The feature was introduced for visualization purposes and is no longer used. Improper use of that feature can lead to unphysical results. It also makes testing the thermostats a lot more difficult. --- doc/sphinx/advanced_methods.rst | 4 +- doc/sphinx/installation.rst | 4 +- doc/sphinx/system_setup.rst | 4 +- samples/chamber_game.py | 11 +-- samples/drude_bmimpf6.py | 2 +- src/core/Particle.hpp | 2 - src/core/particle_data.cpp | 7 -- src/core/particle_data.hpp | 7 -- src/core/thermostats/brownian_inline.hpp | 97 ++----------------- src/core/thermostats/langevin_inline.hpp | 15 ++- src/python/espressomd/drude_helpers.py | 3 +- src/python/espressomd/particle_data.pxd | 3 - src/python/espressomd/particle_data.pyx | 23 ----- testsuite/python/brownian_dynamics.py | 11 +-- testsuite/python/brownian_dynamics_stats.py | 4 +- testsuite/python/langevin_thermostat.py | 11 +-- testsuite/python/langevin_thermostat_stats.py | 31 +----- .../python/mass-and-rinertia_per_particle.py | 91 ++--------------- testsuite/python/thermostats_common.py | 65 ++++++------- 19 files changed, 68 insertions(+), 327 deletions(-) diff --git a/doc/sphinx/advanced_methods.rst b/doc/sphinx/advanced_methods.rst index 0804952a408..adc43773aa0 100644 --- a/doc/sphinx/advanced_methods.rst +++ b/doc/sphinx/advanced_methods.rst @@ -1498,7 +1498,7 @@ charge) to a particle (Drude core) that mimics an electron cloud which can be elongated to create a dynamically inducible dipole. The energetic minimum of the Drude charge can be obtained self-consistently, which requires several iterations of the system's electrostatics and is usually considered -computational expensive. However, with thermalized cold Drude oscillators, the +computationally expensive. However, with thermalized cold Drude oscillators, the distance between Drude charge and core is coupled to a thermostat so that it fluctuates around the SCF solution. This thermostat is kept at a low temperature compared to the global temperature to minimize the heat flow into @@ -1516,7 +1516,7 @@ In |es|, the basic ingredients to simulate such a system are split into three bo The system-wide thermostat has to be applied to the centre of mass and not to the core particle directly. Therefore, the particles have to be excluded from global thermostatting. With ``LANGEVIN_PER_PARTICLE`` enabled, we set the -temperature and friction coefficient of the Drude complex to zero, which allows +friction coefficient of the Drude complex to zero, which allows to still use a global Langevin thermostat for non-polarizable particles. As the Drude charge should not alter the *charge* or *mass* of the Drude diff --git a/doc/sphinx/installation.rst b/doc/sphinx/installation.rst index 20c197557c9..257b7f280f0 100644 --- a/doc/sphinx/installation.rst +++ b/doc/sphinx/installation.rst @@ -315,10 +315,10 @@ General features additional degrees of freedom, which for example means that the kinetic energy changes at constant temperature is twice as large. -- ``LANGEVIN_PER_PARTICLE`` Allows to choose the Langevin temperature and friction coefficient +- ``LANGEVIN_PER_PARTICLE`` Allows to choose the Langevin friction coefficient per particle. -- ``BROWNIAN_PER_PARTICLE`` Allows to choose the Brownian temperature and friction coefficient +- ``BROWNIAN_PER_PARTICLE`` Allows to choose the Brownian friction coefficient per particle. - ``ROTATIONAL_INERTIA`` diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index b6bf671bad6..b3b644eee47 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -276,7 +276,7 @@ or, with feature ``PARTICLE_ANISOTROPY`` compiled in, as the three eigenvalues of the respective friction coefficient tensor. This is enables the simulation of the anisotropic diffusion of anisotropic colloids (rods, etc.). -Using the Langevin thermostat, it is possible to set a temperature and a +Using the Langevin thermostat, it is possible to set a friction coefficient for every particle individually via the feature ``LANGEVIN_PER_PARTICLE``. Consult the reference of the ``part`` command (chapter :ref:`Setting up particles`) for information on how to achieve this. @@ -427,7 +427,7 @@ where ``gamma`` (hereinafter :math:`\gamma`) is a viscous friction coefficient. In terms of the Python interface and setup, the Brownian thermostat is very similar to the :ref:`Langevin thermostat`. The feature ``BROWNIAN_PER_PARTICLE`` is used to control the per-particle -temperature and the friction coefficient setup. The major differences are +friction coefficient setup. The major differences are its internal integrator implementation and other temporal constraints. The integrator is still a symplectic Velocity Verlet-like one. It is implemented via a viscous drag part and a random walk of both the position and diff --git a/samples/chamber_game.py b/samples/chamber_game.py index 55ae8f40891..14b4e14ee0c 100644 --- a/samples/chamber_game.py +++ b/samples/chamber_game.py @@ -61,13 +61,11 @@ # PARAMETERS # PHYSICS -temperature_snake = 0.0 gamma_snake_head = 1.0 gamma_snake_bead = 15.0 -temperature_bubbles = 10000.0 -temp_l = temperature_bubbles -temp_r = temperature_bubbles +temp_l = 10000.0 +temp_r = temp_l temp_max = 1e5 gamma_bubbles = 0.5 @@ -180,7 +178,6 @@ type=snake_head_type, fix=[False, False, True], mass=snake_head_mass, - temp=temperature_snake, gamma=gamma_snake_head) else: system.part.add( @@ -192,7 +189,6 @@ type=snake_bead_type, fix=[False, False, True], mass=snake_bead_mass, - temp=temperature_snake, gamma=gamma_snake_bead) # NB INTER @@ -262,7 +258,6 @@ type=bubble_type, fix=[False, False, True], mass=bubble_mass, - temp=temperature_bubbles, gamma=gamma_bubbles) testid = len(system.part) - 1 n += 1 @@ -453,10 +448,8 @@ def T_to_g(temp): Nl = len(pl) Nr = len(pr) for p in pl: - p.temp = temp_l p.gamma = T_to_g(temp_l) for p in pr: - p.temp = temp_r p.gamma = T_to_g(temp_r) w = visualizer.specs['window_size'] diff --git a/samples/drude_bmimpf6.py b/samples/drude_bmimpf6.py index 2b2767e07a0..b8e357a9c59 100644 --- a/samples/drude_bmimpf6.py +++ b/samples/drude_bmimpf6.py @@ -215,7 +215,7 @@ def combination_rule_sigma(rule, sig1, sig2): system.part.add( id=rid, type=types["BMIM_COM"], pos=pos_com, mass=masses["BMIM_COM"], rinertia=[646.284, 585.158, 61.126], - temp=0, gamma=0, rotation=[1, 1, 1]) + gamma=0, rotation=[1, 1, 1]) cation_com_ids.append(rid) com_id = rid rid += 1 diff --git a/src/core/Particle.hpp b/src/core/Particle.hpp index 2ea25f7c521..6f61c3e93e7 100644 --- a/src/core/Particle.hpp +++ b/src/core/Particle.hpp @@ -147,7 +147,6 @@ struct ParticleProperties { #endif /* VIRTUAL_SITES */ #if defined(LANGEVIN_PER_PARTICLE) || defined(BROWNIAN_PER_PARTICLE) - double T = -1.; #ifndef PARTICLE_ANISOTROPY double gamma = -1.; #else @@ -221,7 +220,6 @@ struct ParticleProperties { #endif /* VIRTUAL_SITES */ #if defined(LANGEVIN_PER_PARTICLE) || defined(BROWNIAN_PER_PARTICLE) - ar &T; ar γ #ifdef ROTATION ar &gamma_rot; diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp index f9ef2dd94a5..6a3e6625369 100644 --- a/src/core/particle_data.cpp +++ b/src/core/particle_data.cpp @@ -122,7 +122,6 @@ using UpdatePropertyMessage = boost::variant #endif #endif #if defined(LANGEVIN_PER_PARTICLE) || defined(BROWNIAN_PER_PARTICLE) - , UpdateProperty #ifndef PARTICLE_ANISOTROPY , UpdateProperty #else @@ -867,9 +866,6 @@ void set_particle_torque_lab(int part, const Utils::Vector3d &torque_lab) { #endif #if defined(LANGEVIN_PER_PARTICLE) || defined(BROWNIAN_PER_PARTICLE) -void set_particle_temperature(int part, double T) { - mpi_update_particle_property(part, T); -} #ifndef PARTICLE_ANISOTROPY void set_particle_gamma(int part, double gamma) { @@ -1264,9 +1260,6 @@ void pointer_to_gamma_rot(Particle const *p, double const *&res) { } #endif // ROTATION -void pointer_to_temperature(Particle const *p, double const *&res) { - res = &(p->p.T); -} #endif // LANGEVIN_PER_PARTICLE || BROWNIAN_PER_PARTICLE #ifdef ENGINE diff --git a/src/core/particle_data.hpp b/src/core/particle_data.hpp index 0d07c4cea93..f2f078b931e 100644 --- a/src/core/particle_data.hpp +++ b/src/core/particle_data.hpp @@ -243,12 +243,6 @@ void set_particle_vs_relative(int part, int vs_relative_to, double vs_distance, #endif #if defined(LANGEVIN_PER_PARTICLE) || defined(BROWNIAN_PER_PARTICLE) -/** Call only on the master node: set particle temperature. - * @param part the particle. - * @param T its new temperature. - */ -void set_particle_temperature(int part, double T); - /** Call only on the master node: set particle frictional coefficient. * @param part the particle. * @param gamma its new frictional coefficient. @@ -393,7 +387,6 @@ void pointer_to_fix(Particle const *p, const uint8_t *&res); #if defined(LANGEVIN_PER_PARTICLE) || defined(BROWNIAN_PER_PARTICLE) void pointer_to_gamma(Particle const *p, double const *&res); -void pointer_to_temperature(Particle const *p, double const *&res); #ifdef ROTATION void pointer_to_gamma_rot(Particle const *p, double const *&res); #endif diff --git a/src/core/thermostats/brownian_inline.hpp b/src/core/thermostats/brownian_inline.hpp index ec76b95521d..e0fc9348da0 100644 --- a/src/core/thermostats/brownian_inline.hpp +++ b/src/core/thermostats/brownian_inline.hpp @@ -170,41 +170,16 @@ inline Utils::Vector3d bd_random_walk(BrownianThermostat const &brownian, // skip the translation thermalizing for virtual sites unless enabled if (p.p.is_virtual && !thermo_virtual) return {}; - // first, set defaults - Thermostat::GammaType sigma_pos = brownian.sigma_pos; - // Override defaults if per-particle values for T and gamma are given + Thermostat::GammaType sigma_pos = brownian.sigma_pos; #ifdef BROWNIAN_PER_PARTICLE + // override default if particle-specific gamma if (p.p.gamma >= Thermostat::GammaType{}) { - // Is a particle-specific temperature also specified? - if (p.p.T >= 0.) { - if (p.p.T > 0.0) { - sigma_pos = BrownianThermostat::sigma(p.p.T, p.p.gamma); - } else { - // just an indication of the infinity - sigma_pos = Thermostat::GammaType{}; - } - } else - // default temperature but particle-specific gamma - if (temperature > 0.0) { + if (temperature > 0.0) { sigma_pos = BrownianThermostat::sigma(temperature, p.p.gamma); } else { sigma_pos = Thermostat::GammaType{}; } - } // particle-specific gamma - else { - // No particle-specific gamma, but is there particle-specific temperature - if (p.p.T >= 0.) { - if (p.p.T > 0.0) { - sigma_pos = BrownianThermostat::sigma(p.p.T, brownian.gamma); - } else { - // just an indication of the infinity - sigma_pos = Thermostat::GammaType{}; - } - } else { - // default values for both - sigma_pos = brownian.sigma_pos; - } } #endif // BROWNIAN_PER_PARTICLE @@ -270,22 +245,6 @@ inline Utils::Vector3d bd_random_walk_vel(BrownianThermostat const &brownian, // skip the translation thermalizing for virtual sites unless enabled if (p.p.is_virtual && !thermo_virtual) return {}; - // Just a square root of kT, see eq. (10.2.17) and comments in 2 paragraphs - // afterwards, Pottier2010 - double sigma_vel; - - // Override defaults if per-particle values for T and gamma are given -#ifdef BROWNIAN_PER_PARTICLE - // Is a particle-specific temperature specified? - if (p.p.T >= 0.) { - sigma_vel = BrownianThermostat::sigma(p.p.T); - } else { - sigma_vel = brownian.sigma_vel; - } -#else - // defaults - sigma_vel = brownian.sigma_vel; -#endif // BROWNIAN_PER_PARTICLE auto const noise = Random::noise_gaussian( brownian.rng_counter(), brownian.rng_seed(), p.identity()); @@ -302,7 +261,7 @@ inline Utils::Vector3d bd_random_walk_vel(BrownianThermostat const &brownian, // (14.31) of Schlick2010. A difference is the mass factor to the friction // tensor. The noise is Gaussian according to the convention at p. 237 // (last paragraph), Pottier2010. - velocity[j] += sigma_vel * noise[j] / sqrt(p.p.mass); + velocity[j] += brownian.sigma_vel * noise[j] / sqrt(p.p.mass); } } return velocity; @@ -398,40 +357,15 @@ bd_drag_vel_rot(Thermostat::GammaType const &brownian_gamma_rotation, inline Utils::Quaternion bd_random_walk_rot(BrownianThermostat const &brownian, Particle const &p, double dt) { - // first, set defaults - Thermostat::GammaType sigma_pos = brownian.sigma_pos_rotation; - // Override defaults if per-particle values for T and gamma are given + Thermostat::GammaType sigma_pos = brownian.sigma_pos_rotation; #ifdef BROWNIAN_PER_PARTICLE + // override default if particle-specific gamma if (p.p.gamma_rot >= Thermostat::GammaType{}) { - // Is a particle-specific temperature also specified? - if (p.p.T >= 0.) { - if (p.p.T > 0.0) { - sigma_pos = BrownianThermostat::sigma(p.p.T, p.p.gamma_rot); - } else { - // just an indication of the infinity - sigma_pos = {}; - } - } else if (temperature > 0.) { - // Default temperature but particle-specific gamma + if (temperature > 0.) { sigma_pos = BrownianThermostat::sigma(temperature, p.p.gamma_rot); } else { - // just an indication of the infinity - sigma_pos = {}; - } - } // particle-specific gamma - else { - // No particle-specific gamma, but is there particle-specific temperature - if (p.p.T >= 0.) { - if (p.p.T > 0.0) { - sigma_pos = BrownianThermostat::sigma(p.p.T, brownian.gamma_rotation); - } else { - // just an indication of the infinity - sigma_pos = {}; - } - } else { - // Default values for both - sigma_pos = brownian.sigma_pos_rotation; + sigma_pos = {}; // just an indication of the infinity } } #endif // BROWNIAN_PER_PARTICLE @@ -472,20 +406,7 @@ bd_random_walk_rot(BrownianThermostat const &brownian, Particle const &p, */ inline Utils::Vector3d bd_random_walk_vel_rot(BrownianThermostat const &brownian, Particle const &p) { - double sigma_vel; - - // Override defaults if per-particle values for T and gamma are given -#ifdef BROWNIAN_PER_PARTICLE - // Is a particle-specific temperature specified? - if (p.p.T >= 0.) { - sigma_vel = BrownianThermostat::sigma(p.p.T); - } else { - sigma_vel = brownian.sigma_vel_rotation; - } -#else - // set defaults - sigma_vel = brownian.sigma_vel_rotation; -#endif // BROWNIAN_PER_PARTICLE + auto const sigma_vel = brownian.sigma_vel_rotation; Utils::Vector3d domega{}; auto const noise = Random::noise_gaussian( diff --git a/src/core/thermostats/langevin_inline.hpp b/src/core/thermostats/langevin_inline.hpp index 273b3615552..30dc31dee6e 100644 --- a/src/core/thermostats/langevin_inline.hpp +++ b/src/core/thermostats/langevin_inline.hpp @@ -48,17 +48,15 @@ friction_thermo_langevin(LangevinThermostat const &langevin, Particle const &p, } // Determine prefactors for the friction and the noise term - // first, set defaults Thermostat::GammaType pref_friction = langevin.pref_friction; Thermostat::GammaType pref_noise = langevin.pref_noise; - // Override defaults if per-particle values for T and gamma are given #ifdef LANGEVIN_PER_PARTICLE - if (p.p.gamma >= Thermostat::GammaType{} or p.p.T >= 0.) { - auto const kT = p.p.T >= 0. ? p.p.T : temperature; + // override default if particle-specific gamma + if (p.p.gamma >= Thermostat::GammaType{}) { auto const gamma = p.p.gamma >= Thermostat::GammaType{} ? p.p.gamma : langevin.gamma; pref_friction = -gamma; - pref_noise = LangevinThermostat::sigma(kT, time_step, gamma); + pref_noise = LangevinThermostat::sigma(temperature, time_step, gamma); } #endif // LANGEVIN_PER_PARTICLE @@ -108,15 +106,14 @@ friction_thermo_langevin_rotation(LangevinThermostat const &langevin, auto pref_friction = -langevin.gamma_rotation; auto pref_noise = langevin.pref_noise_rotation; - // Override defaults if per-particle values for T and gamma are given #ifdef LANGEVIN_PER_PARTICLE - if (p.p.gamma_rot >= Thermostat::GammaType{} or p.p.T >= 0.) { - auto const kT = p.p.T >= 0. ? p.p.T : temperature; + // override default if particle-specific gamma + if (p.p.gamma_rot >= Thermostat::GammaType{}) { auto const gamma = p.p.gamma_rot >= Thermostat::GammaType{} ? p.p.gamma_rot : langevin.gamma_rotation; pref_friction = -gamma; - pref_noise = LangevinThermostat::sigma(kT, time_step, gamma); + pref_noise = LangevinThermostat::sigma(temperature, time_step, gamma); } #endif // LANGEVIN_PER_PARTICLE diff --git a/src/python/espressomd/drude_helpers.py b/src/python/espressomd/drude_helpers.py index 705d81357cc..52f588b1e42 100644 --- a/src/python/espressomd/drude_helpers.py +++ b/src/python/espressomd/drude_helpers.py @@ -76,7 +76,7 @@ def add_drude_particle_to_core(system, harmonic_bond, thermalized_bond, gamma_off = 0.0 system.part.add(id=id_drude, pos=p_core.pos, type=type_drude, - q=q_drude, mass=mass_drude, temp=0, gamma=gamma_off) + q=q_drude, mass=mass_drude, gamma=gamma_off) if verbose: print( @@ -88,7 +88,6 @@ def add_drude_particle_to_core(system, harmonic_bond, thermalized_bond, p_core.add_bond((harmonic_bond, id_drude)) p_core.add_bond((thermalized_bond, id_drude)) - p_core.temp = 0.0 p_core.gamma = gamma_off if type_drude in drude_dict and not ( diff --git a/src/python/espressomd/particle_data.pxd b/src/python/espressomd/particle_data.pxd index 521b99a91e2..390fe9003af 100644 --- a/src/python/espressomd/particle_data.pxd +++ b/src/python/espressomd/particle_data.pxd @@ -130,9 +130,6 @@ cdef extern from "particle_data.hpp": void pointer_to_virtual(const particle * P, const bint * & res) IF LANGEVIN_PER_PARTICLE or BROWNIAN_PER_PARTICLE: - void set_particle_temperature(int part, double T) - void pointer_to_temperature(const particle * p, const double * & res) - IF PARTICLE_ANISOTROPY: void set_particle_gamma(int part, Vector3d gamma) ELSE: diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx index 8a946b34fad..a7e4000b4cb 100644 --- a/src/python/espressomd/particle_data.pyx +++ b/src/python/espressomd/particle_data.pyx @@ -969,29 +969,6 @@ cdef class ParticleHandle: self.particle_data, gamma_rot) return gamma_rot[0] - property temp: - """ - Particle's temperature in the Langevin and Brownian thermostats. - - temp: :obj:`float` - - .. note:: - This needs the feature ``LANGEVIN_PER_PARTICLE`` or - ``BROWNIAN_PER_PARTICLE``. - - """ - - def __set__(self, _temp): - check_type_or_throw_except( - _temp, 1, float, "temp has to be a float.") - set_particle_temperature(self._id, _temp) - - def __get__(self): - self.update_particle_data() - cdef const double * temp = NULL - pointer_to_temperature(self.particle_data, temp) - return temp[0] - IF ROTATION: property rotation: """ diff --git a/testsuite/python/brownian_dynamics.py b/testsuite/python/brownian_dynamics.py index c14dd5fb9ab..a341749242d 100644 --- a/testsuite/python/brownian_dynamics.py +++ b/testsuite/python/brownian_dynamics.py @@ -40,7 +40,7 @@ def tearDown(self): self.system.thermostat.turn_off() self.system.integrator.set_vv() - def check_rng(self, per_particle_gamma=False, per_particle_temp=False): + def check_rng(self, per_particle_gamma=False): """Test for RNG consistency.""" kT = 1.1 @@ -51,17 +51,14 @@ def reset_particle(): p = system.part.add(pos=[0, 0, 0]) if espressomd.has_features("ROTATION"): p.rotation = [1, 1, 1] - if per_particle_gamma or per_particle_temp: - assert espressomd.has_features("BROWNIAN_PER_PARTICLE") if per_particle_gamma: + assert espressomd.has_features("BROWNIAN_PER_PARTICLE") if espressomd.has_features("PARTICLE_ANISOTROPY"): p.gamma = 3 * [gamma / 2] else: p.gamma = gamma / 2 if espressomd.has_features("ROTATION"): p.gamma_rot = p.gamma * 1.5 - if per_particle_temp: - p.temp = 2 * kT return p system = self.system @@ -128,9 +125,7 @@ def test_01__rng(self): @utx.skipIfMissingFeatures("BROWNIAN_PER_PARTICLE") def test_01__rng_per_particle(self): """Test for RNG consistency.""" - self.check_rng(False, True) - self.check_rng(True, False) - self.check_rng(True, True) + self.check_rng(True) @utx.skipIfMissingFeatures("VIRTUAL_SITES") def test_07__virtual(self): diff --git a/testsuite/python/brownian_dynamics_stats.py b/testsuite/python/brownian_dynamics_stats.py index d5433a24d50..fe9b61d66cd 100644 --- a/testsuite/python/brownian_dynamics_stats.py +++ b/testsuite/python/brownian_dynamics_stats.py @@ -88,14 +88,12 @@ def test_vel_dist_per_particle(self): kT = 0.9 gamma = 3.2 gamma2 = 4.3 - kT2 = 1.5 system.thermostat.set_brownian(kT=kT, gamma=gamma, seed=41) loops = 200 v_minmax = 5 bins = 4 error_tol = 0.012 - self.check_per_particle( - N, kT, kT2, gamma2, loops, v_minmax, bins, error_tol) + self.check_per_particle(N, kT, gamma2, loops, v_minmax, bins, error_tol) def setup_diff_mass_rinertia(self, p): if espressomd.has_features("MASS"): diff --git a/testsuite/python/langevin_thermostat.py b/testsuite/python/langevin_thermostat.py index e1e033053ae..7dbbd32e154 100644 --- a/testsuite/python/langevin_thermostat.py +++ b/testsuite/python/langevin_thermostat.py @@ -40,7 +40,7 @@ def tearDown(self): self.system.thermostat.turn_off() self.system.integrator.set_vv() - def check_rng(self, per_particle_gamma=False, per_particle_temp=False): + def check_rng(self, per_particle_gamma=False): """Test for RNG consistency.""" kT = 1.1 @@ -51,17 +51,14 @@ def reset_particle(): p = system.part.add(pos=[0, 0, 0]) if espressomd.has_features("ROTATION"): p.rotation = [1, 1, 1] - if per_particle_gamma or per_particle_temp: - assert espressomd.has_features("LANGEVIN_PER_PARTICLE") if per_particle_gamma: + assert espressomd.has_features("LANGEVIN_PER_PARTICLE") if espressomd.has_features("PARTICLE_ANISOTROPY"): p.gamma = 3 * [gamma / 2] else: p.gamma = gamma / 2 if espressomd.has_features("ROTATION"): p.gamma_rot = p.gamma * 1.5 - if per_particle_temp: - p.temp = 2 * kT return p system = self.system @@ -143,9 +140,7 @@ def test_01__rng(self): @utx.skipIfMissingFeatures("LANGEVIN_PER_PARTICLE") def test_01__rng_per_particle(self): """Test for RNG consistency.""" - self.check_rng(False, True) - self.check_rng(True, False) - self.check_rng(True, True) + self.check_rng(True) def test_02__friction_trans(self): """Tests the translational friction-only part of the thermostat.""" diff --git a/testsuite/python/langevin_thermostat_stats.py b/testsuite/python/langevin_thermostat_stats.py index bfe4d865797..eef5f4e70c5 100644 --- a/testsuite/python/langevin_thermostat_stats.py +++ b/testsuite/python/langevin_thermostat_stats.py @@ -87,14 +87,12 @@ def test_vel_dist_per_particle(self): kT = 0.9 gamma = 3.2 gamma2 = 4.3 - kT2 = 1.5 system.thermostat.set_langevin(kT=kT, gamma=gamma, seed=41) loops = 300 v_minmax = 5 bins = 4 error_tol = 0.016 - self.check_per_particle( - N, kT, kT2, gamma2, loops, v_minmax, bins, error_tol) + self.check_per_particle(N, kT, gamma2, loops, v_minmax, bins, error_tol) def setup_diff_mass_rinertia(self, p): if espressomd.has_features("MASS"): @@ -140,8 +138,6 @@ def test_06__diffusion(self): gamma_rot_a = [4.2, 1, 1.2] # If we have langevin per particle: - # per particle kT - per_part_kT = 1.6 # Translation per_part_gamma = 1.63 # Rotational @@ -166,22 +162,6 @@ def test_06__diffusion(self): if espressomd.has_features("ROTATION"): p_gamma.gamma_rot = per_part_gamma_rot_i - p_kT = system.part.add(pos=(0, 0, 0)) - self.setup_diff_mass_rinertia(p_kT) - p_kT.temp = per_part_kT - - p_both = system.part.add(pos=(0, 0, 0)) - self.setup_diff_mass_rinertia(p_both) - p_both.temp = per_part_kT - if espressomd.has_features("PARTICLE_ANISOTROPY"): - p_both.gamma = per_part_gamma, per_part_gamma, per_part_gamma - if espressomd.has_features("ROTATION"): - p_both.gamma_rot = per_part_gamma_rot_a - else: - p_both.gamma = per_part_gamma - if espressomd.has_features("ROTATION"): - p_both.gamma_rot = per_part_gamma_rot_i - # Thermostat setup if espressomd.has_features("ROTATION"): if espressomd.has_features("PARTICLE_ANISOTROPY"): @@ -207,8 +187,6 @@ def test_06__diffusion(self): all_particles = [p_global] if espressomd.has_features("LANGEVIN_PER_PARTICLE"): all_particles.append(p_gamma) - all_particles.append(p_kT) - all_particles.append(p_both) # linear vel vel_obs = ParticleVelocities(ids=system.part[:].id) @@ -241,9 +219,6 @@ def test_06__diffusion(self): self.verify_diffusion(p_global, corr_vel, kT, gamma) if espressomd.has_features("LANGEVIN_PER_PARTICLE"): self.verify_diffusion(p_gamma, corr_vel, kT, per_part_gamma) - self.verify_diffusion(p_kT, corr_vel, per_part_kT, gamma) - self.verify_diffusion( - p_both, corr_vel, per_part_kT, per_part_gamma) # Rotation if espressomd.has_features("ROTATION"): @@ -261,10 +236,6 @@ def test_06__diffusion(self): if espressomd.has_features("LANGEVIN_PER_PARTICLE"): self.verify_diffusion( p_gamma, corr_omega, kT, eff_per_part_gamma_rot) - self.verify_diffusion( - p_kT, corr_omega, per_part_kT, eff_gamma_rot) - self.verify_diffusion(p_both, corr_omega, - per_part_kT, eff_per_part_gamma_rot) def test_08__noise_correlation(self): """Checks that the Langevin noise is uncorrelated""" diff --git a/testsuite/python/mass-and-rinertia_per_particle.py b/testsuite/python/mass-and-rinertia_per_particle.py index 38bd60a07f2..21d94a3f754 100644 --- a/testsuite/python/mass-and-rinertia_per_particle.py +++ b/testsuite/python/mass-and-rinertia_per_particle.py @@ -235,7 +235,7 @@ def check_dissipation(self, n): """ - tol = 1.2E-3 + tol = 1.3E-3 for _ in range(100): self.system.integrator.run(2) for i in range(n): @@ -384,23 +384,6 @@ def set_particle_specific_gamma(self, n): if espressomd.has_features("ROTATION"): self.system.part[ind].gamma_rot = self.gamma_rot_p[k, :] - def set_particle_specific_temperature(self, n): - """ - Set the particle-specific temperature. - - Parameters - ---------- - n : :obj:`int` - Number of particles of the each type. There are 2 types. - - """ - - for k in range(2): - self.halfkT_p_validate[k] = self.kT_p[k] / 2.0 - for i in range(n): - ind = i + k * n - self.system.part[ind].temp = self.kT_p[k] - def set_diffusivity_tran(self): """ Set the translational diffusivity to validate further. @@ -412,7 +395,7 @@ def set_diffusivity_tran(self): self.D_tran_p_validate[k, :] = 2.0 * \ self.halfkT_p_validate[k] / self.gamma_tran_p_validate[k, :] - # Test case 0.0: no particle specific values / dissipation only + # Test case 0.0: no particle specific gamma / dissipation only def test_case_00(self): # Each of 2 kind of particles will be represented by n instances: n = 1 @@ -424,7 +407,7 @@ def test_case_00(self): # Actual integration and validation run self.check_dissipation(n) - # Test case 0.1: no particle specific values / fluctuation & dissipation + # Test case 0.1: no particle specific gamma / fluctuation & dissipation def test_case_01(self): # Each of 2 kind of particles will be represented by n instances: n = 500 @@ -437,7 +420,7 @@ def test_case_01(self): # Actual integration and validation run self.check_fluctuation_dissipation(n) - # Test case 1.0: particle specific gamma but not temperature / dissipation + # Test case 1.0: particle specific gamma / dissipation # only def test_case_10(self): # Each of 2 kind of particles will be represented by n instances: @@ -451,7 +434,7 @@ def test_case_10(self): # Actual integration and validation run self.check_dissipation(n) - # Test case 1.1: particle specific gamma but not temperature / fluctuation + # Test case 1.1: particle specific gamma / fluctuation # & dissipation def test_case_11(self): # Each of 2 kind of particles will be represented by n instances: @@ -466,67 +449,7 @@ def test_case_11(self): # Actual integration and validation run self.check_fluctuation_dissipation(n) - # Test case 2.0: particle specific temperature but not gamma / dissipation - # only - def test_case_20(self): - # Each of 2 kind of particles will be represented by n instances: - n = 1 - self.dissipation_param_setup(n) - self.set_langevin_global_defaults() - # The test case-specific thermostat and per-particle parameters - self.system.thermostat.set_langevin( - kT=self.kT, gamma=self.gamma_global, seed=42) - self.set_particle_specific_temperature(n) - # Actual integration and validation run - self.check_dissipation(n) - - # Test case 2.1: particle specific temperature but not gamma / fluctuation - # & dissipation - def test_case_21(self): - # Each of 2 kind of particles will be represented by n instances: - n = 500 - self.fluctuation_dissipation_param_setup(n) - self.set_langevin_global_defaults() - # The test case-specific thermostat and per-particle parameters - self.system.thermostat.set_langevin( - kT=self.kT, gamma=self.gamma_global, seed=42) - self.set_particle_specific_temperature(n) - self.set_diffusivity_tran() - # Actual integration and validation run - self.check_fluctuation_dissipation(n) - - # Test case 3.0: both particle specific gamma and temperature / - # dissipation only - def test_case_30(self): - # Each of 2 kind of particles will be represented by n instances: - n = 1 - self.dissipation_param_setup(n) - self.set_langevin_global_defaults() - # The test case-specific thermostat and per-particle parameters - self.system.thermostat.set_langevin( - kT=self.kT, gamma=self.gamma_global, seed=42) - self.set_particle_specific_gamma(n) - self.set_particle_specific_temperature(n) - # Actual integration and validation run - self.check_dissipation(n) - - # Test case 3.1: both particle specific gamma and temperature / - # fluctuation & dissipation - def test_case_31(self): - # Each of 2 kind of particles will be represented by n instances: - n = 500 - self.fluctuation_dissipation_param_setup(n) - self.set_langevin_global_defaults() - # The test case-specific thermostat and per-particle parameters - self.system.thermostat.set_langevin( - kT=self.kT, gamma=self.gamma_global, seed=42) - self.set_particle_specific_gamma(n) - self.set_particle_specific_temperature(n) - self.set_diffusivity_tran() - # Actual integration and validation run - self.check_fluctuation_dissipation(n) - - # Test case 4.0: no particle specific values / rotational specific global + # Test case 4.0: no particle specific gamma / rotational specific global # thermostat / dissipation only def test_case_40(self): # Each of 2 kind of particles will be represented by n instances: @@ -542,7 +465,7 @@ def test_case_40(self): # Actual integration and validation run self.check_dissipation(n) - # Test case 4.1: no particle specific values / rotational specific global + # Test case 4.1: no particle specific gamma / rotational specific global # thermostat / fluctuation & dissipation def test_case_41(self): # Each of 2 kind of particles will be represented by n instances: diff --git a/testsuite/python/thermostats_common.py b/testsuite/python/thermostats_common.py index ea8c35d1c14..6df58369616 100644 --- a/testsuite/python/thermostats_common.py +++ b/testsuite/python/thermostats_common.py @@ -88,6 +88,15 @@ def check_global(self, N, kT, steps, v_minmax, # Warmup system.integrator.run(20) + vel_obs = ParticleVelocities(ids=system.part[:].id) + vel_acc = TimeSeries(obs=vel_obs) + system.auto_update_accumulators.add(vel_acc) + + if espressomd.has_features("ROTATION"): + omega_obs = ParticleBodyAngularVelocities(ids=system.part[:].id) + omega_acc = TimeSeries(obs=omega_obs) + system.auto_update_accumulators.add(omega_acc) + # Sampling v_stored = np.zeros((steps, N, 3)) omega_stored = np.zeros((steps, N, 3)) @@ -97,14 +106,15 @@ def check_global(self, N, kT, steps, v_minmax, if espressomd.has_features("ROTATION"): omega_stored[i] = system.part[:].omega_body - self.check_velocity_distribution( - v_stored.reshape((-1, 3)), v_minmax, n_bins, error_tol, kT) + vel = vel_acc.time_series().reshape((-1, 3)) + self.check_velocity_distribution(vel, v_minmax, n_bins, error_tol, kT) if espressomd.has_features("ROTATION"): + omega = omega_acc.time_series().reshape((-1, 3)) self.check_velocity_distribution( - omega_stored.reshape((-1, 3)), v_minmax, n_bins, error_tol, kT) + omega, v_minmax, n_bins, error_tol, kT) - def check_per_particle(self, N, kT_global, kT_local, gamma_local, - steps, v_minmax, n_bins, error_tol): + def check_per_particle(self, N, kT, gamma_local, steps, v_minmax, n_bins, + error_tol): """ Test Langevin/Brownian dynamics velocity distribution with global and particle-specific kT/gamma. Covers all combinations of particle @@ -114,10 +124,8 @@ def check_per_particle(self, N, kT_global, kT_local, gamma_local, ---------- N : :obj:`int` Number of particles. - kT_global : :obj:`float` + kT : :obj:`float` Global temperature. - kT_local : :obj:`float` - Per-particle temperature. gamma_local : :obj:`float` Per-particle gamma. steps : :obj:`int` @@ -137,45 +145,28 @@ def check_per_particle(self, N, kT_global, kT_local, gamma_local, if espressomd.has_features("PARTICLE_ANISOTROPY"): gamma_local = 3 * [gamma_local] - # Set different kT on 2nd half of particles - system.part[N // 2:].temp = kT_local - # Set different gamma on half of the particles (overlap over both kTs) - system.part[N // 4:3 * N // 4].gamma = gamma_local + # Set different gamma on 2nd half of particles + system.part[N // 2:].gamma = gamma_local system.integrator.run(50) - vel_obs_global = ParticleVelocities(ids=system.part[:N // 2].id) - vel_obs_local = ParticleVelocities(ids=system.part[N // 2:].id) - vel_acc_global = TimeSeries(obs=vel_obs_global) - vel_acc_local = TimeSeries(obs=vel_obs_local) - system.auto_update_accumulators.add(vel_acc_global) - system.auto_update_accumulators.add(vel_acc_local) + vel_obs = ParticleVelocities(ids=system.part[:].id) + vel_acc = TimeSeries(obs=vel_obs) + system.auto_update_accumulators.add(vel_acc) if espressomd.has_features("ROTATION"): - omega_obs_global = ParticleBodyAngularVelocities( - ids=system.part[:N // 2].id) - omega_obs_local = ParticleBodyAngularVelocities( - ids=system.part[N // 2:].id) - omega_acc_global = TimeSeries(obs=omega_obs_global) - omega_acc_local = TimeSeries(obs=omega_obs_local) - system.auto_update_accumulators.add(omega_acc_global) - system.auto_update_accumulators.add(omega_acc_local) + omega_obs = ParticleBodyAngularVelocities(ids=system.part[:].id) + omega_acc = TimeSeries(obs=omega_obs) + system.auto_update_accumulators.add(omega_acc) system.integrator.run(steps) - vel_global = vel_acc_global.time_series().reshape((-1, 3)) - vel_local = vel_acc_local.time_series().reshape((-1, 3)) - self.check_velocity_distribution( - vel_global, v_minmax, n_bins, error_tol, kT_global) - self.check_velocity_distribution( - vel_local, v_minmax, n_bins, error_tol, kT_local) + vel = vel_acc.time_series().reshape((-1, 3)) + self.check_velocity_distribution(vel, v_minmax, n_bins, error_tol, kT) if espressomd.has_features("ROTATION"): - omega_global = omega_acc_global.time_series().reshape((-1, 3)) - omega_local = omega_acc_local.time_series().reshape((-1, 3)) - self.check_velocity_distribution( - omega_global, v_minmax, n_bins, error_tol, kT_global) + omega = omega_acc.time_series().reshape((-1, 3)) self.check_velocity_distribution( - omega_local, v_minmax, n_bins, error_tol, kT_local) + omega, v_minmax, n_bins, error_tol, kT) def check_noise_correlation(self, kT, steps, delta): """Test the Langevin/Brownian noise is uncorrelated.