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

@@ -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

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``

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

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']

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

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;

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<double, &Prop::T>
 #ifndef PARTICLE_ANISOTROPY
         , UpdateProperty<double, &Prop::gamma>
 #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<double, &ParticleProperties::T>(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

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

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<RNGSalt::BROWNIAN_INC>(
       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<double>
 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<RNGSalt::BROWNIAN_ROT_WALK>(

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
 

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 (

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: