From 314115186f79749bda8bfc3280cb37dfa03fea48 Mon Sep 17 00:00:00 2001 From: Pascal Hebbeker Date: Thu, 15 Nov 2018 15:45:11 +0100 Subject: [PATCH 1/9] use Maxwell-Boltzmann velocity distribution in reaction ensemble --- src/core/reaction_ensemble.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 79bb7bbf701..ecab247e950 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 From a611867bea5c92052701db193a59e7e653069802 Mon Sep 17 00:00:00 2001 From: Pascal Hebbeker Date: Tue, 20 Nov 2018 09:54:32 +0100 Subject: [PATCH 2/9] Support neutral charges in the reaction ensemble Fixes espressomd/espresso#2374 --- src/python/espressomd/reaction_ensemble.pyx | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/python/espressomd/reaction_ensemble.pyx b/src/python/espressomd/reaction_ensemble.pyx index 627e88b313a..e72d57c39af 100644 --- a/src/python/espressomd/reaction_ensemble.pyx +++ b/src/python/espressomd/reaction_ensemble.pyx @@ -226,6 +226,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 +240,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: From 30ce98c9e0289b6394e2b9f5cf654bc7071c9f1b Mon Sep 17 00:00:00 2001 From: Michael Lahnert Date: Tue, 20 Nov 2018 12:09:03 +0100 Subject: [PATCH 3/9] Fix Simple Pore --- src/core/shapes/SimplePore.cpp | 129 ++++++++++++++++++++------------- 1 file changed, 80 insertions(+), 49 deletions(-) diff --git a/src/core/shapes/SimplePore.cpp b/src/core/shapes/SimplePore.cpp index 796810cf9c3..96b5f3ba0b9 100644 --- a/src/core/shapes/SimplePore.cpp +++ b/src/core/shapes/SimplePore.cpp @@ -30,56 +30,87 @@ namespace Shapes { * * @returns The distance vector from the surface in the cylinder system. */ -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 */ - return {m_rad - r, 0}; - } else if ((r >= c_r) && (z >= m_half_length)) { - /* Wall section */ - return {0, m_half_length - z}; - } else { - /* Smoothing area */ - /* Vector to center of torus segment */ - auto const dr = c_r - r; - auto const dz = c_z - z; - - /* Rescale to surface */ - auto const d = std::sqrt(dr * dr + dz * dz); - auto const fac = (d - m_smoothing_rad) / d; - - return {fac * dr, fac * dz}; - } -} - -void SimplePore::calculate_dist(const Vector3d &pos, double *dist, - double *vec) const { - /* Coordinate transform to cylinder coords - with origin at m_center. */ - Vector3d const c_dist = pos - m_center; - auto const z = e_z * c_dist; - auto const r_vec = c_dist - z * e_z; - auto const r = r_vec.norm(); - - /* If exactly on the axis, chose e_r orthogonal - to e_z. */ - auto const e_r = (r == 0) ? e_r_axis : r_vec / r; - - /* The pore has mirror symmetry in z with regard to - the center in the {r,z} system. We calculate always - for the z > 0 case, and flip the result if appropriate. */ - double dr, dz; - std::tie(dr, dz) = dist_half_pore(r, std::abs(z)); - - if (z <= 0.0) { - dz *= -1; + std::pair SimplePore::dist_half_pore(double r, double z) const { + assert(z >= 0.0); + assert(r >= 0.0); + + /* + * 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 (((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 */ + /* Vector to center of torus segment */ + auto const dr = c_r - r; + auto const dz = c_z - z; + + /* Rescale to surface */ + auto const d = std::sqrt(dr * dr + dz * dz); + auto const fac = (d - m_smoothing_rad) / d; + + return {fac * dr, fac * dz}; + } } - *dist = std::sqrt(dr * dr + dz * dz); - for (int i = 0; i < 3; i++) { - vec[i] = -dr * e_r[i] + -dz * e_z[i]; + void SimplePore::calculate_dist(const Vector3d &pos, double *dist, + double *vec) const { + /* Coordinate transform to cylinder coords + with origin at m_center. */ + Vector3d const c_dist = pos - m_center; + auto const z = e_z * c_dist; + auto const r_vec = c_dist - z * e_z; + auto const r = r_vec.norm(); + + /* If exactly on the axis, chose e_r orthogonal + to e_z. */ + auto const e_r = (r == 0) ? e_r_axis : r_vec / r; + + /* The pore has mirror symmetry in z with regard to + the center in the {r,z} system. We calculate always + for the z > 0 case, and flip the result if appropriate. */ + 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) * side; + for (int i = 0; i < 3; i++) { + vec[i] = -dr * e_r[i] + -dz * e_z[i]; + } } -} } // namespace Shapes From 529e827b59b13670b21cab9147e562b058a8cd8a Mon Sep 17 00:00:00 2001 From: Michael Lahnert Date: Tue, 20 Nov 2018 12:23:17 +0100 Subject: [PATCH 4/9] Apply coding style patch --- src/core/shapes/SimplePore.cpp | 152 ++++++++++++++++----------------- 1 file changed, 76 insertions(+), 76 deletions(-) diff --git a/src/core/shapes/SimplePore.cpp b/src/core/shapes/SimplePore.cpp index 96b5f3ba0b9..3fc7f0bd80d 100644 --- a/src/core/shapes/SimplePore.cpp +++ b/src/core/shapes/SimplePore.cpp @@ -30,87 +30,87 @@ namespace Shapes { * * @returns The distance vector from the surface in the cylinder system. */ - std::pair SimplePore::dist_half_pore(double r, double z) const { - assert(z >= 0.0); - assert(r >= 0.0); - - /* - * 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 (((z >= c_z) && (r >= c_r)) || +std::pair SimplePore::dist_half_pore(double r, double z) const { + assert(z >= 0.0); + assert(r >= 0.0); + + /* + * 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 (((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 */ - /* Vector to center of torus segment */ - auto const dr = c_r - r; - auto const dz = c_z - z; - - /* Rescale to surface */ - auto const d = std::sqrt(dr * dr + dz * dz); - auto const fac = (d - m_smoothing_rad) / d; - - return {fac * dr, fac * dz}; - } + /* Wall section and outer cylinder */ + return {0, m_half_length - z}; + } else { + /* Smoothing area */ + /* Vector to center of torus segment */ + auto const dr = c_r - r; + auto const dz = c_z - z; + + /* Rescale to surface */ + auto const d = std::sqrt(dr * dr + dz * dz); + auto const fac = (d - m_smoothing_rad) / d; + + return {fac * dr, fac * dz}; } - - void SimplePore::calculate_dist(const Vector3d &pos, double *dist, - double *vec) const { - /* Coordinate transform to cylinder coords - with origin at m_center. */ - Vector3d const c_dist = pos - m_center; - auto const z = e_z * c_dist; - auto const r_vec = c_dist - z * e_z; - auto const r = r_vec.norm(); - - /* If exactly on the axis, chose e_r orthogonal - to e_z. */ - auto const e_r = (r == 0) ? e_r_axis : r_vec / r; - - /* The pore has mirror symmetry in z with regard to - the center in the {r,z} system. We calculate always - for the z > 0 case, and flip the result if appropriate. */ - 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; - } +} + +void SimplePore::calculate_dist(const Vector3d &pos, double *dist, + double *vec) const { + /* Coordinate transform to cylinder coords + with origin at m_center. */ + Vector3d const c_dist = pos - m_center; + auto const z = e_z * c_dist; + auto const r_vec = c_dist - z * e_z; + auto const r = r_vec.norm(); + + /* If exactly on the axis, chose e_r orthogonal + to e_z. */ + auto const e_r = (r == 0) ? e_r_axis : r_vec / r; + + /* The pore has mirror symmetry in z with regard to + the center in the {r,z} system. We calculate always + for the z > 0 case, and flip the result if appropriate. */ + 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; - } + if (z <= 0.0) { + dz *= -1; + } - *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]; - } + *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]; } +} } // namespace Shapes From 2f933fef6e929f7973f10e685ff899f8461b6a70 Mon Sep 17 00:00:00 2001 From: Pascal Hebbeker Date: Mon, 26 Nov 2018 14:41:11 +0100 Subject: [PATCH 5/9] note the missing mass dependency of RE in the doc --- src/python/espressomd/reaction_ensemble.pyx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/python/espressomd/reaction_ensemble.pyx b/src/python/espressomd/reaction_ensemble.pyx index 627e88b313a..3d30154e7de 100644 --- a/src/python/espressomd/reaction_ensemble.pyx +++ b/src/python/espressomd/reaction_ensemble.pyx @@ -16,6 +16,9 @@ cdef class ReactionAlgorithm(object): reaction algorithm by setting the standard pressure, temperature, and the exclusion radius. + Note: When setting the velocities of the particles according the the + Maxwell distribution the mass of all particles is assumed to equal 1. + Parameters ---------- From a28b747f62c54785acabb0f127a1e25e6c416d91 Mon Sep 17 00:00:00 2001 From: Pascal Hebbeker Date: Mon, 26 Nov 2018 18:18:10 +0100 Subject: [PATCH 6/9] set particle velocities in the RE MC move Draw a random velocity according to the Maxwell-Boltzmann distribution. Fixes #1770 --- src/core/reaction_ensemble.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index ecab247e950..01874ab1e76 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -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]; + // we use mass=1 for all particles, think about adapting this + vel[0] = std::sqrt(temperature) * gaussian_random(); + vel[1] = std::sqrt(temperature) * gaussian_random(); + vel[2] = std::sqrt(temperature) * 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()); From 670a4df64d872d7b1e992388b6f17bf4fce44f9b Mon Sep 17 00:00:00 2001 From: Pascal Hebbeker Date: Tue, 27 Nov 2018 16:14:43 +0100 Subject: [PATCH 7/9] account for the particle mass in the RE MC move --- src/core/reaction_ensemble.cpp | 8 ++++---- src/python/espressomd/reaction_ensemble.pyx | 5 +++-- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 01874ab1e76..7ce6711f3da 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -832,10 +832,10 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type( // change particle position new_pos = get_random_position_in_box(); double vel[3]; - // we use mass=1 for all particles, think about adapting this - vel[0] = std::sqrt(temperature) * gaussian_random(); - vel[1] = std::sqrt(temperature) * gaussian_random(); - vel[2] = std::sqrt(temperature) * gaussian_random(); + 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 diff --git a/src/python/espressomd/reaction_ensemble.pyx b/src/python/espressomd/reaction_ensemble.pyx index 3d30154e7de..dbe44529cfe 100644 --- a/src/python/espressomd/reaction_ensemble.pyx +++ b/src/python/espressomd/reaction_ensemble.pyx @@ -16,8 +16,9 @@ cdef class ReactionAlgorithm(object): reaction algorithm by setting the standard pressure, temperature, and the exclusion radius. - Note: When setting the velocities of the particles according the the - Maxwell distribution the mass of all particles is assumed to equal 1. + Note: When creating particles the velocities the new particles are set + according the Maxwell distribution. In this step the mass of the new particle + is assumed to equal 1. Parameters From 33bf7840b69b5cfafea5567b436b243ac24d3e47 Mon Sep 17 00:00:00 2001 From: Pascal Hebbeker Date: Fri, 30 Nov 2018 11:30:17 +0100 Subject: [PATCH 8/9] fix spelling of the documentation --- src/python/espressomd/reaction_ensemble.pyx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/python/espressomd/reaction_ensemble.pyx b/src/python/espressomd/reaction_ensemble.pyx index dbe44529cfe..0bfcd655275 100644 --- a/src/python/espressomd/reaction_ensemble.pyx +++ b/src/python/espressomd/reaction_ensemble.pyx @@ -16,9 +16,9 @@ cdef class ReactionAlgorithm(object): reaction algorithm by setting the standard pressure, temperature, and the exclusion radius. - Note: When creating particles the velocities the new particles are set - according the Maxwell distribution. In this step the mass of the new particle - is assumed to equal 1. + 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 From d48fc009e7f71e4c19ae1465a642b59e135f446d Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Fri, 30 Nov 2018 17:21:36 +0100 Subject: [PATCH 9/9] testsuite: Added test for orientation of SimplePore. --- testsuite/python/simple_pore.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) 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.