Merge pull request #2377 from thepith/maxwell_re

use correct Maxwell-Boltzmann velocity distribution in reaction ensemble

src/core/reaction_ensemble.cpp

@@ -717,9 +717,9 @@ int ReactionAlgorithm::create_particle(int desired_type) {
   // for components
   double vel[3];
   // we use mass=1 for all particles, think about adapting this
-  vel[0] = std::pow(2 * PI * temperature, -3.0 / 2.0) * gaussian_random();
-  vel[1] = std::pow(2 * PI * temperature, -3.0 / 2.0) * gaussian_random();
-  vel[2] = std::pow(2 * PI * temperature, -3.0 / 2.0) * gaussian_random();
+  vel[0] = std::sqrt(temperature) * gaussian_random();
+  vel[1] = std::sqrt(temperature) * gaussian_random();
+  vel[2] = std::sqrt(temperature) * gaussian_random();
 #ifdef ELECTROSTATICS
   double charge = charges_of_types[desired_type];
 #endif
@@ -831,6 +831,12 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type(
     p_id = p_id_s_changed_particles[i];
     // change particle position
     new_pos = get_random_position_in_box();
+    double vel[3];
+    auto const &p = get_particle_data(p_id);
+    vel[0] = std::sqrt(temperature / p.p.mass) * gaussian_random();
+    vel[1] = std::sqrt(temperature / p.p.mass) * gaussian_random();
+    vel[2] = std::sqrt(temperature / p.p.mass) * gaussian_random();
+    set_particle_v(p_id, vel);
     // new_pos=get_random_position_in_box_enhanced_proposal_of_small_radii();
     // //enhanced proposal of small radii
     place_particle(p_id, new_pos.data());

src/python/espressomd/reaction_ensemble.pyx

@@ -16,6 +16,10 @@ cdef class ReactionAlgorithm(object):
     reaction algorithm by setting the standard pressure, temperature, and the
     exclusion radius.
 
+    Note: When creating particles the velocities of the new particles are set
+    according the Maxwell-Boltzmann distribution. In this step the mass of the
+    new particle is assumed to equal 1.
+
 
     Parameters
     ----------