diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 79bb7bbf701..7ce6711f3da 100644 --- a/src/core/reaction_ensemble.cpp +++ b/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()); diff --git a/src/core/shapes/SimplePore.cpp b/src/core/shapes/SimplePore.cpp index 796810cf9c3..3fc7f0bd80d 100644 --- a/src/core/shapes/SimplePore.cpp +++ b/src/core/shapes/SimplePore.cpp @@ -34,11 +34,26 @@ std::pair SimplePore::dist_half_pore(double r, double z) const { assert(z >= 0.0); assert(r >= 0.0); - if (z <= c_z) { - /* Cylinder section */ + /* + * We have to find the line that splits area 1 (r determines distance) from + * area 2 (z determines distance) inside pore. In area 3 we have to consider + * z and r to determine the distance. + * + * | x + * | 2 x + * | x + * _|_ x 1 + * \| ^ r + * 3 |-------| + * | z <- + */ + + if ((z <= c_z) && (r <= (c_z + c_r - z))) { + /* Cylinder section, inner */ return {m_rad - r, 0}; - } else if ((r >= c_r) && (z >= m_half_length)) { - /* Wall section */ + } else if (((z >= c_z) && (r >= c_r)) || + ((z <= c_z) && (r > (c_z + c_r - z)))) { + /* Wall section and outer cylinder */ return {0, m_half_length - z}; } else { /* Smoothing area */ @@ -73,11 +88,27 @@ void SimplePore::calculate_dist(const Vector3d &pos, double *dist, double dr, dz; std::tie(dr, dz) = dist_half_pore(r, std::abs(z)); + double side = -1; + if (((dz == 0) && (r <= m_rad)) || // cylinder section + ((dr == 0) && (std::abs(z) > m_half_length))) { // || + side = 1; + } else { + // smoothing area + if (std::abs(z) >= c_z) { + double angle = std::asin((std::abs(z) - c_z) / m_smoothing_rad); + double dist_offset = + m_smoothing_rad - (std::cos(angle) * m_smoothing_rad); + if (m_half_length < std::abs(z) || r <= (m_rad + dist_offset)) { + side = 1; + } + } + } + if (z <= 0.0) { dz *= -1; } - *dist = std::sqrt(dr * dr + dz * dz); + *dist = std::sqrt(dr * dr + dz * dz) * side; for (int i = 0; i < 3; i++) { vec[i] = -dr * e_r[i] + -dz * e_z[i]; } diff --git a/src/python/espressomd/reaction_ensemble.pyx b/src/python/espressomd/reaction_ensemble.pyx index 627e88b313a..639985fec36 100644 --- a/src/python/espressomd/reaction_ensemble.pyx +++ b/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 ---------- @@ -226,6 +230,11 @@ cdef class ReactionAlgorithm(object): "No dictionary for relation between types and default charges provided.") #check electroneutrality of the provided reaction if(self._params["check_for_electroneutrality"]): + charges = np.array(list(self._params["default_charges"].values())) + if(np.count_nonzero(charges) == 0): + # all partices have zero charge + # no need to check electroneutrality + return total_charge_change = 0.0 for i in range(len(self._params["reactant_coefficients"])): type_here = self._params["reactant_types"][i] @@ -235,7 +244,6 @@ cdef class ReactionAlgorithm(object): type_here = self._params["product_types"][j] total_charge_change += self._params["product_coefficients"][ j] * self._params["default_charges"][type_here] - charges = np.array(list(self._params["default_charges"].values())) min_abs_nonzero_charge = np.min( np.abs(charges[np.nonzero(charges)[0]])) if abs(total_charge_change) / min_abs_nonzero_charge > 1e-10: diff --git a/testsuite/python/simple_pore.py b/testsuite/python/simple_pore.py index d33c85e1696..18cf31c05c0 100644 --- a/testsuite/python/simple_pore.py +++ b/testsuite/python/simple_pore.py @@ -31,7 +31,17 @@ "Features not available, skipping test!") class SimplePoreConstraint(ut.TestCase): - def test(self): + def test_orientation(self): + pore = SimplePore( + axis=[1., 0., 0.], radius=2., smoothing_radius=.1, length=2., center=[5., 5., 5.]) + + d, _ = pore.call_method("calc_distance", position=[.0, .0, .0]) + self.assertGreater(d, 0.) + + d, _ = pore.call_method("calc_distance", position=[5., 5., .0]) + self.assertLess(d, 0.) + + def test_stability(self): s = espressomd.System(box_l=[1.0, 1.0, 1.0]) s.seed = s.cell_system.get_state()['n_nodes'] * [1234] box_yz = 15.