From e1fd461c429b84ec99753a92498a3b750d3f4c74 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Thu, 12 Dec 2019 14:24:32 +0100 Subject: [PATCH] Fix particle-id based observables with non-contigeous particle numbering --- src/core/observables/ComForce.hpp | 6 +- .../observables/ParticleAngularVelocities.hpp | 3 +- .../observables/ParticleBodyVelocities.hpp | 4 +- src/python/espressomd/observables.py | 2 +- testsuite/python/observables.py | 111 ++++++++++++------ 5 files changed, 81 insertions(+), 45 deletions(-) diff --git a/src/core/observables/ComForce.hpp b/src/core/observables/ComForce.hpp index 59aa56112c6..727f6df74cb 100644 --- a/src/core/observables/ComForce.hpp +++ b/src/core/observables/ComForce.hpp @@ -34,9 +34,9 @@ class ComForce : public PidObservable { for (int i : ids()) { if (partCfg[i].p.is_virtual) continue; - res[0] += partCfg[i].f.f[0] * partCfg[i].p.mass; - res[1] += partCfg[i].f.f[1] * partCfg[i].p.mass; - res[2] += partCfg[i].f.f[2] * partCfg[i].p.mass; + res[0] += partCfg[i].f.f[0]; + res[1] += partCfg[i].f.f[1]; + res[2] += partCfg[i].f.f[2]; } return res; }; diff --git a/src/core/observables/ParticleAngularVelocities.hpp b/src/core/observables/ParticleAngularVelocities.hpp index 89072247484..9b47aa64a46 100644 --- a/src/core/observables/ParticleAngularVelocities.hpp +++ b/src/core/observables/ParticleAngularVelocities.hpp @@ -34,8 +34,9 @@ class ParticleAngularVelocities : public PidObservable { std::vector res(n_values()); for (int i = 0; i < ids().size(); i++) { #ifdef ROTATION + auto const id = ids()[i]; const Utils::Vector3d omega = - convert_vector_body_to_space(partCfg[i], partCfg[i].m.omega); + convert_vector_body_to_space(partCfg[id], partCfg[id].m.omega); res[3 * i + 0] = omega[0]; res[3 * i + 1] = omega[1]; res[3 * i + 2] = omega[2]; diff --git a/src/core/observables/ParticleBodyVelocities.hpp b/src/core/observables/ParticleBodyVelocities.hpp index 5ec2d3b05e8..c46e61a40e2 100644 --- a/src/core/observables/ParticleBodyVelocities.hpp +++ b/src/core/observables/ParticleBodyVelocities.hpp @@ -35,9 +35,9 @@ class ParticleBodyVelocities : public PidObservable { for (int i = 0; i < ids().size(); i++) { #ifdef ROTATION - double RMat[9]; + auto const id = ids()[i]; const Utils::Vector3d vel_body = - convert_vector_space_to_body(partCfg[i], partCfg[i].m.v); + convert_vector_space_to_body(partCfg[id], partCfg[id].m.v); res[3 * i + 0] = vel_body[0]; res[3 * i + 1] = vel_body[1]; diff --git a/src/python/espressomd/observables.py b/src/python/espressomd/observables.py index b5ef982797b..1404aefe35c 100644 --- a/src/python/espressomd/observables.py +++ b/src/python/espressomd/observables.py @@ -29,7 +29,7 @@ class ComForce(Observable): """Calculates the total force on particles with given ids. - Note that virtual sites are not included since they do not have a meaningful mass. + Note that virtual sites are not included since forces on them do not enter the equation of motion directly. Output format: :math:`\\left(\\sum_i f^x_i, \\sum_i f^y_i, \\sum_i f^z_i\\right)` diff --git a/testsuite/python/observables.py b/testsuite/python/observables.py index 059dab180b8..40bbcb8d55f 100644 --- a/testsuite/python/observables.py +++ b/testsuite/python/observables.py @@ -19,45 +19,54 @@ import unittest_decorators as utx import espressomd import numpy as np -from numpy.random import random, randint import espressomd.observables -def calc_com_x(system, x): - masses = system.part[:].mass +def calc_com_x(system, x, id_list): + """Mass-weighted average, skipping virtual sites""" + masses = system.part[id_list].mass - # Virtual sites are excluded since they do not have meaningful mass - if espressomd.has_features("VIRTUAL_SITES"): - for i, p in enumerate(system.part): - if p.virtual: - masses[i] = 0. + # Filter out virtual particles by using mass=0 for them + virtual = system.part[id_list].virtual + for i in range(len(masses)): + if virtual[i]: + masses[i] = 0. com_x = np.average( - getattr(system.part[:], x), weights=masses, axis=0) + getattr(system.part[id_list], x), weights=masses, axis=0) return com_x class Observables(ut.TestCase): - N_PART = 1000 + N_PART = 200 # Handle for espresso system system = espressomd.System(box_l=[10.0, 10.0, 10.0]) - system.seed = system.cell_system.get_state()['n_nodes'] * [1234] - - def setUp(self): - if not self.system.part: - for i in range(self.N_PART): - self.system.part.add(pos=random(3) * 10, v=random(3), id=i) - if espressomd.has_features(["MASS"]): - self.system.part[i].mass = random() - if espressomd.has_features(["DIPOLES"]): - self.system.part[i].dip = random(3) - if espressomd.has_features(["ROTATION"]): - self.system.part[i].omega_lab = random(3) - if espressomd.has_features("ELECTROSTATICS"): - self.system.part[i].q = (1 if i % 2 == 0 else -1) - - if espressomd.has_features("VIRTUAL_SITES"): - self.system.part[randint(self.N_PART)].virtual = True + system.part.add( + id=np.arange(3, 3 + 2 * N_PART, 2), + pos=np.random.random((N_PART, 3)) * system.box_l, + v=np.random.random((N_PART, 3)) * 3.2 - 1, + f=np.random.random((N_PART, 3))) + + if espressomd.has_features(["MASS"]): + system.part[:].mass = np.random.random(N_PART) + + if espressomd.has_features(["DIPOLES"]): + system.part[:].dip = np.random.random((N_PART, 3)) - .3 + + if espressomd.has_features(["ROTATION"]): + system.part[:].omega_body = np.random.random((N_PART, 3)) - .5 + system.part[:].torque_lab = np.random.random((N_PART, 3)) - .5 + system.part[:].quat = np.random.random((N_PART, 4)) + + if espressomd.has_features("DIPOLES"): + system.part[:].dipm = np.random.random(N_PART) + 2 + + if espressomd.has_features("ELECTROSTATICS"): + system.part[:].q = np.random.random(N_PART) + + if espressomd.has_features("VIRTUAL_SITES"): + p = system.part[system.part[:].id[8]] + p.virtual = True def generate_test_for_pid_observable( _obs_name, _pprop_name, _agg_type=None): @@ -72,9 +81,19 @@ def func(self): # This code is run at the execution of the generated function. # It will use the state of the variables in the outer function, # which was there, when the outer function was called + # Randomly pick a subset of the particles + id_list = sorted( + np.random.choice( + self.system.part[:].id, + size=int( + self.N_PART * .9), + replace=False)) + for id in id_list: + assert(self.system.part.exists(id)) + # Get data from particles - id_list = range(self.N_PART) part_data = getattr(self.system.part[id_list], pprop_name) + # Reshape and aggregate to linear array if len(part_data.shape) > 1: if agg_type == "average": @@ -82,19 +101,22 @@ def func(self): if agg_type == "sum": part_data = np.sum(part_data, 0) if agg_type == 'com': - part_data = calc_com_x(self.system, pprop_name) - part_data = part_data.flatten() + part_data = calc_com_x(self.system, pprop_name, id_list) # Data from observable observable = obs_name(ids=id_list) - obs_data = observable.calculate() - self.assertEqual(observable.n_values(), len(part_data)) + obs_data = np.reshape(observable.calculate(), part_data.shape) + + # Check + self.assertEqual(observable.n_values(), len(part_data.flatten())) + np.testing.assert_equal(id_list, observable.ids) + np.testing.assert_array_almost_equal( obs_data, part_data, err_msg="Data did not agree for observable " + str(obs_name) + " and particle property " + - pprop_name, decimal=9) + pprop_name, decimal=11) return func @@ -108,8 +130,6 @@ def func(self): espressomd.observables.ComPosition, 'pos', 'com') test_com_velocity = generate_test_for_pid_observable( espressomd.observables.ComVelocity, 'v', 'com') - test_com_force = generate_test_for_pid_observable( - espressomd.observables.ComForce, 'f', 'com') if espressomd.has_features(["DIPOLES"]): test_mag_dip = generate_test_for_pid_observable( @@ -124,7 +144,7 @@ def func(self): @utx.skipIfMissingFeatures(['ROTATION']) def test_particle_body_velocities(self): obs = espressomd.observables.ParticleBodyVelocities( - ids=range(self.N_PART)) + ids=self.system.part[:].id) obs_data = obs.calculate() part_data = np.array([p.convert_vector_space_to_body(p.v) for p in self.system.part]) @@ -146,7 +166,7 @@ def test_stress_tensor(self): @utx.skipIfMissingFeatures('ELECTROSTATICS') def test_current(self): obs_data = espressomd.observables.Current( - ids=range(self.N_PART)).calculate() + ids=self.system.part[:].id).calculate() part_data = self.system.part[:].q.dot(self.system.part[:].v) self.assertEqual(espressomd.observables.Current( ids=range(self.N_PART)).n_values(), len(part_data.flatten())) @@ -155,13 +175,28 @@ def test_current(self): @utx.skipIfMissingFeatures('ELECTROSTATICS') def test_dipolemoment(self): - obs = espressomd.observables.DipoleMoment(ids=range(self.N_PART)) + obs = espressomd.observables.DipoleMoment(ids=self.system.part[:].id) obs_data = obs.calculate() part_data = self.system.part[:].q.dot(self.system.part[:].pos) self.assertEqual(obs.n_values(), len(part_data.flatten())) np.testing.assert_array_almost_equal( obs_data, part_data, err_msg="Data did not agree for observable 'DipoleMoment'", decimal=9) + def test_com_force(self): + id_list = sorted( + np.random.choice( + self.system.part[:].id, + size=int( + self.N_PART * .9), + replace=False)) + + particles = self.system.part.select( + lambda p: p.id in id_list and not p.virtual) + + np.testing.assert_allclose( + np.sum(particles.f, axis=0), + espressomd.observables.ComForce(ids=id_list).calculate()) + if __name__ == "__main__": ut.main()