From b6aa78eb4f7164b87e2dc7975156a7984f71da94 Mon Sep 17 00:00:00 2001 From: Christoph Lohrmann Date: Tue, 7 Dec 2021 18:02:03 +0100 Subject: [PATCH 1/2] rework particle access --- .../charged_system/charged_system.ipynb | 8 +- .../ferrofluid/ferrofluid_part1.ipynb | 6 +- .../ferrofluid/ferrofluid_part2.ipynb | 14 +- .../ferrofluid/ferrofluid_part3.ipynb | 6 +- .../lennard_jones/lennard_jones.ipynb | 39 ++- .../raspberry_electrophoresis.ipynb | 8 +- samples/chamber_game.py | 12 +- samples/dancing.py | 2 +- samples/electrophoresis.py | 11 +- samples/gibbs_ensemble/gibbs_ensemble.py | 16 +- samples/immersed_boundary/addBending.py | 2 +- samples/immersed_boundary/addSoft.py | 2 +- samples/immersed_boundary/addVolCons.py | 4 +- samples/immersed_boundary/writeVTK.py | 4 +- samples/lj-demo.py | 7 +- samples/load_checkpoint.py | 2 +- samples/rigid_body.py | 8 +- samples/slice_input.py | 43 +-- samples/visualization_elc.py | 5 +- samples/visualization_interactive.py | 9 +- src/python/espressomd/drude_helpers.py | 4 +- src/python/espressomd/io/writer/vtf.py | 12 +- src/python/espressomd/particle_data.pyx | 44 ++- src/python/espressomd/polymer.pyx | 6 +- src/python/object_in_fluid/oif_classes.py | 2 +- testsuite/python/accumulator_mean_variance.py | 2 +- testsuite/python/accumulator_time_series.py | 2 +- testsuite/python/analyze_chains.py | 55 +-- testsuite/python/analyze_distance.py | 31 +- testsuite/python/analyze_distribution.py | 6 +- testsuite/python/analyze_energy.py | 10 +- testsuite/python/analyze_gyration_tensor.py | 2 +- testsuite/python/analyze_mass_related.py | 11 +- testsuite/python/array_properties.py | 98 +++--- testsuite/python/auto_exclusions.py | 12 +- testsuite/python/cluster_analysis.py | 6 +- testsuite/python/collision_detection.py | 61 ++-- testsuite/python/comfixed.py | 7 +- testsuite/python/constraint_shape_based.py | 2 +- testsuite/python/coulomb_cloud_wall.py | 2 +- .../python/coulomb_cloud_wall_duplicated.py | 2 +- testsuite/python/coulomb_mixed_periodicity.py | 2 +- testsuite/python/coulomb_tuning.py | 2 +- testsuite/python/dawaanr-and-bh-gpu.py | 8 +- testsuite/python/dawaanr-and-dds-gpu.py | 8 +- testsuite/python/dds-and-bh-gpu.py | 12 +- testsuite/python/dipolar_direct_summation.py | 27 +- .../dipolar_mdlc_p3m_scafacos_p2nfft.py | 29 +- testsuite/python/dipolar_p3m.py | 21 +- testsuite/python/domain_decomposition.py | 9 +- testsuite/python/dpd.py | 20 +- testsuite/python/dpd_stats.py | 18 +- testsuite/python/elc_vs_analytic.py | 2 +- testsuite/python/electrostaticInteractions.py | 36 +- testsuite/python/engine_langevin.py | 2 +- testsuite/python/exclusions.py | 80 ++--- testsuite/python/force_cap.py | 2 +- testsuite/python/galilei.py | 24 +- testsuite/python/h5md.py | 12 +- testsuite/python/ibm.py | 21 +- testsuite/python/integrator_npt.py | 10 +- .../python/integrator_steepest_descent.py | 64 ++-- testsuite/python/interactions_bond_angle.py | 37 +- testsuite/python/interactions_bonded.py | 6 +- testsuite/python/interactions_dihedral.py | 4 +- testsuite/python/interactions_non-bonded.py | 4 +- testsuite/python/langevin_thermostat_stats.py | 4 +- testsuite/python/lb_get_u_at_pos.py | 2 +- testsuite/python/lb_stats.py | 4 +- testsuite/python/lb_thermostat.py | 2 +- testsuite/python/lj.py | 5 +- .../python/mass-and-rinertia_per_particle.py | 32 +- testsuite/python/mdanalysis.py | 15 +- testsuite/python/mmm1d.py | 2 +- testsuite/python/mpiio.py | 34 +- testsuite/python/nsquare.py | 7 +- testsuite/python/observable_chain.py | 24 +- testsuite/python/observable_cylindrical.py | 4 +- testsuite/python/observable_cylindricalLB.py | 2 +- testsuite/python/observable_profile.py | 8 +- testsuite/python/observables.py | 39 ++- testsuite/python/oif_volume_conservation.py | 9 +- testsuite/python/p3m_tuning_exceptions.py | 4 +- testsuite/python/pair_criteria.py | 8 +- testsuite/python/particle.py | 103 +++--- testsuite/python/particle_slice.py | 331 +++++++++--------- testsuite/python/polymer_diamond.py | 7 +- testsuite/python/pressure.py | 5 +- testsuite/python/random_pairs.py | 3 +- testsuite/python/rdf.py | 21 +- testsuite/python/rescale.py | 14 +- .../python/rotational-diffusion-aniso.py | 8 +- testsuite/python/scafacos_dipoles_1d_2d.py | 6 +- testsuite/python/scafacos_interface.py | 10 +- testsuite/python/tabulated.py | 6 +- testsuite/python/test_checkpoint.py | 23 +- testsuite/python/tests_common.py | 4 +- testsuite/python/thermalized_bond.py | 39 ++- testsuite/python/thermostats_common.py | 29 +- testsuite/python/thole.py | 14 +- testsuite/python/virtual_sites_relative.py | 42 +-- testsuite/python/writevtf.py | 14 +- .../samples/test_MDAnalysisIntegration.py | 9 +- .../scripts/samples/test_visualization_elc.py | 3 +- 104 files changed, 1009 insertions(+), 932 deletions(-) diff --git a/doc/tutorials/charged_system/charged_system.ipynb b/doc/tutorials/charged_system/charged_system.ipynb index ed56909fef7..460427eda7f 100644 --- a/doc/tutorials/charged_system/charged_system.ipynb +++ b/doc/tutorials/charged_system/charged_system.ipynb @@ -223,7 +223,7 @@ " ROD_CHARGE_DENS, N_rod_beads, ROD_TYPE)\n", "\n", "# check that the particle setup was done correctly\n", - "assert abs(sum(system.part[:].q)) < 1e-10\n", + "assert abs(sum(system.part.all().q)) < 1e-10\n", "assert np.all(system.part.select(type=ROD_TYPE).fix)" ] }, @@ -303,7 +303,7 @@ "\n", " # Initialize integrator to obtain initial forces\n", " system.integrator.run(0)\n", - " maxforce = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", + " maxforce = np.max(np.linalg.norm(system.part.all().f, axis=1))\n", " energy = system.analysis.energy()['total']\n", "\n", " i = 0\n", @@ -311,7 +311,7 @@ " prev_maxforce = maxforce\n", " prev_energy = energy\n", " system.integrator.run(sd_params['emstep'])\n", - " maxforce = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", + " maxforce = np.max(np.linalg.norm(system.part.all().f, axis=1))\n", " relforce = np.abs((maxforce - prev_maxforce) / prev_maxforce)\n", " energy = system.analysis.energy()['total']\n", " relener = np.abs((energy - prev_energy) / prev_energy)\n", @@ -944,7 +944,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.5" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/doc/tutorials/ferrofluid/ferrofluid_part1.ipynb b/doc/tutorials/ferrofluid/ferrofluid_part1.ipynb index 8779d3f05ba..6a2f2370cd4 100644 --- a/doc/tutorials/ferrofluid/ferrofluid_part1.ipynb +++ b/doc/tutorials/ferrofluid/ferrofluid_part1.ipynb @@ -714,7 +714,7 @@ " system.integrator.run(100)\n", "\n", " # Save current system state as a plot\n", - " x_data, y_data = system.part[:].pos_folded[:, 0], system.part[:].pos_folded[:, 1]\n", + " x_data, y_data = system.part.all().pos_folded[:, 0], system.part.all().pos_folded[:, 1]\n", " ax.figure.canvas.draw()\n", " part.set_data(x_data, y_data)\n", " print(\"progress: {:3.0f}%\".format((i + 1) * 100. / LOOPS), end=\"\\r\")\n", @@ -905,7 +905,7 @@ "plt.ylim(0, BOX_SIZE)\n", "plt.xlabel('x-position', fontsize=20)\n", "plt.ylabel('y-position', fontsize=20)\n", - "plt.plot(system.part[:].pos_folded[:, 0], system.part[:].pos_folded[:, 1], 'o')\n", + "plt.plot(system.part.all().pos_folded[:, 0], system.part.all().pos_folded[:, 1], 'o')\n", "plt.show()" ] }, @@ -1026,7 +1026,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.9" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/doc/tutorials/ferrofluid/ferrofluid_part2.ipynb b/doc/tutorials/ferrofluid/ferrofluid_part2.ipynb index eee066cd979..ce3266f047f 100644 --- a/doc/tutorials/ferrofluid/ferrofluid_part2.ipynb +++ b/doc/tutorials/ferrofluid/ferrofluid_part2.ipynb @@ -152,7 +152,7 @@ "pos = box_size * np.hstack((np.random.random((N_PART, 2)), np.zeros((N_PART, 1))))\n", "\n", "# Add particles\n", - "system.part.add(pos=pos, rotation=N_PART * [(1, 1, 1)], dip=dip, fix=N_PART * [(0, 0, 1)])\n", + "particles = system.part.add(pos=pos, rotation=N_PART * [(True, True, True)], dip=dip, fix=N_PART * [(False, False, True)])\n", "\n", "# Remove overlap between particles by means of the steepest descent method\n", "system.integrator.set_steepest_descent(\n", @@ -289,7 +289,7 @@ "plt.ylim(0, box_size)\n", "plt.xlabel('x-position', fontsize=20)\n", "plt.ylabel('y-position', fontsize=20)\n", - "plt.plot(system.part[:].pos_folded[:, 0], system.part[:].pos_folded[:, 1], 'o')\n", + "plt.plot(particles.pos_folded[:, 0], particles.pos_folded[:, 1], 'o')\n", "plt.show()" ] }, @@ -351,7 +351,7 @@ " system.integrator.run(50)\n", "\n", " # Save current system state as a plot\n", - " xdata, ydata = system.part[:].pos_folded[:, 0], system.part[:].pos_folded[:, 1]\n", + " xdata, ydata = particles.pos_folded[:, 0], particles.pos_folded[:, 1]\n", " ax.figure.canvas.draw()\n", " part.set_data(xdata, ydata)\n", " print(\"progress: {:3.0f}%\".format(i + 1), end=\"\\r\")\n", @@ -498,7 +498,7 @@ "outputs": [], "source": [ "# remove all particles\n", - "system.part[:].remove()\n", + "system.part.clear()\n", "system.thermostat.turn_off()\n", "\n", "# Random dipole moments\n", @@ -514,7 +514,7 @@ "pos = box_size * np.hstack((np.random.random((N_PART, 2)), np.zeros((N_PART, 1))))\n", "\n", "# Add particles\n", - "system.part.add(pos=pos, rotation=N_PART * [(1, 1, 1)], dip=dip, fix=N_PART * [(0, 0, 1)])\n", + "particles = system.part.add(pos=pos, rotation=N_PART * [(True, True, True)], dip=dip, fix=N_PART * [(False, False, True)])\n", "\n", "# Remove overlap between particles by means of the steepest descent method\n", "system.integrator.set_steepest_descent(f_max=0, gamma=0.1, max_displacement=0.05)\n", @@ -570,7 +570,7 @@ "source": [ "```python\n", "import espressomd.observables\n", - "dipm_tot = espressomd.observables.MagneticDipoleMoment(ids=system.part[:].id)\n", + "dipm_tot = espressomd.observables.MagneticDipoleMoment(ids=particles.id)\n", "```" ] }, @@ -1041,7 +1041,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.9" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/doc/tutorials/ferrofluid/ferrofluid_part3.ipynb b/doc/tutorials/ferrofluid/ferrofluid_part3.ipynb index cf653cf82e4..37c822ec07b 100644 --- a/doc/tutorials/ferrofluid/ferrofluid_part3.ipynb +++ b/doc/tutorials/ferrofluid/ferrofluid_part3.ipynb @@ -215,7 +215,7 @@ "pos = box_size * np.random.random((N, 3))\n", "\n", "# Add particles\n", - "system.part.add(pos=pos, rotation=N * [(1, 1, 1)], dip=dip)\n", + "particles = system.part.add(pos=pos, rotation=N * [(True, True, True)], dip=dip)\n", "\n", "# Remove overlap between particles by means of the steepest descent method\n", "system.integrator.set_steepest_descent(\n", @@ -275,7 +275,7 @@ "outputs": [], "source": [ "import espressomd.observables\n", - "dipm_tot_calc = espressomd.observables.MagneticDipoleMoment(ids=system.part[:].id)" + "dipm_tot_calc = espressomd.observables.MagneticDipoleMoment(ids=particles.id)" ] }, { @@ -670,7 +670,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.8" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/doc/tutorials/lennard_jones/lennard_jones.ipynb b/doc/tutorials/lennard_jones/lennard_jones.ipynb index 82e4f4ce0e3..53193d98aad 100644 --- a/doc/tutorials/lennard_jones/lennard_jones.ipynb +++ b/doc/tutorials/lennard_jones/lennard_jones.ipynb @@ -311,7 +311,7 @@ "\n", "### Placing and accessing particles\n", "\n", - "Particles in the simulation can be added and accessed via the part property of the System class. Individual particles are referred to by an integer id, e.g., system.part[0]. If id is unspecified, an unused particle id is automatically assigned. It is also possible to use common python iterators and slicing operations to add or access several particles at once.\n", + "Particles in the simulation can be added and accessed via the part property of the System class. Individual particles should be referred to by the particle handle returned upon creation. You can also retrieve them by an integer id, e.g. system.part.by_id(0). If id is unspecified when creating a particle, an unused particle id is automatically assigned. It is also possible to use common python iterators and slicing operations to add or access several particles at once.\n", "\n", "Particles can be grouped into several types, so that, e.g., a binary fluid can be simulated. Particle types are identified by integer ids, which are set via the particles' type attribute. If it is not specified, zero is implied." ] @@ -325,11 +325,11 @@ "source": [ "\n", "\n", - "* Create N_PART particles at random positions.\n", + "* Create N_PART particles at random positions, store their handles in a variable called ``particles``.\n", "\n", " Use [system.part.add()](https://espressomd.github.io/doc/espressomd.html#espressomd.particle_data.ParticleList).\n", " \n", - " Either write a loop or use an (N_PART x 3) array for positions.\n", + " Use an (N_PART x 3) array for positions.\n", " Use np.random.random() to generate random numbers." ] }, @@ -340,7 +340,7 @@ }, "source": [ "```python\n", - "system.part.add(type=[0] * N_PART, pos=np.random.random((N_PART, 3)) * system.box_l)\n", + "particles = system.part.add(type=[0] * N_PART, pos=np.random.random((N_PART, 3)) * system.box_l)\n", "```" ] }, @@ -375,16 +375,23 @@ "outputs": [], "source": [ "# Access position of a single particle\n", - "print(\"position of particle with id 0:\", system.part[0].pos)\n", + "print(\"position of particle with id 0:\", system.part.by_id(0).pos)\n", "\n", "# Iterate over the first five particles for the purpose of demonstration.\n", - "# For accessing all particles, use a slice: system.part[:]\n", - "for i in range(5):\n", - " print(\"id\", i, \"position:\", system.part[i].pos)\n", - " print(\"id\", i, \"velocity:\", system.part[i].v)\n", + "first_five = system.part.by_ids(range(5))\n", + "for p in first_five:\n", + " print(\"id\", p.id, \"position:\", p.pos)\n", + " print(\"id\", p.id, \"velocity:\", p.v)\n", "\n", - "# Obtain all particle positions\n", - "cur_pos = system.part[:].pos" + "# Obtain particle positions for the particles created until now\n", + "cur_pos = particles.pos" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can also get all particles using ``system.part.all()``, but ``particles`` already contains all particles that are in the simulation so far." ] }, { @@ -400,7 +407,7 @@ "metadata": {}, "outputs": [], "source": [ - "print(system.part[0])" + "print(system.part.by_id(0))" ] }, { @@ -540,12 +547,12 @@ "\n", "# Initialize integrator to obtain initial forces\n", "system.integrator.run(0)\n", - "old_force = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", + "old_force = np.max(np.linalg.norm(system.part.all().f, axis=1))\n", "\n", "\n", "while system.time / system.time_step < MAX_STEPS:\n", " system.integrator.run(EM_STEP)\n", - " force = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", + " force = np.max(np.linalg.norm(system.part.all().f, axis=1))\n", " rel_force = np.abs((force - old_force) / old_force)\n", " print(f'rel. force change: {rel_force:.2e}')\n", " if rel_force < F_TOL:\n", @@ -983,7 +990,7 @@ }, "source": [ "```python\n", - "rdf_obs = espressomd.observables.RDF(ids1=system.part[:].id, min_r=R_MIN, max_r=R_MAX, n_r_bins=N_BINS)\n", + "rdf_obs = espressomd.observables.RDF(ids1=system.part.all().id, min_r=R_MIN, max_r=R_MAX, n_r_bins=N_BINS)\n", "rdf_acc = espressomd.accumulators.MeanVarianceCalculator(obs=rdf_obs, delta_N=steps_per_subsample)\n", "system.auto_update_accumulators.add(rdf_acc)\n", "```" @@ -1211,7 +1218,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.5" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/doc/tutorials/raspberry_electrophoresis/raspberry_electrophoresis.ipynb b/doc/tutorials/raspberry_electrophoresis/raspberry_electrophoresis.ipynb index f9b56c9c3a6..4d7c313240d 100644 --- a/doc/tutorials/raspberry_electrophoresis/raspberry_electrophoresis.ipynb +++ b/doc/tutorials/raspberry_electrophoresis/raspberry_electrophoresis.ipynb @@ -283,7 +283,7 @@ "for _ in range(100):\n", " system.integrator.run(50)\n", " constrain_surface_particles()\n", - " force_max = np.max(np.linalg.norm(system.part[:].f, axis=1))\n", + " force_max = np.max(np.linalg.norm(system.part.all().f, axis=1))\n", " logging.info(f\"maximal force: {force_max:.1f}\")\n", " if force_max < 10.:\n", " break\n", @@ -385,7 +385,7 @@ "source": [ "```python\n", "# Calculate the center of mass position (com) and the moment of inertia (momI) of the colloid\n", - "com = np.average(surface_parts.pos, 0) # system.part[:].pos returns an n-by-3 array\n", + "com = np.average(surface_parts.pos, 0) # surface_parts.pos returns an n-by-3 array\n", "momI = 0\n", "for p in surface_parts:\n", " momI += np.power(np.linalg.norm(com - p.pos), 2)\n", @@ -474,7 +474,7 @@ "outputs": [], "source": [ "# Check charge neutrality\n", - "assert np.abs(np.sum(system.part[:].q)) < 1E-10" + "assert np.abs(np.sum(system.part.all().q)) < 1E-10" ] }, { @@ -620,7 +620,7 @@ "metadata": {}, "outputs": [], "source": [ - "system.part[:].v = (0, 0, 0)" + "system.part.all().v = (0, 0, 0)" ] }, { diff --git a/samples/chamber_game.py b/samples/chamber_game.py index 42d04855e98..3a192fac16b 100644 --- a/samples/chamber_game.py +++ b/samples/chamber_game.py @@ -266,7 +266,7 @@ new_part.remove() n -= 1 -p_bubbles = np.where(system.part[:].type == bubble_type)[0] +p_bubbles = system.part.select(type=bubble_type) # TEMP CHANGE PARTICLES bpos = [np.random.random() * (pore_xl - snake_head_sigma * 4) + @@ -291,7 +291,7 @@ max_displacement=0.01) system.integrator.run(10000) system.integrator.set_vv() -p_startpos = system.part[:].pos +p_startpos = system.part.all().pos # THERMOSTAT system.thermostat.set_langevin(kT=temperature, gamma=gamma, seed=42) @@ -344,7 +344,7 @@ def set_particle_force(): def restart(): - system.part[:].pos = p_startpos + system.part.all().pos = p_startpos system.galilei.kill_particle_motion() system.galilei.kill_particle_forces() @@ -358,7 +358,7 @@ def explode(): if not exploding: exploding = True expl_time = time.time() - for p in system.part[p_bubbles]: + for p in p_bubbles: dv = p.pos - p_head.pos lv = np.linalg.norm(dv) if lv < expl_range: @@ -478,14 +478,14 @@ def T_to_g(temp): temp_r -= dtemp if temp_r < 0: temp_r = 0.0 - for p in system.part[p_bubbles]: + for p in p_bubbles: if p.pos[0] > pore_xr: p.v = [0, 0, 0] else: temp_l -= dtemp if temp_l < 0: temp_l = 0.0 - for p in system.part[p_bubbles]: + for p in p_bubbles: if p.pos[0] < pore_xl: p.v = [0, 0, 0] diff --git a/samples/dancing.py b/samples/dancing.py index a6dde2e38fd..ca7be18894a 100644 --- a/samples/dancing.py +++ b/samples/dancing.py @@ -59,7 +59,7 @@ gravity = espressomd.constraints.Gravity(g=[0, -1, 0]) system.constraints.add(gravity) -obs = espressomd.observables.ParticlePositions(ids=system.part[:].id) +obs = espressomd.observables.ParticlePositions(ids=system.part.all().id) acc = espressomd.accumulators.TimeSeries(obs=obs, delta_N=1) system.auto_update_accumulators.add(acc) acc.update() diff --git a/samples/electrophoresis.py b/samples/electrophoresis.py index def5c850c30..b0abdd665dc 100644 --- a/samples/electrophoresis.py +++ b/samples/electrophoresis.py @@ -111,8 +111,9 @@ q=np.resize((1, -1), N_IONS), type=np.resize((1, 2), N_IONS)) -logging.info(f"particle types: {system.part[:].type}\n") -logging.info(f"total charge: {np.sum(system.part[:].q)}") +all_partcls = system.part.all() +logging.info(f"particle types: {all_partcls.type}\n") +logging.info(f"total charge: {np.sum(all_partcls.q)}") # warm-up integration ############################################################### @@ -138,8 +139,8 @@ # apply external force (external electric field) ############################################################# n_part = len(system.part) -system.part[:].ext_force = np.dstack( - (system.part[:].q * np.ones(n_part) * E_FIELD, np.zeros(n_part), np.zeros(n_part)))[0] +all_partcls.ext_force = np.dstack( + (all_partcls.q * np.ones(n_part) * E_FIELD, np.zeros(n_part), np.zeros(n_part)))[0] # equilibration ############################################################# @@ -169,7 +170,7 @@ if i % 100 == 0: logging.info(f"\rsampling: {i:4d}") system.integrator.run(N_INT_STEPS) - pos[i] = system.part[:N_MONOMERS].pos + pos[i] = system.part.by_ids(range(N_MONOMERS)).pos logging.info("\nsampling finished!\n") diff --git a/samples/gibbs_ensemble/gibbs_ensemble.py b/samples/gibbs_ensemble/gibbs_ensemble.py index 6eb26962ef1..8f0b991704c 100644 --- a/samples/gibbs_ensemble/gibbs_ensemble.py +++ b/samples/gibbs_ensemble/gibbs_ensemble.py @@ -61,7 +61,8 @@ def __init__(self, @property def density(self): - return len(self.system.part[:].id) / float(np.prod(self.system.box_l)) + return len(self.system.part.all().id) / \ + float(np.prod(self.system.box_l)) @property def energy(self): @@ -76,12 +77,12 @@ def end_simulation(self): def send_info(self): self.pipe.send([MessageID.INFO, self.system.box_l[0], - len(self.system.part[:].id)]) + len(self.system.part.all().id)]) def move_particle(self, DX_MAX): """ Move random particle inside the box """ - random_particle_id = np.random.choice(self.system.part[:].id) - self.old_particle = self.system.part[random_particle_id] + random_particle_id = np.random.choice(self.system.part.all().id) + self.old_particle = self.system.part.by_id(random_particle_id) self.old_pos = self.old_particle.pos self.old_particle.pos = self.old_pos + \ (0.5 - np.random.random(size=3)) * DX_MAX @@ -105,9 +106,10 @@ def revert_add_particle(self): def remove_particle(self): """ Remove a random particle """ - random_particle_id = np.random.choice(self.system.part[:].id) - self.old_particle = self.system.part[random_particle_id].to_dict() - self.system.part[self.old_particle["id"]].remove() + random_particle_id = np.random.choice(self.system.part.all().id) + self.old_particle = self.system.part.by_id( + random_particle_id).to_dict() + self.system.part.by_id(self.old_particle["id"]).remove() self.send_energy() def revert_remove_particle(self): diff --git a/samples/immersed_boundary/addBending.py b/samples/immersed_boundary/addBending.py index 4b1b2a53f16..c39d23dcab3 100644 --- a/samples/immersed_boundary/addBending.py +++ b/samples/immersed_boundary/addBending.py @@ -35,4 +35,4 @@ def AddBending(system, kb): tribend = espressomd.interactions.IBM_Tribend( ind1=id1, ind2=id2, ind3=id3, ind4=id4, kb=kb, refShape="Initial") system.bonded_inter.add(tribend) - system.part[id1].add_bond((tribend, id2, id3, id4)) + system.part.by_id(id1).add_bond((tribend, id2, id3, id4)) diff --git a/samples/immersed_boundary/addSoft.py b/samples/immersed_boundary/addSoft.py index c8397da7537..0beaa96b118 100644 --- a/samples/immersed_boundary/addSoft.py +++ b/samples/immersed_boundary/addSoft.py @@ -48,4 +48,4 @@ def AddSoft(system, comX, comY, comZ, k1, k2): ind1=id1, ind2=id2, ind3=id3, elasticLaw="Skalak", maxDist=5, k1=k1, k2=k2) system.bonded_inter.add(tri) - system.part[id1].add_bond((tri, id2, id3)) + system.part.by_id(id1).add_bond((tri, id2, id3)) diff --git a/samples/immersed_boundary/addVolCons.py b/samples/immersed_boundary/addVolCons.py index 864976f6581..d4817c5a5dc 100644 --- a/samples/immersed_boundary/addVolCons.py +++ b/samples/immersed_boundary/addVolCons.py @@ -26,5 +26,5 @@ def AddVolCons(system, kV): system.bonded_inter.add(volCons) # loop over particles and add - for i in range(len(system.part)): - system.part[i].add_bond((volCons,)) + for p in system.part: + p.add_bond((volCons,)) diff --git a/samples/immersed_boundary/writeVTK.py b/samples/immersed_boundary/writeVTK.py index 719cb6902f3..7f652dd1a22 100644 --- a/samples/immersed_boundary/writeVTK.py +++ b/samples/immersed_boundary/writeVTK.py @@ -31,5 +31,5 @@ def WriteVTK(system, outFile): fp.write(f"POINTS {numPoints} floats\n") # points, positions - for i in range(0, len(system.part)): - fp.write(f"{' '.join(map(str, system.part[i].pos_folded))}\n") + for p in system.part: + fp.write(f"{' '.join(map(str, p.pos_folded))}\n") diff --git a/samples/lj-demo.py b/samples/lj-demo.py index efa23ab4ea4..46173f31cbf 100644 --- a/samples/lj-demo.py +++ b/samples/lj-demo.py @@ -453,9 +453,8 @@ def main_loop(): new_box = np.ones(3) * controls.volume**(1. / 3.) if np.any(np.array(system.box_l) != new_box): - for i in range(len(system.part)): - system.part[i].pos = system.part[i].pos * \ - new_box / system.box_l[0] + for p in system.part: + p.pos *= new_box / system.box_l[0] print("volume changed") system.force_cap = lj_cap system.box_l = new_box @@ -468,7 +467,7 @@ def main_loop(): system.force_cap = lj_cap elif new_part < len(system.part): for i in range(new_part, len(system.part)): - system.part[i].remove() + system.part.by_id(i).remove() print("particles removed") plt1_x_data = plot1.get_xdata() diff --git a/samples/load_checkpoint.py b/samples/load_checkpoint.py index 8e1167efd0c..480aae15d3c 100644 --- a/samples/load_checkpoint.py +++ b/samples/load_checkpoint.py @@ -55,7 +55,7 @@ # test "system.part" print("\n### system.part test ###") -print(f"system.part[:].pos = {system.part[:].pos}") +print(f"system.part.all().pos = {system.part.all().pos}") # test "system.thermostat" print("\n### system.thermostat test ###") diff --git a/samples/rigid_body.py b/samples/rigid_body.py index 036d65f2962..a70f7ecb9af 100644 --- a/samples/rigid_body.py +++ b/samples/rigid_body.py @@ -60,7 +60,7 @@ class ParticleTypes(enum.IntEnum): type=ParticleTypes.BRANCH.value) system.part.add(pos=center - (n + 1) * direction, type=ParticleTypes.BRANCH.value) - +all_partcls = system.part.all() center_of_mass = system.analysis.center_of_mass( p_type=ParticleTypes.BRANCH.value) @@ -68,12 +68,12 @@ class ParticleTypes(enum.IntEnum): # if using multiple nodes, we need to change min_global_cut to the largest # separation -max_inter = np.max(np.linalg.norm(system.part[:].pos - center_of_mass, axis=1)) +max_inter = np.max(np.linalg.norm(all_partcls.pos - center_of_mass, axis=1)) system.min_global_cut = max_inter principal_moments, principal_axes = espressomd.rotation.diagonalized_inertia_tensor( - system.part[:].pos, system.part[:].mass) + all_partcls.pos, all_partcls.mass) # in this simple case, the cluster has principal axes aligned with the box print(f"Principal moments: {principal_moments}") print(f"Principal axes tensor:\n{principal_axes}") @@ -94,7 +94,7 @@ def rotate_vector(vector, axis, angle): principal_moments, principal_axes = espressomd.rotation.diagonalized_inertia_tensor( - system.part[:].pos, system.part[:].mass) + all_partcls.pos, all_partcls.mass) # after rotating the whole object the principal axes changed print("After rotating:") print(f"Principal moments: {principal_moments}") diff --git a/samples/slice_input.py b/samples/slice_input.py index 6b7367ad165..1c62a0bc8b0 100644 --- a/samples/slice_input.py +++ b/samples/slice_input.py @@ -57,32 +57,33 @@ pos_list = np.random.random((n_part, 3)) * system.box_l type_list = np.ones(n_part, dtype=int) -system.part.add(id=id_list, pos=pos_list, type=type_list) +partcls = system.part.add(id=id_list, pos=pos_list, type=type_list) +p0p1 = system.part.by_ids([0, 1]) +print("TYPE\n%s" % partcls.type) +p0p1.type = [3, 3] +print("TYPE_NEW\n%s" % partcls.type) -print("TYPE\n%s" % system.part[:].type) -system.part[0:2].type = [3, 3] -print("TYPE_NEW\n%s" % system.part[:].type) +print("POS\n%s" % partcls.pos) +system.part.by_ids(range(5)).pos = [[1, 1, 1], [2, 2, 2], [ + 3, 3, 3], [4, 4, 4], [5, 5, 5]] +print("POS_NEW\n%s" % partcls.pos) -print("POS\n%s" % system.part[:].pos) -system.part[:5].pos = [[1, 1, 1], [2, 2, 2], [3, 3, 3], [4, 4, 4], [5, 5, 5]] -print("POS_NEW\n%s" % system.part[:].pos) +print("V\n%s" % partcls.v) +p0p1.v = [[1, 2, 3], [2, 3, 4]] +print("V_NEW\n%s" % partcls.v) -print("V\n%s" % system.part[:].v) -system.part[:2].v = [[1, 2, 3], [2, 3, 4]] -print("V_NEW\n%s" % system.part[:].v) - -print("F\n%s" % system.part[:].f) -system.part[:2].f = [[3, 4, 5], [4, 5, 6]] -print("F_NEW\n%s" % system.part[:].f) +print("F\n%s" % partcls.f) +p0p1.f = [[3, 4, 5], [4, 5, 6]] +print("F_NEW\n%s" % partcls.f) if espressomd.has_features(["MASS"]): - print("MASS\n%s" % system.part[:].mass) - system.part[:2].mass = [2, 3] - print("MASS_NEW\n%s" % system.part[:].mass) + print("MASS\n%s" % partcls.mass) + p0p1.mass = [2, 3] + print("MASS_NEW\n%s" % partcls.mass) if espressomd.has_features(["ELECTROSTATICS"]): - print("Q\n%s" % system.part[:].q) - system.part[::2].q = np.ones(n_part // 2) - system.part[1::2].q = -np.ones(n_part // 2) - print("Q_NEW\n%s" % system.part[:].q) + print("Q\n%s" % partcls.q) + system.part.by_ids(range(0, n_part, 2)).q = np.ones(n_part // 2) + system.part.by_ids(range(1, n_part, 2)).q = -np.ones(n_part // 2) + print("Q_NEW\n%s" % partcls.q) diff --git a/samples/visualization_elc.py b/samples/visualization_elc.py index b38107b517c..3038448c723 100644 --- a/samples/visualization_elc.py +++ b/samples/visualization_elc.py @@ -46,13 +46,12 @@ system.time_step = 0.02 system.cell_system.skin = 0.4 -system.periodicity = [1, 1, 1] +system.periodicity = 3 * [True] qion = 1 for i in range(300): rpos = np.random.random(3) * box_l - system.part.add(pos=rpos, type=0, q=qion) - qion *= -1 + system.part.add(pos=rpos, type=0, q=(-1)**i * qion) system.constraints.add(shape=espressomd.shapes.Wall( dist=0, normal=[0, 0, 1]), particle_type=1) diff --git a/samples/visualization_interactive.py b/samples/visualization_interactive.py index 54e6076a99e..d8bc2b0e658 100644 --- a/samples/visualization_interactive.py +++ b/samples/visualization_interactive.py @@ -38,8 +38,7 @@ system.cell_system.skin = 3.0 N = 50 -for i in range(N): - system.part.add(pos=[0, 0, 0]) +partcls = system.part.add(pos=N * [[0, 0, 0]]) system.thermostat.set_langevin(kT=1.0, gamma=1.0, seed=42) @@ -47,9 +46,9 @@ def spin(): - system.part[:].pos = [[box_l * 0.5, box_l * (i + 1) / (N + 2), box_l * 0.5] - for i in range(N)] - system.part[:].v = [ + partcls.pos = [[box_l * 0.5, box_l * (i + 1) / (N + 2), box_l * 0.5] + for i in range(N)] + partcls.v = [ [np.sin(10.0 * i / N) * 20, 0, np.cos(10.0 * i / N) * 20] for i in range(N)] diff --git a/src/python/espressomd/drude_helpers.py b/src/python/espressomd/drude_helpers.py index 9e508ab98bb..f481bb4babb 100644 --- a/src/python/espressomd/drude_helpers.py +++ b/src/python/espressomd/drude_helpers.py @@ -232,8 +232,8 @@ def setup_and_add_drude_exclusion_bonds(self, system, verbose=False): for drude_id in self.drude_id_list: core_id = self.core_id_from_drude_id[drude_id] - pd = system.part[drude_id] - pc = system.part[core_id] + pd = system.part.by_id(drude_id) + pc = system.part.by_id(core_id) bond = self.drude_dict[pd.type]["subtr_sr_bonds_drude-core"] pc.add_bond((bond, drude_id)) if verbose: diff --git a/src/python/espressomd/io/writer/vtf.py b/src/python/espressomd/io/writer/vtf.py index de39249ec1a..529d9f7ed02 100644 --- a/src/python/espressomd/io/writer/vtf.py +++ b/src/python/espressomd/io/writer/vtf.py @@ -66,12 +66,14 @@ def writevsf(system, fp, types='all'): fp.write(f"unitcell {' '.join(map(str, system.box_l))}\n") for pid, vtf_id, in vtf_index.items(): + partcl = system.part.by_id(pid) fp.write( - f"atom {vtf_id} radius 1 name {system.part[pid].type} type {system.part[pid].type} \n") + f"atom {vtf_id} radius 1 name {partcl.type} type {partcl.type} \n") for pid, vtf_id, in vtf_index.items(): - for b in system.part[pid].bonds: - if system.part[b[1]].id in vtf_index: - fp.write(f"bond {vtf_id}:{vtf_index[system.part[b[1]].id]}\n") + for b in system.part.by_id(pid).bonds: + if system.part.by_id(b[1]).id in vtf_index: + fp.write( + f"bond {vtf_id}:{vtf_index[system.part.by_id(b[1]).id]}\n") def writevcf(system, fp, types='all'): @@ -91,4 +93,4 @@ def writevcf(system, fp, types='all'): vtf_index = vtf_pid_map(system, types) fp.write("\ntimestep indexed\n") for pid, vtf_id, in vtf_index.items(): - fp.write(f"{vtf_id} {' '.join(map(str, system.part[pid].pos))}\n") + fp.write(f"{vtf_id} {' '.join(map(str, system.part.by_id(pid).pos))}\n") diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx index c3d2d0f6168..f467eaf1b2d 100644 --- a/src/python/espressomd/particle_data.pyx +++ b/src/python/espressomd/particle_data.pyx @@ -1565,7 +1565,7 @@ cdef class _ParticleSliceImpl: pl = ParticleList() for i in self.id_selection: if pl.exists(i): - res += f"{pl[i]}, " + res += f"{pl.by_id(i)}, " # Remove final comma return f"ParticleSlice([{res[:-2]}])" @@ -1658,20 +1658,28 @@ class ParticleSlice(_ParticleSliceImpl): cdef class ParticleList: """ - Provides access to the particles via ``[i]``, where ``i`` is the particle - id. Returns a :class:`espressomd.particle_data.ParticleHandle` object. + Provides access to the particles. """ - # Retrieve a particle + def by_id(self, id): + """ + Access a particle by its integer id. + """ + return ParticleHandle(id) - def __getitem__(self, key): - # Single particle id results in a ParticleHandle, everything else - # in a ParticleSlice - if is_valid_type(key, int): - return ParticleHandle(key) - else: - return ParticleSlice(key) + def by_ids(self, ids): + """ + Get a slice of particles by their integer ids. + """ + return ParticleSlice(ids) + + def all(self): + """ + Get a slice containing all particles. + """ + all_ids = get_particle_ids() + return self.by_ids(all_ids) # __getstate__ and __setstate__ define the pickle interaction def __getstate__(self): @@ -1709,7 +1717,7 @@ cdef class ParticleList: self._place_new_particle(params[particle_number]) IF EXCLUSIONS: for pid in exclusions: - self[pid].exclusions = exclusions[pid] + self.by_id(pid).exclusions = exclusions[pid] def __len__(self): return get_n_part() @@ -1787,7 +1795,7 @@ dip is sufficient as the length of the vector defines the scalar dipole moment." Setting dip overwrites the rotation of the particle around the dipole axis. \ Set quat and scalar dipole moment (dipm) instead.") - # The ParticleList[]-getter ist not valid yet, as the particle + # The ParticleList can not be used yet, as the particle # doesn't yet exist. Hence, the setting of position has to be # done here. the code is from the pos:property of ParticleHandle check_type_or_throw_except( @@ -1801,9 +1809,9 @@ Set quat and scalar dipole moment (dipm) instead.") del P["id"] if P != {}: - self[pid].update(P) + self.by_id(pid).update(P) - return self[pid] + return self.by_id(pid) def _place_new_particles(self, Ps): # Check if all entries have the same length @@ -1826,14 +1834,14 @@ Set quat and scalar dipole moment (dipm) instead.") self._place_new_particle(P) # Return slice of added particles - return self[Ps["id"]] + return self.by_ids(Ps["id"]) # Iteration over all existing particles def __iter__(self): ids = get_particle_ids() for i in ids: - yield self[i] + yield self.by_id(i) def exists(self, idx): if is_valid_type(idx, int): @@ -1954,7 +1962,7 @@ Set quat and scalar dipole moment (dipm) instead.") for i in ids: for j in ids[i + 1:]: - yield (self[i], self[j]) + yield (self.by_id(i), self.by_id(j)) def select(self, *args, **kwargs): """Generates a particle slice by filtering particles via a user-defined criterion diff --git a/src/python/espressomd/polymer.pyx b/src/python/espressomd/polymer.pyx index 719ca533fbe..40e5ae68536 100644 --- a/src/python/espressomd/polymer.pyx +++ b/src/python/espressomd/polymer.pyx @@ -289,12 +289,12 @@ def setup_diamond_polymer(system=None, bond=None, MPC=0, # add bonds along the chain if not no_bonds: if j == 1: - system.part[node_ids[start_node]].add_bond( + system.part.by_id(node_ids[start_node]).add_bond( (bond, current_id)) else: - system.part[current_id].add_bond( + system.part.by_id(current_id).add_bond( (bond, current_id - 1)) if j == MPC: - system.part[node_ids[end_node]].add_bond( + system.part.by_id(node_ids[end_node]).add_bond( (bond, current_id)) current_id += 1 diff --git a/src/python/object_in_fluid/oif_classes.py b/src/python/object_in_fluid/oif_classes.py index c6c1d8e54dc..314595a4838 100644 --- a/src/python/object_in_fluid/oif_classes.py +++ b/src/python/object_in_fluid/oif_classes.py @@ -489,7 +489,7 @@ def copy(self, origin=None, particle_type=-1, particle_mass=1.0, new_part_id = len(self.system.part) self.system.part.add( pos=tmp_pos, type=particle_type, mass=particle_mass, mol_id=particle_type) - new_part = self.system.part[new_part_id] + new_part = self.system.part.by_id(new_part_id) new_part_point = PartPoint(new_part, len(mesh.points), new_part_id) mesh.points.append(new_part_point) for edge in self.edges: diff --git a/testsuite/python/accumulator_mean_variance.py b/testsuite/python/accumulator_mean_variance.py index 78d7fd1ae4c..0ebfdc4bedd 100644 --- a/testsuite/python/accumulator_mean_variance.py +++ b/testsuite/python/accumulator_mean_variance.py @@ -58,7 +58,7 @@ def test_accumulator(self): positions = np.copy(system.box_l) * np.random.random((10, N_PART, 3)) for pos in positions: - system.part[:].pos = pos + system.part.all().pos = pos acc.update() self.assertEqual(acc.get_params()['obs'], obs) diff --git a/testsuite/python/accumulator_time_series.py b/testsuite/python/accumulator_time_series.py index 265f5eb0ebe..cd9c0377927 100644 --- a/testsuite/python/accumulator_time_series.py +++ b/testsuite/python/accumulator_time_series.py @@ -58,7 +58,7 @@ def test_accumulator(self): positions = np.copy(system.box_l) * np.random.random((10, N_PART, 3)) for pos in positions: - system.part[:].pos = pos + system.part.all().pos = pos acc.update() self.assertEqual(acc.get_params()['obs'], obs) diff --git a/testsuite/python/analyze_chains.py b/testsuite/python/analyze_chains.py index 652737668f7..655a6e7de48 100644 --- a/testsuite/python/analyze_chains.py +++ b/testsuite/python/analyze_chains.py @@ -24,7 +24,7 @@ @utx.skipIfMissingFeatures("LENNARD_JONES") class AnalyzeChain(ut.TestCase): - system = espressomd.System(box_l=[1.0, 1.0, 1.0]) + system = espressomd.System(box_l=3 * [1.0]) np.random.seed(1234) num_poly = 2 @@ -34,7 +34,7 @@ class AnalyzeChain(ut.TestCase): def setUpClass(cls): box_l = 20.0 # start with a small box - cls.system.box_l = np.array([box_l, box_l, box_l]) + cls.system.box_l = np.array(3 * [box_l]) cls.system.cell_system.set_n_square(use_verlet_lists=False) fene = espressomd.interactions.FeneBond(k=30, d_r_max=2) cls.system.bonded_inter.add(fene) @@ -43,28 +43,27 @@ def setUpClass(cls): beads_per_chain=cls.num_mono, seed=42) for p in positions: for ndx, m in enumerate(p): - part_id = len(cls.system.part) - cls.system.part.add(id=part_id, pos=m) + partcl = cls.system.part.add(pos=m) if ndx > 0: - cls.system.part[part_id].add_bond((fene, part_id - 1)) + partcl.add_bond((fene, partcl.id - 1)) # bring two polymers to opposite corners: # far in cell centre, but mirror images are close - head_id = 0 + pol_0_parts = cls.system.part.by_ids(range(0, cls.num_mono)) + cm = np.mean(pol_0_parts.pos, axis=0) + pol_0_parts.pos += (- cm + cls.system.box_l) + head_id = cls.num_mono tail_id = head_id + cls.num_mono - cm = np.mean(cls.system.part[head_id:tail_id].pos, axis=0) - cls.system.part[head_id:tail_id].pos = cls.system.part[ - head_id:tail_id].pos - cm + cls.system.box_l - head_id = cls.num_mono + 1 - tail_id = head_id + cls.num_mono - cm = np.mean(cls.system.part[head_id:tail_id].pos, axis=0) - cls.system.part[head_id:tail_id].pos -= cm + pol_1_parts = cls.system.part.by_ids(range(head_id, tail_id)) + cm = np.mean(pol_1_parts.pos, axis=0) + pol_1_parts.pos -= cm # python version of the espresso core function, # does not check mirror distances def calc_re(self): - head_id = np.arange(0, self.num_poly * self.num_mono, self.num_mono) - tail_id = head_id + self.num_mono - 1 - dist = self.system.part[head_id].pos - self.system.part[tail_id].pos + head_ids = np.arange(0, self.num_poly * self.num_mono, self.num_mono) + tail_ids = head_ids + self.num_mono - 1 + dist = self.system.part.by_ids( + head_ids).pos - self.system.part.by_ids(tail_ids).pos dist = np.sum(dist**2, axis=-1) return np.mean(np.sqrt(dist)), np.std( np.sqrt(dist)), np.mean(dist), np.std(dist) @@ -72,12 +71,12 @@ def calc_re(self): # python version of the espresso core function, # does not check mirror distances def calc_rg(self): - head_id = np.arange(0, self.num_poly * self.num_mono, self.num_mono) - tail_id = head_id + self.num_mono - 1 + head_ids = np.arange(0, self.num_poly * self.num_mono, self.num_mono) + tail_ids = head_ids + self.num_mono rg2 = [] for p in range(self.num_poly): rg2.append( - np.var(self.system.part[head_id[p]:tail_id[p] + 1].pos, axis=0)) + np.var(self.system.part.by_ids(range(head_ids[p], tail_ids[p])).pos, axis=0)) rg2 = np.array(rg2) rg2 = np.sum(rg2, axis=1) return np.mean(np.sqrt(rg2)), np.std( @@ -86,11 +85,16 @@ def calc_rg(self): # python version of the espresso core function, # does not check mirror distances def calc_rh(self): - head_id = np.arange(0, self.num_poly * self.num_mono, self.num_mono) - tail_id = head_id + self.num_mono - 1 + head_ids = np.arange(0, self.num_poly * self.num_mono, self.num_mono) + tail_ids = head_ids + self.num_mono - 1 rh = [] for p in range(self.num_poly): - r = np.array(self.system.part[head_id[p]:tail_id[p] + 1].pos) + r = np.array( + self.system.part.by_ids( + range( + head_ids[p], + tail_ids[p] + + 1)).pos) # this generates indices for all i 0.5 * self.system.box_l, self.system.box_l - dist, dist) @@ -56,27 +56,27 @@ def nbhood(self, pos, r_catch): # python version of the espresso core function, using pos def dist_to_pos(self, pos): - dist = np.fabs(self.system.part[:].pos - pos) + dist = np.fabs(self.system.part.all().pos - pos) # check smaller distances via PBC dist = np.where( dist > 0.5 * self.system.box_l, self.system.box_l - dist, dist) - dist = np.sum(dist**2, axis=-1) - return np.sqrt(np.min(dist)) + dist = np.linalg.norm(dist, axis=-1) + return np.min(dist) # python version of the espresso core function, using id def dist_to_id(self, id): dist = np.fabs( - np.delete(self.system.part[:].pos, id, axis=0) - self.system.part[id].pos) + np.delete(self.system.part.all().pos, id, axis=0) - self.system.part.by_id(id).pos) # check smaller distances via PBC dist = np.where( dist > 0.5 * self.system.box_l, self.system.box_l - dist, dist) - dist = np.sum(dist**2, axis=1) - return np.sqrt(np.min(dist)) + dist = np.linalg.norm(dist, axis=1) + return np.min(dist) def test_min_dist(self): # try five times for _ in range(5): - self.system.part[:].pos = np.random.random( + self.system.part.all().pos = np.random.random( (len(self.system.part), 3)) * BOX_L self.assertAlmostEqual(self.system.analysis.min_dist(), self.min_dist(), @@ -89,7 +89,7 @@ def test_min_dist_empty(self): def test_nbhood(self): # try five times for i in range(1, 10, 2): - self.system.part[:].pos = np.random.random( + self.system.part.all().pos = np.random.random( (len(self.system.part), 3)) * BOX_L np.testing.assert_allclose( self.system.analysis.nbhood([i, i, i], i * 2), @@ -99,7 +99,7 @@ def test_distance_to_pos(self): parts = self.system.part # try five times for i in range(5): - self.system.part[:].pos = np.random.random( + self.system.part.all().pos = np.random.random( (len(self.system.part), 3)) * BOX_L self.assertAlmostEqual( np.min([self.system.distance([i, i, i], p.pos) @@ -114,16 +114,17 @@ def test_distance_to_part(self): parts = self.system.part # try five times for i in range(5): - self.system.part[:].pos = np.random.random( + part_i = parts.by_id(i) + self.system.part.all().pos = np.random.random( (len(self.system.part), 3)) * BOX_L self.assertAlmostEqual( - np.min([self.system.distance(parts[i], p.pos) + np.min([self.system.distance(part_i, p.pos) for p in parts if p.id != i]), self.dist_to_id(i), delta=1e-10) self.assertEqual( - np.min([self.system.distance(parts[i], p.pos) + np.min([self.system.distance(part_i, p.pos) for p in parts if p.id != i]), - np.min([self.system.distance(p.pos, parts[i]) for p in parts if p.id != i])) + np.min([self.system.distance(p.pos, part_i) for p in parts if p.id != i])) if __name__ == "__main__": diff --git a/testsuite/python/analyze_distribution.py b/testsuite/python/analyze_distribution.py index 85a1bd5bece..ec91669a460 100644 --- a/testsuite/python/analyze_distribution.py +++ b/testsuite/python/analyze_distribution.py @@ -32,7 +32,7 @@ def setUpClass(cls): # start with a small box cls.system.box_l = np.array([box_l, box_l, box_l]) cls.system.cell_system.set_n_square(use_verlet_lists=False) - cls.system.part.add( + cls.partcls = cls.system.part.add( pos=np.outer(np.random.random(cls.num_part), cls.system.box_l)) def calc_min_distribution(self, bins): @@ -46,9 +46,9 @@ def calc_min_distribution(self, bins): # test system.analysis.distribution(), all the same particle types def test_distribution_lin(self): # increase PBC to remove mirror images - old_pos = self.system.part[:].pos.copy() + old_pos = self.partcls.pos.copy() self.system.box_l = self.system.box_l * 2. - self.system.part[:].pos = old_pos + self.partcls.pos = old_pos r_min = 0.0 r_max = 100.0 r_bins = 100 diff --git a/testsuite/python/analyze_energy.py b/testsuite/python/analyze_energy.py index 29ccadffb37..6e1e50cc589 100644 --- a/testsuite/python/analyze_energy.py +++ b/testsuite/python/analyze_energy.py @@ -53,7 +53,7 @@ def tearDown(self): self.system.part.clear() def test_kinetic(self): - p0, p1 = self.system.part[:] + p0, p1 = self.system.part.all() p0.pos = [1, 2, 2] p1.pos = [5, 2, 2] # single moving particle @@ -73,7 +73,7 @@ def test_kinetic(self): self.assertAlmostEqual(energy["non_bonded"], 0., delta=1e-7) def test_non_bonded(self): - p0, p1 = self.system.part[:] + p0, p1 = self.system.part.all() p0.pos = [1, 2, 2] p1.pos = [2, 2, 2] energy = self.system.analysis.energy() @@ -95,7 +95,7 @@ def test_non_bonded(self): + energy["non_bonded", 1, 1], energy["total"], delta=1e-7) def test_bonded(self): - p0, p1 = self.system.part[:] + p0, p1 = self.system.part.all() p0.pos = [1, 2, 2] p1.pos = [3, 2, 2] # single bond @@ -122,7 +122,7 @@ def test_bonded(self): self.assertAlmostEqual(energy["non_bonded"], 0., delta=1e-7) def test_all(self): - p0, p1 = self.system.part[:] + p0, p1 = self.system.part.all() p0.pos = [1, 2, 2] p1.pos = [2, 2, 2] p0.v = [3, 4, 5] @@ -153,7 +153,7 @@ def test_all(self): @utx.skipIfMissingFeatures(["ELECTROSTATICS", "P3M"]) def test_electrostatics(self): - p0, p1 = self.system.part[:] + p0, p1 = self.system.part.all() p0.pos = [1, 2, 2] p1.pos = [3, 2, 2] p0.q = 1 diff --git a/testsuite/python/analyze_gyration_tensor.py b/testsuite/python/analyze_gyration_tensor.py index 304bf4151a5..523c49db236 100644 --- a/testsuite/python/analyze_gyration_tensor.py +++ b/testsuite/python/analyze_gyration_tensor.py @@ -60,7 +60,7 @@ def test_gyration_tensor(self): res = self.system.analysis.gyration_tensor( p_type=[self.type_stick, self.type_cube]) rg = self.system.analysis.calc_rg( - chain_start=0, number_of_chains=1, chain_length=len(self.system.part[:]))[0] + chain_start=0, number_of_chains=1, chain_length=len(self.system.part.all()))[0] # check eigenvectors np.testing.assert_allclose( np.abs(res['eva0'][1]), [0., 0., 1.], atol=1e-6) diff --git a/testsuite/python/analyze_mass_related.py b/testsuite/python/analyze_mass_related.py index 3d09569ed16..1b97b19698e 100644 --- a/testsuite/python/analyze_mass_related.py +++ b/testsuite/python/analyze_mass_related.py @@ -44,14 +44,15 @@ def setUpClass(cls): if espressomd.has_features("VIRTUAL_SITES"): cls.system.part.add( pos=np.random.random((10, 3)) * cls.box_l, type=[0] * 10, virtual=[True] * 10) - cls.system.part[:].v = np.random.random( - (len(cls.system.part), 3)) + 0.1 + all_partcls = cls.system.part.all() + all_partcls.v = np.random.random( + (len(all_partcls), 3)) + 0.1 if espressomd.has_features("MASS"): - cls.system.part[:].mass = 0.5 + \ - np.random.random(len(cls.system.part)) + all_partcls.mass = 0.5 + \ + np.random.random(len(all_partcls)) def i_tensor(self, ids): - pslice = self.system.part[ids] + pslice = self.system.part.by_ids(ids) I = np.zeros((3, 3)) diff --git a/testsuite/python/array_properties.py b/testsuite/python/array_properties.py index 6e450f1cffa..bc337f5f99a 100644 --- a/testsuite/python/array_properties.py +++ b/testsuite/python/array_properties.py @@ -99,7 +99,7 @@ class ArrayPropertyTest(ArrayCommon): system.box_l = [12.0, 12.0, 12.0] system.time_step = 0.01 system.cell_system.skin = 0.01 - system.part.add(pos=[0, 0, 0]) + partcl = system.part.add(pos=[0, 0, 0]) def setUp(self): self.system.box_l = [12.0, 12.0, 12.0] @@ -112,77 +112,77 @@ def assert_copy_is_writable(self, array): self.assertTrue(cpy.flags.writeable) def test_common(self): - self.assert_operator_usage_raises(self.system.part[0].pos) - self.assert_operator_usage_raises(self.system.part[0].v) - self.assert_operator_usage_raises(self.system.part[0].f) - self.assert_operator_usage_raises(self.system.part[0].pos_folded) + self.assert_operator_usage_raises(self.partcl.pos) + self.assert_operator_usage_raises(self.partcl.v) + self.assert_operator_usage_raises(self.partcl.f) + self.assert_operator_usage_raises(self.partcl.pos_folded) self.assert_operator_usage_raises(self.system.box_l) - check_array_writable(self.system.part[0].pos) - check_array_writable(self.system.part[0].v) - check_array_writable(self.system.part[0].f) + check_array_writable(self.partcl.pos) + check_array_writable(self.partcl.v) + check_array_writable(self.partcl.f) check_array_writable(self.system.box_l) - self.assert_copy_is_writable(self.system.part[0].pos) - self.assert_copy_is_writable(self.system.part[0].v) - self.assert_copy_is_writable(self.system.part[0].f) - self.assert_copy_is_writable(self.system.part[0].pos_folded) + self.assert_copy_is_writable(self.partcl.pos) + self.assert_copy_is_writable(self.partcl.v) + self.assert_copy_is_writable(self.partcl.f) + self.assert_copy_is_writable(self.partcl.pos_folded) self.assert_copy_is_writable(self.system.box_l) @utx.skipIfMissingFeatures(["ROTATION"]) def test_rotation(self): - self.assert_operator_usage_raises(self.system.part[0].omega_lab) - self.assert_operator_usage_raises(self.system.part[0].quat) - self.assert_operator_usage_raises(self.system.part[0].rotation) - self.assert_operator_usage_raises(self.system.part[0].omega_body) - self.assert_operator_usage_raises(self.system.part[0].torque_lab) + self.assert_operator_usage_raises(self.partcl.omega_lab) + self.assert_operator_usage_raises(self.partcl.quat) + self.assert_operator_usage_raises(self.partcl.rotation) + self.assert_operator_usage_raises(self.partcl.omega_body) + self.assert_operator_usage_raises(self.partcl.torque_lab) if espressomd.has_features("EXTERNAL_FORCES"): - self.assert_operator_usage_raises(self.system.part[0].ext_torque) + self.assert_operator_usage_raises(self.partcl.ext_torque) - check_array_writable(self.system.part[0].quat) - check_array_writable(self.system.part[0].omega_lab) - check_array_writable(self.system.part[0].rotation) - check_array_writable(self.system.part[0].omega_body) - check_array_writable(self.system.part[0].torque_lab) + check_array_writable(self.partcl.quat) + check_array_writable(self.partcl.omega_lab) + check_array_writable(self.partcl.rotation) + check_array_writable(self.partcl.omega_body) + check_array_writable(self.partcl.torque_lab) if espressomd.has_features("EXTERNAL_FORCES"): - check_array_writable(self.system.part[0].ext_torque) + check_array_writable(self.partcl.ext_torque) - self.assert_copy_is_writable(self.system.part[0].omega_lab) - self.assert_copy_is_writable(self.system.part[0].quat) - self.assert_copy_is_writable(self.system.part[0].rotation) - self.assert_copy_is_writable(self.system.part[0].omega_body) - self.assert_copy_is_writable(self.system.part[0].torque_lab) + self.assert_copy_is_writable(self.partcl.omega_lab) + self.assert_copy_is_writable(self.partcl.quat) + self.assert_copy_is_writable(self.partcl.rotation) + self.assert_copy_is_writable(self.partcl.omega_body) + self.assert_copy_is_writable(self.partcl.torque_lab) if espressomd.has_features("EXTERNAL_FORCES"): - self.assert_copy_is_writable(self.system.part[0].ext_torque) + self.assert_copy_is_writable(self.partcl.ext_torque) @utx.skipIfMissingFeatures(["ROTATIONAL_INERTIA"]) def test_rotational_inertia(self): - self.assert_operator_usage_raises(self.system.part[0].rinertia) - check_array_writable(self.system.part[0].rinertia) - self.assert_copy_is_writable(self.system.part[0].rinertia) + self.assert_operator_usage_raises(self.partcl.rinertia) + check_array_writable(self.partcl.rinertia) + self.assert_copy_is_writable(self.partcl.rinertia) @utx.skipIfMissingFeatures(["EXTERNAL_FORCES"]) def test_external_forces(self): - self.assert_operator_usage_raises(self.system.part[0].ext_force) - self.assert_operator_usage_raises(self.system.part[0].fix) + self.assert_operator_usage_raises(self.partcl.ext_force) + self.assert_operator_usage_raises(self.partcl.fix) - check_array_writable(self.system.part[0].ext_force) - check_array_writable(self.system.part[0].fix) + check_array_writable(self.partcl.ext_force) + check_array_writable(self.partcl.fix) - self.assert_copy_is_writable(self.system.part[0].ext_force) - self.assert_copy_is_writable(self.system.part[0].fix) + self.assert_copy_is_writable(self.partcl.ext_force) + self.assert_copy_is_writable(self.partcl.fix) @utx.skipIfMissingFeatures(["ROTATION", "THERMOSTAT_PER_PARTICLE", "PARTICLE_ANISOTROPY"]) def test_rot_aniso(self): - self.assert_operator_usage_raises(self.system.part[0].gamma_rot) + self.assert_operator_usage_raises(self.partcl.gamma_rot) - check_array_writable(self.system.part[0].gamma_rot) + check_array_writable(self.partcl.gamma_rot) - self.assert_copy_is_writable(self.system.part[0].gamma_rot) + self.assert_copy_is_writable(self.partcl.gamma_rot) def test_lb(self): lbf = espressomd.lb.LBFluid(agrid=0.5, dens=1, visc=1, tau=0.01) @@ -196,23 +196,23 @@ def test_lb(self): @utx.skipIfMissingFeatures(["THERMOSTAT_PER_PARTICLE", "PARTICLE_ANISOTROPY"]) def test_thermostat_per_particle_aniso(self): - self.assert_operator_usage_raises(self.system.part[0].gamma) + self.assert_operator_usage_raises(self.partcl.gamma) - check_array_writable(self.system.part[0].gamma) + check_array_writable(self.partcl.gamma) - self.assert_copy_is_writable(self.system.part[0].gamma) + self.assert_copy_is_writable(self.partcl.gamma) @utx.skipIfMissingFeatures(["DIPOLES"]) def test_dipoles(self): - self.assert_operator_usage_raises(self.system.part[0].dip) + self.assert_operator_usage_raises(self.partcl.dip) - check_array_writable(self.system.part[0].dip) + check_array_writable(self.partcl.dip) - self.assert_copy_is_writable(self.system.part[0].dip) + self.assert_copy_is_writable(self.partcl.dip) @utx.skipIfMissingFeatures(["EXCLUSIONS"]) def test_exclusions(self): - self.assert_operator_usage_raises(self.system.part[0].exclusions) + self.assert_operator_usage_raises(self.partcl.exclusions) def test_partial_periodic(self): self.assert_operator_usage_raises(self.system.periodicity) diff --git a/testsuite/python/auto_exclusions.py b/testsuite/python/auto_exclusions.py index 183775ac42e..d236d6d2f9b 100644 --- a/testsuite/python/auto_exclusions.py +++ b/testsuite/python/auto_exclusions.py @@ -38,21 +38,21 @@ def test_linear(self): s.part.add(id=i, pos=[0, 0, 0]) for i in range(9): - s.part[i].add_bond((bond, i + 1)) + s.part.by_id(i).add_bond((bond, i + 1)) s.auto_exclusions(1) for p in range(1, 9): - excl = s.part[p].exclusions + excl = s.part.by_id(p).exclusions self.assertEqual(len(excl), 2) self.assertIn(p - 1, excl) self.assertIn(p + 1, excl) - excl = s.part[0].exclusions + excl = s.part.by_id(0).exclusions self.assertEqual(len(excl), 1) self.assertIn(1, excl) - excl = s.part[9].exclusions + excl = s.part.by_id(9).exclusions self.assertEqual(len(excl), 1) self.assertIn(8, excl) @@ -65,12 +65,12 @@ def test_ring(self): s.part.add(id=i, pos=[0, 0, 0]) for i in range(10): - s.part[i].add_bond((bond, (i + 1) % 10)) + s.part.by_id(i).add_bond((bond, (i + 1) % 10)) s.auto_exclusions(2) for p in range(10): - excl = s.part[p].exclusions + excl = s.part.by_id(p).exclusions self.assertEqual(len(excl), 4) for i in range(1, 3): self.assertIn((p - i) % 10, excl) diff --git a/testsuite/python/cluster_analysis.py b/testsuite/python/cluster_analysis.py index 416cef9b876..89945dfa1e9 100644 --- a/testsuite/python/cluster_analysis.py +++ b/testsuite/python/cluster_analysis.py @@ -134,12 +134,12 @@ def test_single_cluster_analysis(self): # Longest distance ref_dist = self.system.distance( - self.system.part[0], self.system.part[len(self.system.part) - 1]) + self.system.part.by_id(0), self.system.part.by_id(len(self.system.part) - 1)) self.assertAlmostEqual(c.longest_distance(), ref_dist, delta=1E-8) # Radius of gyration rg = 0. - com_particle = self.system.part[len(self.system.part) // 2] + com_particle = self.system.part.by_id(len(self.system.part) // 2) for p in c.particles(): rg += self.system.distance(p, com_particle)**2 rg /= len(self.system.part) @@ -187,7 +187,7 @@ def test_analysis_for_bonded_particles(self): # Check particle to cluster id mapping, once by ParticleHandle, once by # id self.assertEqual(self.cs.cid_for_particle( - self.system.part[0]), self.cs.cluster_ids()[0]) + self.system.part.by_id(0)), self.cs.cluster_ids()[0]) self.assertEqual(self.cs.cid_for_particle(1), self.cs.cluster_ids()[0]) # Unknown method should return None diff --git a/testsuite/python/collision_detection.py b/testsuite/python/collision_detection.py index 57403fff2d4..b4ab3a5bbdb 100644 --- a/testsuite/python/collision_detection.py +++ b/testsuite/python/collision_detection.py @@ -80,13 +80,13 @@ def test_bind_centers(self): system.collision_detection.set_params(mode="off") system.part.clear() - system.part.add(pos=(0, 0, 0), id=0) - system.part.add(pos=(0.1, 0, 0), id=1) - system.part.add(pos=(0.1, 0.3, 0), id=2) + p0 = system.part.add(pos=(0, 0, 0), id=0) + p1 = system.part.add(pos=(0.1, 0, 0), id=1) + p2 = system.part.add(pos=(0.1, 0.3, 0), id=2) system.integrator.run(1) - self.assertEqual(system.part[0].bonds, ()) - self.assertEqual(system.part[1].bonds, ()) - self.assertEqual(system.part[2].bonds, ()) + self.assertEqual(p0.bonds, ()) + self.assertEqual(p1.bonds, ()) + self.assertEqual(p2.bonds, ()) # Check that it cannot be activated system.collision_detection.set_params( @@ -96,14 +96,14 @@ def test_bind_centers(self): bond0 = ((system.bonded_inter[0], 1),) bond1 = ((system.bonded_inter[0], 0),) self.assertTrue( - system.part[0].bonds == bond0 or system.part[1].bonds == bond1) - self.assertEqual(system.part[2].bonds, ()) + p0.bonds == bond0 or p1.bonds == bond1) + self.assertEqual(p2.bonds, ()) # Check that no additional bonds appear system.integrator.run(1) self.assertTrue( - system.part[0].bonds == bond0 or system.part[1].bonds == bond1) - self.assertEqual(system.part[2].bonds, ()) + p0.bonds == bond0 or p1.bonds == bond1) + self.assertEqual(p2.bonds, ()) # Check turning it off system.collision_detection.set_params(mode="off") @@ -140,7 +140,7 @@ def run_test_bind_at_point_of_collision_for_pos(self, *positions): # Check that nothing explodes when the particles are moved. # In particular for parallel simulations system.thermostat.set_langevin(kT=0, gamma=0.01, seed=42) - system.part[:].v = [0.05, 0.01, 0.15] + system.part.all().v = [0.05, 0.01, 0.15] system.integrator.run(3000) self.verify_state_after_bind_at_poc(expected_np) @@ -165,19 +165,19 @@ def verify_state_after_bind_at_poc(self, expected_np): # Bond type self.assertEqual(p.bonds[0][0], self.H2) # get partner - p2 = system.part[p.bonds[0][1]] + p2 = system.part.by_id(p.bonds[0][1]) # Is that really a vs self.assertTrue(p2.virtual) # Get base particles - base_p1 = system.part[p.vs_relative[0]] - base_p2 = system.part[p2.vs_relative[0]] + base_p1 = system.part.by_id(p.vs_relative[0]) + base_p2 = system.part.by_id(p2.vs_relative[0]) # Take note of accounted-for particles for _p in p, p2, base_p1, base_p2: parts_not_accounted_for.remove(_p.id) self.verify_bind_at_poc_pair(base_p1, base_p2, p, p2) # Check particle that did not take part in collision. self.assertEqual(len(parts_not_accounted_for), 1) - p = system.part[parts_not_accounted_for[0]] + p = system.part.by_id(parts_not_accounted_for[0]) self.assertFalse(p.virtual) self.assertEqual(p.bonds, ()) parts_not_accounted_for.remove(p.id) @@ -216,7 +216,7 @@ def verify_bind_at_poc_pair(self, p1, p2, vs1, vs2): dist_centers = np.copy(p2.pos - p1.pos) else: dist_centers = p1.pos - p2.pos - expected_pos = system.part[rel_to].pos_folded + \ + expected_pos = system.part.by_id(rel_to).pos_folded + \ system.collision_detection.vs_placement * dist_centers dist = expected_pos - p.pos_folded dist -= np.round(dist / system.box_l) * system.box_l @@ -318,8 +318,8 @@ def test_bind_at_point_of_collision_random(self): for vs_pair in vs_pairs: # Get corresponding non-virtual particles base_particles = tuple(sorted( - [system.part[vs_pair[0]].vs_relative[0], - system.part[vs_pair[1]].vs_relative[0]])) + [system.part.by_id(vs_pair[0]).vs_relative[0], + system.part.by_id(vs_pair[1]).vs_relative[0]])) # Is there a corresponding bond? self.assertIn(base_particles, bonds) @@ -366,7 +366,7 @@ def run_test_glue_to_surface_for_pos(self, *positions): # Check that nothing explodes, when the particles are moved. # In particular for parallel simulations system.thermostat.set_langevin(kT=0, gamma=0.01, seed=42) - system.part[:].v = [0.05, 0.01, 0.15] + system.part.all().v = [0.05, 0.01, 0.15] system.integrator.run(3000) self.verify_state_after_glue_to_surface(expected_np) @@ -387,7 +387,7 @@ def verify_state_after_glue_to_surface(self, expected_np): self.assertEqual(p.bonds, ()) # Get base particles - base_p = system.part[p.vs_relative[0]] + base_p = system.part.by_id(p.vs_relative[0]) # Get bound particle # There is a bond between the base particle and the bound particle @@ -396,7 +396,7 @@ def verify_state_after_glue_to_surface(self, expected_np): p2 = None if len(base_p.bonds) == 1: self.assertEqual(base_p.bonds[0][0], self.H) - p2 = system.part[base_p.bonds[0][1]] + p2 = system.part.by_id(base_p.bonds[0][1]) else: # We need to go through all particles to find it for candidate in system.part: @@ -414,7 +414,7 @@ def verify_state_after_glue_to_surface(self, expected_np): self.verify_glue_to_surface_pair(base_p, p, p2) # Check particle that did not take part in collision. self.assertEqual(len(parts_not_accounted_for), 1) - p = system.part[parts_not_accounted_for[0]] + p = system.part.by_id(parts_not_accounted_for[0]) self.assertFalse(p.virtual) self.assertEqual(p.type, self.other_type) self.assertEqual(p.bonds, ()) @@ -546,16 +546,16 @@ def test_glue_to_surface_random(self): # Bonds to virtual sites: if bond[0] == self.H2: self.assertEqual( - system.part[bond[1]].type, + system.part.by_id(bond[1]).type, self.part_type_vs) else: self.assertEqual( - system.part[bond[1]].type, + system.part.by_id(bond[1]).type, self.part_type_to_attach_vs_to) elif p.type == self.part_type_to_attach_vs_to: self.assertEqual(bond[0], self.H) self.assertEqual( - system.part[bond[1]].type, + system.part.by_id(bond[1]).type, self.part_type_after_glueing) else: raise Exception( @@ -639,7 +639,7 @@ def test_bind_three_particles(self): system.integrator.run(1, recalc_forces=True) self.verify_triangle_binding(cutoff, system.bonded_inter[2], res) system.cell_system.set_n_square() - system.part[:].bonds = () + system.part.all().bonds = () system.integrator.run(1, recalc_forces=True) self.verify_triangle_binding(cutoff, system.bonded_inter[2], res) system.time_step = self.time_step @@ -653,7 +653,8 @@ def verify_triangle_binding(self, distance, first_bond, angle_res): expected_pairs = [] for i in range(n): for j in range(i + 1, n, 1): - if system.distance(system.part[i], system.part[j]) <= distance: + if system.distance(system.part.by_id(i), + system.part.by_id(j)) <= distance: expected_pairs.append((i, j)) # Find triangles @@ -664,9 +665,7 @@ def verify_triangle_binding(self, distance, first_bond, angle_res): for j in range(i + 1, n, 1): for k in range(j + 1, n, 1): # Ref to particles - p_i = system.part[i] - p_j = system.part[j] - p_k = system.part[k] + p_i, p_j, p_k = system.part.by_ids([i, j, k]) # Normalized distance vectors d_ij = np.copy(p_j.pos - p_i.pos) @@ -700,7 +699,7 @@ def verify_triangle_binding(self, distance, first_bond, angle_res): found_pairs = [] found_angle_bonds = [] for i in range(n): - for b in system.part[i].bonds: + for b in system.part.by_id(i).bonds: if len(b) == 2: self.assertEqual(b[0]._bond_id, self.H._bond_id) found_pairs.append(tuple(sorted((i, b[1])))) diff --git a/testsuite/python/comfixed.py b/testsuite/python/comfixed.py index 39695e82443..17427d84728 100644 --- a/testsuite/python/comfixed.py +++ b/testsuite/python/comfixed.py @@ -43,12 +43,13 @@ def test(self): v = 3 * [0.] # Make sure that id and type gaps work correctly system.part.add(id=2 * i, pos=r, v=v, type=2 * (i % 2)) + partcls = system.part.all() if espressomd.has_features(["MASS"]): # Avoid masses too small for the time step - system.part[:].mass = 2. * (0.1 + np.random.random(100)) + partcls.mass = 2. * (0.1 + np.random.random(100)) - com_0 = self.com(system.part[:]) + com_0 = self.com(partcls) system.comfixed.types = [0, 2] @@ -56,7 +57,7 @@ def test(self): self.assertEqual(system.comfixed.types, [2, 0]) for i in range(10): - com_i = self.com(system.part[:]) + com_i = self.com(partcls) for j in range(3): self.assertAlmostEqual(com_0[j], com_i[j], places=10) diff --git a/testsuite/python/constraint_shape_based.py b/testsuite/python/constraint_shape_based.py index 4818149c8d1..1f0163b1271 100644 --- a/testsuite/python/constraint_shape_based.py +++ b/testsuite/python/constraint_shape_based.py @@ -612,7 +612,7 @@ def test_spherocylinder(self): theta = 2. * i / float(N) * np.pi for r in radii: pos = r * np.array([np.cos(theta), 0, np.sin(theta)]) - system.part[0].pos = pos + self.box_l / 2.0 + p.pos = pos + self.box_l / 2.0 system.integrator.run(recalc_forces=True, steps=0) energy = system.analysis.energy() self.assertAlmostEqual(energy["total"], np.abs(10. - r)) diff --git a/testsuite/python/coulomb_cloud_wall.py b/testsuite/python/coulomb_cloud_wall.py index 904b65b00f2..c661660f837 100644 --- a/testsuite/python/coulomb_cloud_wall.py +++ b/testsuite/python/coulomb_cloud_wall.py @@ -62,7 +62,7 @@ def compare(self, method_name, energy=True, prefactor=None): # Compare forces and energy now in the system to stored ones # Force - force_diff = np.linalg.norm(self.system.part[:].f / prefactor - self.forces, + force_diff = np.linalg.norm(self.system.part.all().f / prefactor - self.forces, axis=1) self.assertLess( np.mean(force_diff), self.tolerance, diff --git a/testsuite/python/coulomb_cloud_wall_duplicated.py b/testsuite/python/coulomb_cloud_wall_duplicated.py index f5b778f446e..79289a0f360 100644 --- a/testsuite/python/coulomb_cloud_wall_duplicated.py +++ b/testsuite/python/coulomb_cloud_wall_duplicated.py @@ -57,7 +57,7 @@ def compare(self, method_name, energy=True): # Force force_diff = np.linalg.norm( - self.system.part[:].f - self.forces, axis=1) + self.system.part.all().f - self.forces, axis=1) self.assertLess( np.mean(force_diff), self.tolerance, msg="Absolute force difference too large for method " + method_name) diff --git a/testsuite/python/coulomb_mixed_periodicity.py b/testsuite/python/coulomb_mixed_periodicity.py index 1477a7c23a9..d5b276094ea 100644 --- a/testsuite/python/coulomb_mixed_periodicity.py +++ b/testsuite/python/coulomb_mixed_periodicity.py @@ -58,7 +58,7 @@ def compare(self, method_name, energy=True): # Force force_diff = np.linalg.norm( - self.system.part[:].f - self.forces, axis=1) + self.system.part.all().f - self.forces, axis=1) self.assertLessEqual( np.mean(force_diff), self.tolerance_force, "Absolute force difference too large for method " + method_name) diff --git a/testsuite/python/coulomb_tuning.py b/testsuite/python/coulomb_tuning.py index 50e99d5645c..737cb5c4d47 100644 --- a/testsuite/python/coulomb_tuning.py +++ b/testsuite/python/coulomb_tuning.py @@ -50,7 +50,7 @@ def tearDown(self): def compare(self, method_name): # Compare forces now in the system to stored ones difference = np.linalg.norm( - self.system.part[:].f - self.forces, axis=1) + self.system.part.all().f - self.forces, axis=1) self.assertLessEqual( np.mean(difference), self.tolerance, "Absolute force difference too large for method " + method_name) diff --git a/testsuite/python/dawaanr-and-bh-gpu.py b/testsuite/python/dawaanr-and-bh-gpu.py index f14a009fe80..1d0fb3d021d 100644 --- a/testsuite/python/dawaanr-and-bh-gpu.py +++ b/testsuite/python/dawaanr-and-bh-gpu.py @@ -79,8 +79,8 @@ def test(self): self.system.actors.add(dds_cpu) self.system.integrator.run(steps=0, recalc_forces=True) - dawaanr_f = np.copy(self.system.part[:].f) - dawaanr_t = np.copy(self.system.part[:].torque_lab) + dawaanr_f = np.copy(self.system.part.all().f) + dawaanr_t = np.copy(self.system.part.all().torque_lab) dawaanr_e = self.system.analysis.energy()["total"] del dds_cpu @@ -92,8 +92,8 @@ def test(self): self.system.actors.add(bh_gpu) self.system.integrator.run(steps=0, recalc_forces=True) - bhgpu_f = np.copy(self.system.part[:].f) - bhgpu_t = np.copy(self.system.part[:].torque_lab) + bhgpu_f = np.copy(self.system.part.all().f) + bhgpu_t = np.copy(self.system.part.all().torque_lab) bhgpu_e = self.system.analysis.energy()["total"] # compare diff --git a/testsuite/python/dawaanr-and-dds-gpu.py b/testsuite/python/dawaanr-and-dds-gpu.py index e1e17d1da4d..a4cb0b18fb8 100644 --- a/testsuite/python/dawaanr-and-dds-gpu.py +++ b/testsuite/python/dawaanr-and-dds-gpu.py @@ -75,8 +75,8 @@ def test(self): self.system.actors.add(dds_cpu) self.system.integrator.run(steps=0, recalc_forces=True) - dawaanr_f = np.copy(self.system.part[:].f) - dawaanr_t = np.copy(self.system.part[:].torque_lab) + dawaanr_f = np.copy(self.system.part.all().f) + dawaanr_t = np.copy(self.system.part.all().torque_lab) dawaanr_e = self.system.analysis.energy()["total"] del dds_cpu @@ -89,8 +89,8 @@ def test(self): self.system.actors.add(dds_gpu) self.system.integrator.run(steps=0, recalc_forces=True) - ddsgpu_f = np.copy(self.system.part[:].f) - ddsgpu_t = np.copy(self.system.part[:].torque_lab) + ddsgpu_f = np.copy(self.system.part.all().f) + ddsgpu_t = np.copy(self.system.part.all().torque_lab) ddsgpu_e = self.system.analysis.energy()["total"] # compare diff --git a/testsuite/python/dds-and-bh-gpu.py b/testsuite/python/dds-and-bh-gpu.py index 5a44d5c9851..ad6e463ef75 100644 --- a/testsuite/python/dds-and-bh-gpu.py +++ b/testsuite/python/dds-and-bh-gpu.py @@ -43,8 +43,8 @@ def test(self): pf_dds_gpu = 3.524 ratio_dawaanr_bh_gpu = pf_dds_gpu / pf_bh_gpu l = 15 - self.system.box_l = [l, l, l] - self.system.periodicity = [0, 0, 0] + self.system.box_l = 3 * [l] + self.system.periodicity = 3 * [False] self.system.time_step = 1E-4 self.system.cell_system.skin = 0.1 @@ -80,8 +80,8 @@ def test(self): self.system.actors.add(dds_gpu) self.system.integrator.run(steps=0, recalc_forces=True) - dawaanr_f = np.copy(self.system.part[:].f) - dawaanr_t = np.copy(self.system.part[:].torque_lab) + dawaanr_f = np.copy(self.system.part.all().f) + dawaanr_t = np.copy(self.system.part.all().torque_lab) dawaanr_e = self.system.analysis.energy()["total"] del dds_gpu @@ -93,8 +93,8 @@ def test(self): self.system.actors.add(bh_gpu) self.system.integrator.run(steps=0, recalc_forces=True) - bhgpu_f = np.copy(self.system.part[:].f) - bhgpu_t = np.copy(self.system.part[:].torque_lab) + bhgpu_f = np.copy(self.system.part.all().f) + bhgpu_t = np.copy(self.system.part.all().torque_lab) bhgpu_e = self.system.analysis.energy()["total"] # compare diff --git a/testsuite/python/dipolar_direct_summation.py b/testsuite/python/dipolar_direct_summation.py index 202d3dba68d..acb688ef84d 100644 --- a/testsuite/python/dipolar_direct_summation.py +++ b/testsuite/python/dipolar_direct_summation.py @@ -51,8 +51,8 @@ def dds_gpu_data(self): system.integrator.run(steps=0, recalc_forces=True) ref_e = system.analysis.energy()["dipolar"] - ref_f = np.copy(system.part[:].f) - ref_t = np.copy(system.part[:].torque_lab) + ref_f = np.copy(self.particles.f) + ref_t = np.copy(self.particles.torque_lab) system.actors.clear() @@ -66,8 +66,8 @@ def dds_data(self): system.integrator.run(steps=0, recalc_forces=True) ref_e = system.analysis.energy()["dipolar"] - ref_f = np.copy(system.part[:].f) - ref_t = np.copy(system.part[:].torque_lab) + ref_f = np.copy(self.particles.f) + ref_t = np.copy(self.particles.torque_lab) system.actors.clear() @@ -82,8 +82,8 @@ def dds_replica_data(self): system.integrator.run(steps=0, recalc_forces=True) ref_e = system.analysis.energy()["dipolar"] - ref_f = np.copy(system.part[:].f) - ref_t = np.copy(system.part[:].torque_lab) + ref_f = np.copy(self.particles.f) + ref_t = np.copy(self.particles.torque_lab) system.actors.clear() @@ -101,8 +101,8 @@ def fcs_data(self): system.integrator.run(0, recalc_forces=True) ref_e = system.analysis.energy()["dipolar"] - ref_f = np.copy(system.part[:].f) - ref_t = np.copy(system.part[:].torque_lab) + ref_f = np.copy(self.particles.f) + ref_t = np.copy(self.particles.torque_lab) system.actors.clear() @@ -131,8 +131,8 @@ def gen_reference_data(self, filepath_energy=OPEN_BOUNDARIES_REF_ENERGY, dipole_modulus = 1.3 part_pos = np.random.random((N, 3)) * system.box_l part_dip = dipole_modulus * tests_common.random_dipoles(N) - particles = system.part.add(pos=part_pos, dip=part_dip, - rotation=N * [(1, 1, 1)]) + self.particles = system.part.add(pos=part_pos, dip=part_dip, + rotation=N * [(1, 1, 1)]) # minimize system system.non_bonded_inter[0, 0].lennard_jones.set_params( @@ -154,8 +154,8 @@ def gen_reference_data(self, filepath_energy=OPEN_BOUNDARIES_REF_ENERGY, np.save( filepath_arrays, np.hstack( - (particles.pos_folded, - particles.dip, + (self.particles.pos_folded, + self.particles.dip, ref_f, ref_t)), allow_pickle=False) @@ -171,7 +171,8 @@ def check_open_bc(self, method_func=None, energy_tol=None, ref_f = array_data[:, 6:9] ref_t = array_data[:, 9:12] - system.part.add(pos=pos, dip=dip, rotation=[[1, 1, 1]] * len(pos)) + self.particles = system.part.add( + pos=pos, dip=dip, rotation=[[1, 1, 1]] * len(pos)) dds_e, dds_f, dds_t = method_func() self.assertAlmostEqual(dds_e, ref_e, delta=energy_tol) np.testing.assert_allclose(dds_f, ref_f, atol=force_tol) diff --git a/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py b/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py index ebfb0ffcce5..bdcca976a83 100644 --- a/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py +++ b/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py @@ -70,8 +70,8 @@ def test_mdlc(self): # Particles data = np.genfromtxt( tests_common.abspath("data/mdlc_reference_data_forces_torques.dat")) - s.part.add(pos=data[:, 1:4], dip=data[:, 4:7]) - s.part[:].rotation = (1, 1, 1) + partcls = s.part.add(pos=data[:, 1:4], dip=data[:, 4:7]) + partcls.rotation = 3 * [True] p3m = espressomd.magnetostatics.DipolarP3M( prefactor=DIPOLAR_PREFACTOR, mesh=32, accuracy=1E-4) @@ -81,9 +81,9 @@ def test_mdlc(self): s.actors.add(dlc) s.integrator.run(0) err_f = self.vector_error( - s.part[:].f, data[:, 7:10] * DIPOLAR_PREFACTOR) + partcls.f, data[:, 7:10] * DIPOLAR_PREFACTOR) err_t = self.vector_error( - s.part[:].torque_lab, data[:, 10:13] * DIPOLAR_PREFACTOR) + partcls.torque_lab, data[:, 10:13] * DIPOLAR_PREFACTOR) err_e = s.analysis.energy()["dipolar"] - ref_E tol_f = 2E-3 @@ -96,7 +96,8 @@ def test_mdlc(self): # Check if error is thrown when particles enter the MDLC gap # positive direction - s.part[0].pos = [ + p0 = s.part.by_id(0) + p0.pos = [ s.box_l[0] / 2, s.box_l[1] / 2, s.box_l[2] - gap_size / 2] @@ -105,7 +106,7 @@ def test_mdlc(self): with self.assertRaises(Exception): self.integrator.run(2) # negative direction - s.part[0].pos = [s.box_l[0] / 2, s.box_l[1] / 2, -gap_size / 2] + p0.pos = [s.box_l[0] / 2, s.box_l[1] / 2, -gap_size / 2] with self.assertRaises(Exception): self.system.analysis.energy() with self.assertRaises(Exception): @@ -126,8 +127,8 @@ def test_p3m(self): # Particles data = np.genfromtxt( tests_common.abspath("data/p3m_magnetostatics_system.data")) - s.part.add(pos=data[:, 1:4], dip=data[:, 4:7]) - s.part[:].rotation = (1, 1, 1) + partcls = s.part.add(pos=data[:, 1:4], dip=data[:, 4:7]) + partcls.rotation = 3 * [True] p3m = espressomd.magnetostatics.DipolarP3M( prefactor=DIPOLAR_PREFACTOR, mesh=32, accuracy=1E-6, epsilon="metallic") @@ -136,9 +137,9 @@ def test_p3m(self): expected = np.genfromtxt( tests_common.abspath("data/p3m_magnetostatics_expected.data"))[:, 1:] err_f = self.vector_error( - s.part[:].f, expected[:, 0:3] * DIPOLAR_PREFACTOR) + partcls.f, expected[:, 0:3] * DIPOLAR_PREFACTOR) err_t = self.vector_error( - s.part[:].torque_lab, expected[:, 3:6] * DIPOLAR_PREFACTOR) + partcls.torque_lab, expected[:, 3:6] * DIPOLAR_PREFACTOR) ref_E = 5.570 * DIPOLAR_PREFACTOR err_e = s.analysis.energy()["dipolar"] - ref_E @@ -166,8 +167,8 @@ def test_scafacos_dipoles(self): # Particles data = np.genfromtxt( tests_common.abspath("data/p3m_magnetostatics_system.data")) - s.part.add(pos=data[:, 1:4], dip=data[:, 4:7]) - s.part[:].rotation = (1, 1, 1) + partcls = s.part.add(pos=data[:, 1:4], dip=data[:, 4:7]) + partcls.rotation = 3 * [True] scafacos = espressomd.magnetostatics.Scafacos( prefactor=DIPOLAR_PREFACTOR, @@ -187,9 +188,9 @@ def test_scafacos_dipoles(self): expected = np.genfromtxt( tests_common.abspath("data/p3m_magnetostatics_expected.data"))[:, 1:] err_f = self.vector_error( - s.part[:].f, expected[:, 0:3] * DIPOLAR_PREFACTOR) + partcls.f, expected[:, 0:3] * DIPOLAR_PREFACTOR) err_t = self.vector_error( - s.part[:].torque_lab, expected[:, 3:6] * DIPOLAR_PREFACTOR) + partcls.torque_lab, expected[:, 3:6] * DIPOLAR_PREFACTOR) ref_E = 5.570 * DIPOLAR_PREFACTOR err_e = s.analysis.energy()["dipolar"] - ref_E diff --git a/testsuite/python/dipolar_p3m.py b/testsuite/python/dipolar_p3m.py index 65c3959b217..bca0e09c5a3 100644 --- a/testsuite/python/dipolar_p3m.py +++ b/testsuite/python/dipolar_p3m.py @@ -27,12 +27,11 @@ @utx.skipIfMissingFeatures(["DP3M"]) class MagnetostaticsP3M(ut.TestCase): - system = espressomd.System(box_l=[10., 10., 10.]) + system = espressomd.System(box_l=3 * [10.]) def setUp(self): - self.system.box_l = [10., 10., 10.] - self.system.part.add(pos=[4.0, 2.0, 2.0], dip=(1.3, 2.1, -6.0)) - self.system.part.add(pos=[6.0, 2.0, 2.0], dip=(7.3, 6.1, -4.0)) + self.partcls = self.system.part.add(pos=[[4.0, 2.0, 2.0], [6.0, 2.0, 2.0]], + dip=[(1.3, 2.1, -6.0), (7.3, 6.1, -4.0)]) def tearDown(self): self.system.part.clear() @@ -48,7 +47,7 @@ def test_dp3m(self): self.system.time_step = 0.01 prefactor = 1.1 box_vol = self.system.volume() - p1, p2 = self.system.part[:] + p1, p2 = self.partcls dip = np.copy(p1.dip + p2.dip) dp3m_params = {'accuracy': 1e-6, 'mesh': [49, 49, 49], @@ -82,7 +81,7 @@ def test_dp3m(self): # keep current values as reference to check for DP3M dipole correction ref_dp3m_energy_metallic = self.system.analysis.energy()['dipolar'] - ref_dp3m_forces_metallic = np.copy(self.system.part[:].f) + ref_dp3m_forces_metallic = np.copy(self.partcls.f) ref_dp3m_torque_metallic = np.array([ p1.convert_vector_space_to_body(p1.torque_lab), p2.convert_vector_space_to_body(p2.torque_lab)]) @@ -94,8 +93,8 @@ def test_dp3m(self): # keep current values as reference to check for MDLC dipole correction self.system.integrator.run(0, recalc_forces=True) ref_mdlc_energy_metallic = self.system.analysis.energy()['dipolar'] - ref_mdlc_forces_metallic = np.copy(self.system.part[:].f) - ref_mdlc_torque_metallic = np.copy(self.system.part[:].torque_lab) + ref_mdlc_forces_metallic = np.copy(self.partcls.f) + ref_mdlc_torque_metallic = np.copy(self.partcls.torque_lab) self.system.actors.clear() # check non-metallic case @@ -111,7 +110,7 @@ def test_dp3m(self): prefactor=prefactor, epsilon=epsilon, tune=False, **dp3m_params) self.system.actors.add(dp3m) self.system.integrator.run(0, recalc_forces=True) - dp3m_forces = np.copy(self.system.part[:].f) + dp3m_forces = np.copy(self.partcls.f) dp3m_torque = np.array([ p1.convert_vector_space_to_body(p1.torque_lab), p2.convert_vector_space_to_body(p2.torque_lab)]) @@ -127,8 +126,8 @@ def test_dp3m(self): mdlc = espressomd.magnetostatic_extensions.DLC(**mdlc_params) self.system.actors.add(mdlc) self.system.integrator.run(0, recalc_forces=True) - mdlc_forces = np.copy(self.system.part[:].f) - mdlc_torque = np.copy(self.system.part[:].torque_lab) + mdlc_forces = np.copy(self.partcls.f) + mdlc_torque = np.copy(self.partcls.torque_lab) mdlc_energy = self.system.analysis.energy()['dipolar'] np.testing.assert_allclose(mdlc_forces, ref_mdlc_forces, atol=tol) np.testing.assert_allclose(mdlc_torque, ref_mdlc_torque, atol=tol) diff --git a/testsuite/python/domain_decomposition.py b/testsuite/python/domain_decomposition.py index 44e07a44482..741a79ed13a 100644 --- a/testsuite/python/domain_decomposition.py +++ b/testsuite/python/domain_decomposition.py @@ -24,7 +24,7 @@ class DomainDecomposition(ut.TestCase): - system = espressomd.System(box_l=[50.0, 50.0, 50.0]) + system = espressomd.System(box_l=3 * [50.0]) original_node_grid = tuple(system.cell_system.node_grid) def setUp(self): @@ -40,10 +40,11 @@ def check_resort(self): n_part = 2351 # Add the particles on node 0, so that they have to be resorted - self.system.part.add(pos=n_part * [(0, 0, 0)], type=n_part * [1]) + partcls = self.system.part.add( + pos=n_part * [(0, 0, 0)], type=n_part * [1]) # And now change their positions - self.system.part[:].pos = self.system.box_l * \ + partcls.pos = self.system.box_l * \ np.random.random((n_part, 3)) # Add an interacting particle in a corner of the box @@ -63,7 +64,7 @@ def check_resort(self): # Check that we can still access all the particles # This basically checks if part_node and local_particles # is still in a valid state after the particle exchange - self.assertEqual(sum(self.system.part[:].type), n_part) + self.assertEqual(sum(self.system.part.all().type), n_part) # Check that the system is still valid if espressomd.has_features(['LENNARD_JONES']): diff --git a/testsuite/python/dpd.py b/testsuite/python/dpd.py index 64d74c4f23f..03911093574 100644 --- a/testsuite/python/dpd.py +++ b/testsuite/python/dpd.py @@ -125,7 +125,7 @@ def test_const_weight_function(self): system.integrator.run(0) # Outside of both cutoffs, forces should be 0 - for f in system.part[:].f: + for f in system.part.all().f: np.testing.assert_array_equal(np.copy(f), [0., 0., 0.]) # Only trans @@ -167,7 +167,7 @@ def calc_omega(dist, r_cut): system.integrator.run(0) # Outside of both cutoffs, forces should be 0 - for f in system.part[:].f: + for f in system.part.all().f: np.testing.assert_array_equal(np.copy(f), [0., 0., 0.]) # Only trans @@ -221,7 +221,7 @@ def calc_omega(dist, r_cut): p1 = system.part.add(pos=[3, 5, 5], type=0, v=v) # Outside of both cutoffs, forces should be 0 - for f in system.part[:].f: + for f in system.part.all().f: np.testing.assert_array_equal(np.copy(f), [0., 0., 0.]) # Place the particle at different positions to test the parabolic @@ -347,13 +347,13 @@ def calc_stress(dist, vel_diff): trans_weight_function=1, trans_gamma=gamma / 2.0, trans_r_cut=r_cut) pos = system.box_l * np.random.random((n_part, 3)) - system.part.add(pos=pos) + partcls = system.part.add(pos=pos) system.integrator.run(10) for kT in [0., 2.]: system.thermostat.set_dpd(kT=kT, seed=3) # run 1 integration step to get velocities - system.part[:].v = np.zeros((n_part, 3)) + partcls.v = np.zeros((n_part, 3)) system.integrator.run(steps=1) pairs = system.part.pairs() @@ -383,18 +383,18 @@ def test_momentum_conservation(self): system = self.system system.thermostat.set_dpd(kT=1.3, seed=42) - system.part.add(pos=((0, 0, 0), (0.1, 0.1, 0.1), (0.1, 0, 0)), - mass=(1, 2, 3)) + partcls = system.part.add(pos=((0, 0, 0), (0.1, 0.1, 0.1), (0.1, 0, 0)), + mass=(1, 2, 3)) system.non_bonded_inter[0, 0].dpd.set_params( weight_function=1, gamma=gamma, r_cut=r_cut, trans_weight_function=1, trans_gamma=gamma / 2.0, trans_r_cut=r_cut) - momentum = np.matmul(system.part[:].v.T, system.part[:].mass) + momentum = np.matmul(partcls.v.T, partcls.mass) for _ in range(10): system.integrator.run(25) - np.testing.assert_almost_equal(np.sum(system.part[:].f), 3 * [0.]) + np.testing.assert_almost_equal(np.sum(partcls.f), 3 * [0.]) np.testing.assert_allclose( - np.matmul(system.part[:].v.T, system.part[:].mass), momentum, atol=1E-12) + np.matmul(partcls.v.T, partcls.mass), momentum, atol=1E-12) if __name__ == "__main__": diff --git a/testsuite/python/dpd_stats.py b/testsuite/python/dpd_stats.py index 1c9cece51ed..c9d6571e700 100644 --- a/testsuite/python/dpd_stats.py +++ b/testsuite/python/dpd_stats.py @@ -42,14 +42,14 @@ def tearDown(self): self.system.integrator.set_vv() def check_total_zero(self): - v_total = np.sum(self.system.part[:].v, axis=0) + v_total = np.sum(self.system.part.all().v, axis=0) np.testing.assert_allclose(v_total, np.zeros(3), atol=1e-11) def single(self, with_langevin=False): """Test velocity distribution of a dpd fluid with a single type.""" N = 500 system = self.system - system.part.add(pos=system.box_l * np.random.random((N, 3))) + partcls = system.part.add(pos=system.box_l * np.random.random((N, 3))) kT = 2.3 gamma = 1.5 if with_langevin: @@ -63,7 +63,7 @@ def single(self, with_langevin=False): v_stored = np.zeros((loops, N, 3)) for i in range(loops): system.integrator.run(10) - v_stored[i] = system.part[:].v + v_stored[i] = partcls.v v_minmax = 5 bins = 5 error_tol = 0.01 @@ -104,7 +104,7 @@ def test_binary(self): v_stored = np.zeros((loops, N, 3)) for i in range(loops): system.integrator.run(10) - v_stored[i] = system.part[:].v + v_stored[i] = system.part.all().v v_minmax = 5 bins = 5 error_tol = 0.01 @@ -116,7 +116,7 @@ def test_disable(self): N = 200 system = self.system system.time_step = 0.01 - system.part.add(pos=system.box_l * np.random.random((N, 3))) + partcls = system.part.add(pos=system.box_l * np.random.random((N, 3))) kT = 2.3 gamma = 1.5 system.thermostat.set_dpd(kT=kT, seed=42) @@ -129,12 +129,12 @@ def test_disable(self): system.thermostat.turn_off() # Reset velocities - system.part[:].v = [1., 2., 3.] + partcls.v = [1., 2., 3.] system.integrator.run(10) # Check that there was neither noise nor friction - for v in system.part[:].v: + for v in partcls.v: for i in range(3): self.assertEqual(v[i], float(i + 1)) @@ -142,7 +142,7 @@ def test_disable(self): system.thermostat.set_dpd(kT=kT, seed=42) # Reset velocities for faster convergence - system.part[:].v = [0., 0., 0.] + partcls.v = [0., 0., 0.] # Equilibrate system.integrator.run(250) @@ -151,7 +151,7 @@ def test_disable(self): v_stored = np.zeros((loops, N, 3)) for i in range(loops): system.integrator.run(10) - v_stored[i] = system.part[:].v + v_stored[i] = partcls.v v_minmax = 5 bins = 5 error_tol = 0.012 diff --git a/testsuite/python/elc_vs_analytic.py b/testsuite/python/elc_vs_analytic.py index 8080f0cf0a7..7805e8a26d6 100644 --- a/testsuite/python/elc_vs_analytic.py +++ b/testsuite/python/elc_vs_analytic.py @@ -83,7 +83,7 @@ def test_elc(self): elc_results, analytic_results, rtol=0, atol=self.check_accuracy) def scan(self): - p1, p2 = self.system.part[:] + p1, p2 = self.system.part.all() result_array = np.empty((len(self.q), len(self.zPos), 2)) for chargeIndex, charge in enumerate(self.q): p1.q = charge diff --git a/testsuite/python/electrostaticInteractions.py b/testsuite/python/electrostaticInteractions.py index a24d98d3c15..e818bebfbd0 100644 --- a/testsuite/python/electrostaticInteractions.py +++ b/testsuite/python/electrostaticInteractions.py @@ -32,8 +32,8 @@ class ElectrostaticInteractionsTests(ut.TestCase): def setUp(self): self.system.time_step = 0.01 - self.system.part.add(pos=(9.0, 2.0, 2.0), q=1) - self.system.part.add(pos=(11.0, 2.0, 2.0), q=-1) + self.p0 = self.system.part.add(pos=(9.0, 2.0, 2.0), q=1) + self.p1 = self.system.part.add(pos=(11.0, 2.0, 2.0), q=-1) def tearDown(self): self.system.part.clear() @@ -41,7 +41,7 @@ def tearDown(self): def calc_dh_potential(self, r, dh_params): kT = 1.0 - q1, q2 = self.system.part[:].q + q1, q2 = self.system.part.all().q u = np.zeros_like(r) # r 1E-12: # sign undefined for phi = phi0 - force_axis = np.cross(self.system.distance_vec(p, p0), p.f) + force_axis = np.cross( + self.system.distance_vec( + p, self.p0), p.f) self.assertEqual( np.sign(phi_diff), np.sign(sign * np.dot(force_axis, self.axis)), @@ -122,7 +125,7 @@ def run_test(self, bond_instance, force_func, energy_func, phi0): # Total force =0? np.testing.assert_allclose( - np.sum(np.copy(self.system.part[:].f), 0), [0, 0, 0], atol=1E-12) + np.sum(np.copy(self.system.part.all().f), 0), [0, 0, 0], atol=1E-12) # No pressure (isotropic compression preserves angles) self.assertAlmostEqual( @@ -138,16 +141,16 @@ def run_test(self, bond_instance, force_func, energy_func, phi0): # and r_p2 =r_p1 +r_{p1,p2} and r_p3 =r_p1 +r_{p1,p3} # P_ij =1/V (F_p2 r_{p1,p2} + F_p3 r_{p1,p3}) p_tensor_expected = \ - np.outer(p1.f, self.system.distance_vec(p0, p1)) \ - + np.outer(p2.f, self.system.distance_vec(p0, p2)) + np.outer(self.p1.f, self.system.distance_vec(self.p0, self.p1)) \ + + np.outer(self.p2.f, self.system.distance_vec(self.p0, self.p2)) p_tensor_expected /= self.system.volume() np.testing.assert_allclose( self.system.analysis.pressure_tensor()["bonded"], p_tensor_expected, atol=1E-12) # Remove bonds - p0.delete_bond((bond_instance, 1, 2)) - p0.delete_bond((self.harmonic_bond, 1)) + self.p0.delete_bond((bond_instance, 1, 2)) + self.p0.delete_bond((self.harmonic_bond, 1)) def test_angle_harmonic(self): ah_bend = 1. diff --git a/testsuite/python/interactions_bonded.py b/testsuite/python/interactions_bonded.py index 0376b99e0e4..792a66d2e7c 100644 --- a/testsuite/python/interactions_bonded.py +++ b/testsuite/python/interactions_bonded.py @@ -83,7 +83,7 @@ def test_coulomb(self): coulomb_k = 1 q1 = 1 q2 = -1 - p1, p2 = self.system.part[:] + p1, p2 = self.system.part.all() p1.q = q1 p2.q = q2 self.run_test( @@ -98,7 +98,7 @@ def test_coulomb_sr(self): # interactions q1 = 1.2 q2 = -q1 - p1, p2 = self.system.part[:] + p1, p2 = self.system.part.all() p1.q = q1 p2.q = q2 r_cut = 2 @@ -143,7 +143,7 @@ def test_quartic(self): def run_test(self, bond_instance, force_func, energy_func, min_dist, cutoff, test_breakage=False): self.system.bonded_inter.add(bond_instance) - p1, p2 = self.system.part[:] + p1, p2 = self.system.part.all() p1.bonds = ((bond_instance, p2),) # n+1 steps from min_dist to cut, then we remove the cut, because that diff --git a/testsuite/python/interactions_dihedral.py b/testsuite/python/interactions_dihedral.py index 6a93f1dd17b..91c1cbb5c45 100644 --- a/testsuite/python/interactions_dihedral.py +++ b/testsuite/python/interactions_dihedral.py @@ -121,7 +121,7 @@ def dihedral_angle(self, p1, p2, p3, p4): # Test Dihedral Angle def test_dihedral(self): - p0, p1, p2, p3 = self.system.part[:] + p0, p1, p2, p3 = self.system.part.all() dh_k = 1 dh_phase = np.pi / 6 @@ -161,7 +161,7 @@ def test_dihedral(self): # Test Tabulated Dihedral Angle @utx.skipIfMissingFeatures(["TABULATED"]) def test_tabulated_dihedral(self): - p0, p1, p2, p3 = self.system.part[:] + p0, p1, p2, p3 = self.system.part.all() N = 111 d_phi = 2 * np.pi / N diff --git a/testsuite/python/interactions_non-bonded.py b/testsuite/python/interactions_non-bonded.py index bc52acb08d6..ebea83d6c0b 100644 --- a/testsuite/python/interactions_non-bonded.py +++ b/testsuite/python/interactions_non-bonded.py @@ -344,7 +344,7 @@ def get_reference_torque(gb_params, r, dir1, dir2): setup_system(gb_params) - p1, p2 = self.system.part[:] + p1, p2 = self.system.part.all() delta = 1.0e-6 @@ -405,7 +405,7 @@ def run_test(self, name, parameters, force_kernel, getattr(self.system.non_bonded_inter[0, 0], name).set_params( **parameters) - p0, p1 = self.system.part[:] + p0, p1 = self.system.part.all() p1.pos = p0.pos + self.step * n_initial_steps force_parameters = parameters.copy() diff --git a/testsuite/python/langevin_thermostat_stats.py b/testsuite/python/langevin_thermostat_stats.py index 117dcb10265..604ef46f7a9 100644 --- a/testsuite/python/langevin_thermostat_stats.py +++ b/testsuite/python/langevin_thermostat_stats.py @@ -192,7 +192,7 @@ def test_06__diffusion(self): # linear vel vel_obs = espressomd.observables.ParticleVelocities( - ids=system.part[:].id) + ids=system.part.all().id) corr_vel = espressomd.accumulators.Correlator( obs1=vel_obs, tau_lin=10, tau_max=1.4, delta_N=2, corr_operation="componentwise_product", compress1="discard1") @@ -200,7 +200,7 @@ def test_06__diffusion(self): # angular vel if espressomd.has_features("ROTATION"): omega_obs = espressomd.observables.ParticleBodyAngularVelocities( - ids=system.part[:].id) + ids=system.part.all().id) corr_omega = espressomd.accumulators.Correlator( obs1=omega_obs, tau_lin=10, tau_max=1.5, delta_N=2, corr_operation="componentwise_product", compress1="discard1") diff --git a/testsuite/python/lb_get_u_at_pos.py b/testsuite/python/lb_get_u_at_pos.py index 587a01d2cbb..8ae2abaf164 100644 --- a/testsuite/python/lb_get_u_at_pos.py +++ b/testsuite/python/lb_get_u_at_pos.py @@ -77,7 +77,7 @@ def test_get_u_at_pos(self): numpy.testing.assert_allclose( self.interpolated_vels[:-1], self.lb_fluid.get_interpolated_fluid_velocity_at_positions( - self.system.part[:].pos, False)[:-1], + self.system.part.all().pos, False)[:-1], atol=1e-4) diff --git a/testsuite/python/lb_stats.py b/testsuite/python/lb_stats.py index 549e963c715..00253ca0a1b 100644 --- a/testsuite/python/lb_stats.py +++ b/testsuite/python/lb_stats.py @@ -50,10 +50,10 @@ def tearDown(self): def test_mass_momentum_thermostat(self): self.n_col_part = 100 - self.system.part.add(pos=np.random.random( + partcls = self.system.part.add(pos=np.random.random( (self.n_col_part, 3)) * self.system.box_l[0]) if espressomd.has_features("MASS"): - self.system.part[:].mass = 0.1 + np.random.random( + partcls.mass = 0.1 + np.random.random( len(self.system.part)) self.system.thermostat.turn_off() diff --git a/testsuite/python/lb_thermostat.py b/testsuite/python/lb_thermostat.py index 66890536019..81846bb3d77 100644 --- a/testsuite/python/lb_thermostat.py +++ b/testsuite/python/lb_thermostat.py @@ -66,7 +66,7 @@ def test_velocity_distribution(self): v_stored = np.zeros((loops, N, 3)) for i in range(loops): self.system.integrator.run(3) - v_stored[i] = self.system.part[:].v + v_stored[i] = self.system.part.all().v minmax = 5 n_bins = 7 error_tol = 0.01 diff --git a/testsuite/python/lj.py b/testsuite/python/lj.py index f4891a71af0..9631323c22c 100644 --- a/testsuite/python/lj.py +++ b/testsuite/python/lj.py @@ -47,8 +47,9 @@ def setUp(self): self.system.part.add(pos=self.pos) def check(self): - f_diff = np.linalg.norm(self.system.part[:].f - self.forces, axis=1) - max_deviation = np.max(np.abs(self.system.part[:].f - self.forces)) + all_partcls = self.system.part.all() + f_diff = np.linalg.norm(all_partcls.f - self.forces, axis=1) + max_deviation = np.max(np.abs(all_partcls.f - self.forces)) self.assertLess(np.mean(f_diff), 1e-7) self.assertLess(max_deviation, 1e-5) diff --git a/testsuite/python/mass-and-rinertia_per_particle.py b/testsuite/python/mass-and-rinertia_per_particle.py index f84044da377..d3c1d7c5b34 100644 --- a/testsuite/python/mass-and-rinertia_per_particle.py +++ b/testsuite/python/mass-and-rinertia_per_particle.py @@ -151,13 +151,13 @@ def dissipation_param_setup(self, n): for i in range(n): for k in range(2): ind = i + k * n - self.system.part.add( + partcl = self.system.part.add( rotation=(1, 1, 1), pos=(0.0, 0.0, 0.0), id=ind) - self.system.part[ind].v = [1.0, 1.0, 1.0] + partcl.v = [1.0, 1.0, 1.0] if espressomd.has_features("ROTATION"): - self.system.part[ind].omega_body = [1.0, 1.0, 1.0] - self.system.part[ind].mass = self.mass - self.system.part[ind].rinertia = self.J + partcl.omega_body = [1.0, 1.0, 1.0] + partcl.mass = self.mass + partcl.rinertia = self.J def fluctuation_dissipation_param_setup(self, n): """ @@ -214,7 +214,7 @@ def fluctuation_dissipation_param_setup(self, n): part_pos = np.random.random(3) * box part_v = [0.0, 0.0, 0.0] part_omega_body = [0.0, 0.0, 0.0] - self.system.part.add( + partcl = self.system.part.add( rotation=(1, 1, 1), id=ind, mass=self.mass, @@ -222,7 +222,7 @@ def fluctuation_dissipation_param_setup(self, n): pos=part_pos, v=part_v) if espressomd.has_features("ROTATION"): - self.system.part[ind].omega_body = part_omega_body + partcl.omega_body = part_omega_body def check_dissipation(self, n): """ @@ -247,13 +247,13 @@ def check_dissipation(self, n): # Hence, only isotropic gamma_tran_p_validate could be # tested here. self.assertAlmostEqual( - self.system.part[ind].v[j], + self.system.part.by_id(ind).v[j], math.exp(- self.gamma_tran_p_validate[k, j] * self.system.time / self.mass), delta=tol) if espressomd.has_features("ROTATION"): self.assertAlmostEqual( - self.system.part[ind].omega_body[j], + self.system.part.by_id(ind).omega_body[j], math.exp(- self.gamma_rot_p_validate[k, j] * self.system.time / self.J[j]), delta=tol) @@ -288,7 +288,7 @@ def check_fluctuation_dissipation(self, n): for p in range(n): for k in range(2): ind = p + k * n - pos0[ind, :] = self.system.part[ind].pos + pos0[ind, :] = self.system.part.by_id(ind).pos dt0 = self.mass / self.gamma_tran_p_validate loops = 10 @@ -302,11 +302,12 @@ def check_fluctuation_dissipation(self, n): for p in range(n): for k in range(2): ind = p + k * n - v = self.system.part[ind].v + partcl = self.system.part.by_id(ind) + v = partcl.v if espressomd.has_features("ROTATION"): - o = self.system.part[ind].omega_body + o = partcl.omega_body o2[k, :] = o2[k, :] + np.power(o[:], 2) - pos = self.system.part[ind].pos + pos = partcl.pos v2[k, :] = v2[k, :] + np.power(v[:], 2) dr2[k, :] = np.power((pos[:] - pos0[ind, :]), 2) dt = (int_steps * (i + 1) + therm_steps) * \ @@ -380,9 +381,10 @@ def set_particle_specific_gamma(self, n): self.gamma_rot_p_validate[k, :] = self.gamma_rot_p[k, :] for i in range(n): ind = i + k * n - self.system.part[ind].gamma = self.gamma_tran_p[k, :] + partcl = self.system.part.by_id(ind) + partcl.gamma = self.gamma_tran_p[k, :] if espressomd.has_features("ROTATION"): - self.system.part[ind].gamma_rot = self.gamma_rot_p[k, :] + partcl.gamma_rot = self.gamma_rot_p[k, :] def set_diffusivity_tran(self): """ diff --git a/testsuite/python/mdanalysis.py b/testsuite/python/mdanalysis.py index 20648fa8ea4..8a2659f10c4 100644 --- a/testsuite/python/mdanalysis.py +++ b/testsuite/python/mdanalysis.py @@ -50,12 +50,13 @@ class TestMDAnalysis(ut.TestCase): system.bonded_inter.add(bond) system.bonded_inter.add(angle) system.bonded_inter.add(dihe) - system.part[1].add_bond((bond, 0)) - system.part[3].add_bond((angle, 2, 4)) - system.part[6].add_bond((dihe, 5, 7, 8)) + system.part.by_id(1).add_bond((bond, 0)) + system.part.by_id(3).add_bond((angle, 2, 4)) + system.part.by_id(6).add_bond((dihe, 5, 7, 8)) def test_universe(self): system = self.system + partcls = system.part.all() eos = espressomd.MDA_ESP.Stream(system) u = mda.Universe(eos.topology, eos.trajectory) # check atoms @@ -63,13 +64,13 @@ def test_universe(self): np.testing.assert_equal(u.atoms.ids, np.arange(10) + 1) np.testing.assert_equal(u.atoms.types, 5 * ['T0', 'T1']) np.testing.assert_almost_equal( - u.atoms.charges, system.part[:].q, decimal=6) + u.atoms.charges, partcls.q, decimal=6) np.testing.assert_almost_equal( - u.atoms.positions, system.part[:].pos, decimal=6) + u.atoms.positions, partcls.pos, decimal=6) np.testing.assert_almost_equal( - u.atoms.velocities, system.part[:].v, decimal=6) + u.atoms.velocities, partcls.v, decimal=6) np.testing.assert_almost_equal( - u.atoms.forces, system.part[:].f, decimal=6) + u.atoms.forces, partcls.f, decimal=6) # check bonds self.assertEqual(len(u.bonds), 1) self.assertEqual(len(u.angles), 1) diff --git a/testsuite/python/mmm1d.py b/testsuite/python/mmm1d.py index 45a6e529795..1940b88aec5 100644 --- a/testsuite/python/mmm1d.py +++ b/testsuite/python/mmm1d.py @@ -55,7 +55,7 @@ def tearDown(self): self.system.actors.clear() def test_forces(self): - measured_f = np.copy(self.system.part[:].f) + measured_f = np.copy(self.system.part.all().f) np.testing.assert_allclose(measured_f, self.forces_target, atol=self.allowed_error) diff --git a/testsuite/python/mpiio.py b/testsuite/python/mpiio.py index 6f9d1649621..1a2b1f34d1a 100644 --- a/testsuite/python/mpiio.py +++ b/testsuite/python/mpiio.py @@ -28,7 +28,7 @@ import unittest as ut import random import os -import argparse +import dataclasses # Number of particles npart = 10 @@ -40,6 +40,15 @@ filenames = [filename + "." + ext for ext in exts] +@dataclasses.dataclass +class MockParticle: + id: int + type: int + pos: np.ndarray + v: np.ndarray + bonds: list + + def clean_files(): for f in filenames: if os.path.isfile(f): @@ -54,16 +63,15 @@ def randint_different_from(a, b, n): return r -def random_particles(): +def get_random_mock_particles(): """Returns a list of random particle descriptions.""" parts = [] for i in range(npart): - p = argparse.Namespace() - p.id = i - p.type = random.randint(0, 100) - p.pos = np.random.rand(3) - p.v = np.random.rand(3) - p.bonds = [] + p = MockParticle(id=i, + type=random.randint(0, 100), + pos=np.random.rand(3), + v=np.random.rand(3), + bonds=[]) # Up to 20 bonds; otherwise this test will take ages for _ in range(random.randint(0, 20)): btype = random.randint(0, nbonds - 1) @@ -90,16 +98,16 @@ class MPIIOTest(ut.TestCase): for i in range(nbonds): system.bonded_inter[i] = espressomd.interactions.AngleHarmonic( bend=i, phi0=i) - test_particles = random_particles() + test_mock_particles = get_random_mock_particles() def setUp(self): - """Sets up a system from test_particles and prepares environment + """Sets up a system from test_mock_particles and prepares environment for the tests.""" clean_files() # Prior call might not have completed successfully - for p in self.test_particles: + for p in self.test_mock_particles: self.system.part.add(id=p.id, type=p.type, pos=p.pos, v=p.v) for b in p.bonds: - self.system.part[p.id].add_bond(b) + self.system.part.by_id(p.id).add_bond(b) def tearDown(self): clean_files() @@ -112,7 +120,7 @@ def check_files_exist(self): def check_sample_system(self): """Checks the particles in the ESPResSo system "self.s" against the true values in "self.test_particles".""" - for p, q in zip(self.system.part, self.test_particles): + for p, q in zip(self.system.part, self.test_mock_particles): self.assertEqual(p.id, q.id) self.assertEqual(p.type, q.type) np.testing.assert_array_equal(np.copy(p.pos), q.pos) diff --git a/testsuite/python/nsquare.py b/testsuite/python/nsquare.py index 463132422e5..7d5775d66c5 100644 --- a/testsuite/python/nsquare.py +++ b/testsuite/python/nsquare.py @@ -40,10 +40,11 @@ def test_load_balancing(self): n_part_avg = n_part // n_nodes # Add the particles on node 0, so that they have to be resorted - self.system.part.add(pos=n_part * [(0, 0, 0)], type=n_part * [1]) + partcls = self.system.part.add( + pos=n_part * [(0, 0, 0)], type=n_part * [1]) # And now change their positions - self.system.part[:].pos = self.system.box_l * \ + partcls.pos = self.system.box_l * \ np.random.random((n_part, 3)) # Add an interacting particle in a corner of the box @@ -67,7 +68,7 @@ def test_load_balancing(self): # Check that we can still access all the particles # This basically checks if part_node and local_particles # are still in a valid state after the particle exchange - self.assertEqual(sum(self.system.part[:].type), n_part) + self.assertEqual(sum(partcls.type), n_part) # Check that the system is still valid if espressomd.has_features(['LENNARD_JONES']): diff --git a/testsuite/python/observable_chain.py b/testsuite/python/observable_chain.py index bb199e03035..83c9331df10 100644 --- a/testsuite/python/observable_chain.py +++ b/testsuite/python/observable_chain.py @@ -51,6 +51,7 @@ class ObservableTests(ut.TestCase): def setUp(self): for i in range(self.n_parts): self.system.part.add(pos=[1 + i, 1 + i, 1 + i], id=i) + self.partcls = self.system.part.all() def tearDown(self): self.system.part.clear() @@ -74,7 +75,7 @@ def test_ParticleDistances(self): for i in range(1, self.n_parts): pos[i] = pos[i - 1] + np.random.uniform( low=0, high=max_bond_length, size=3) - self.system.part[:].pos = pos + self.partcls.pos = pos # expected values distances = np.linalg.norm(pos[1:] - pos[:-1], axis=1) # observed values @@ -113,7 +114,7 @@ def test_BondAngles(self): for i in range(1, self.n_parts): pos[i] = pos[i - 1] + np.random.uniform( low=0, high=max_bond_length, size=3) - self.system.part[:].pos = pos + self.partcls.pos = pos # expected values v1 = pos[:-2] - pos[1:-1] v2 = pos[2:] - pos[1:-1] @@ -175,7 +176,7 @@ def place_particles(bl, offset): pos[3] = pos[2] + [bl * np.cos(np.pi - phi), bl * np.sin(phi), 0.] pos[4] = pos[3] + [bl, 0., 0.] pos += offset - self.system.part[:].pos = pos + self.partcls.pos = pos return pos pids = list(range(self.n_parts)) @@ -183,17 +184,17 @@ def place_particles(bl, offset): obs_chain = espressomd.observables.BondDihedrals(ids=pids) # test multiple angles, take periodic boundaries into account - p = self.system.part + p0, p4 = self.system.part.by_ids([0, 4]) for bond_length in [0.1, self.box_l / 2.0]: for offset in [1.0, self.box_l / 2.0]: for phi in np.arange(0, np.pi, np.pi / 6): # place particles and keep list of unfolded positions pos = place_particles(bond_length, 3 * [offset]) # rotate the 1st particle - p[0].pos = pos[0] = rotate_particle(*pos[1:4, :][::-1], - phi=phi) + p0.pos = pos[0] = rotate_particle(*pos[1:4, :][::-1], + phi=phi) # rotate the 5th particle - p[4].pos = pos[4] = rotate_particle(*pos[2:5, :], phi=phi) + p4.pos = pos[4] = rotate_particle(*pos[2:5, :], phi=phi) # expected values dih1 = calculate_dihedral(*pos[0:4, :][::-1]) dih2 = calculate_dihedral(*pos[1:5, :]) @@ -219,12 +220,12 @@ def place_particles(bl, offset): def test_CosPersistenceAngles(self): # First test: compare with python implementation self.system.part.clear() - self.system.part.add(pos=np.array( + partcls = self.system.part.add(pos=np.array( [np.linspace(0, self.system.box_l[0], 20)] * 3).T + np.random.random((20, 3))) obs = espressomd.observables.CosPersistenceAngles( - ids=self.system.part[:].id) + ids=partcls.id) np.testing.assert_allclose( - obs.calculate(), cos_persistence_angles(self.system.part[:].pos)) + obs.calculate(), cos_persistence_angles(partcls.pos)) self.system.part.clear() # Second test: place particles with fixed angles and check that the # result of PersistenceAngle.calculate()[i] is i*phi @@ -232,8 +233,9 @@ def test_CosPersistenceAngles(self): for i in range(10): pos = [np.cos(i * delta_phi), np.sin(i * delta_phi), 0.0] self.system.part.add(pos=pos) + new_partcls = self.system.part.all() obs = espressomd.observables.CosPersistenceAngles( - ids=self.system.part[:].id) + ids=new_partcls.id) expected = np.arange(1, 9) * delta_phi np.testing.assert_allclose(obs.calculate(), np.cos(expected)) diff --git a/testsuite/python/observable_cylindrical.py b/testsuite/python/observable_cylindrical.py index 8cf7c14f4fd..f1d12d37893 100644 --- a/testsuite/python/observable_cylindrical.py +++ b/testsuite/python/observable_cylindrical.py @@ -134,7 +134,7 @@ def setup_system_get_np_hist(self): self.cyl_transform_params.center) vel_aligned.append(self.align_with_observable_frame(vel)) self.system.part.add(pos=pos_aligned, v=vel_aligned) - self.params['ids'] = self.system.part[:].id + self.params['ids'] = self.system.part.all().id return np_dens, np_edges @@ -206,7 +206,7 @@ def test_cylindrical_pid_profile_interface(self): params['n_z_bins'] = 8 self.system.part.add(pos=[0, 0, 0], type=0) self.system.part.add(pos=[0, 0, 0], type=1) - params['ids'] = self.system.part[:].id + params['ids'] = self.system.part.all().id observable = espressomd.observables.CylindricalDensityProfile(**params) # check pids np.testing.assert_array_equal(np.copy(observable.ids), params['ids']) diff --git a/testsuite/python/observable_cylindricalLB.py b/testsuite/python/observable_cylindricalLB.py index 0b77f86e049..5c2579cf69b 100644 --- a/testsuite/python/observable_cylindricalLB.py +++ b/testsuite/python/observable_cylindricalLB.py @@ -138,7 +138,7 @@ def setup_system_get_np_hist(self): [0.5]), dtype=int) self.system.part.add(pos=pos_aligned, v=vel_aligned) - self.params['ids'] = self.system.part[:].id + self.params['ids'] = self.system.part.all().id for node, vel in zip(node_aligned, vel_aligned): self.lbf[node].velocity = vel diff --git a/testsuite/python/observable_profile.py b/testsuite/python/observable_profile.py index 2a6baa11fc3..37e432e45b6 100644 --- a/testsuite/python/observable_profile.py +++ b/testsuite/python/observable_profile.py @@ -29,7 +29,7 @@ class ProfileObservablesTest(ut.TestCase): system.part.add(pos=[4.0, 4.0, 6.0], v=[0.0, 0.0, 1.0]) system.part.add(pos=[4.0, 4.0, 6.0], v=[0.0, 0.0, 1.0]) bin_volume = 5.0**3 - kwargs = {'ids': list(system.part[:].id), + kwargs = {'ids': list(system.part.all().id), 'n_x_bins': 2, 'n_y_bins': 3, 'n_z_bins': 4, @@ -47,7 +47,7 @@ def test_density_profile(self): obs_bin_edges = density_profile.bin_edges() obs_bin_centers = density_profile.bin_centers() np_hist, np_edges = tests_common.get_histogram( - np.copy(self.system.part[:].pos), self.kwargs, 'cartesian', + np.copy(self.system.part.all().pos), self.kwargs, 'cartesian', normed=True) np_hist *= len(self.system.part) np.testing.assert_array_almost_equal(obs_data, np_hist) @@ -69,7 +69,7 @@ def test_density_profile(self): def test_force_density_profile(self): density_profile = espressomd.observables.ForceDensityProfile( **self.kwargs) - self.system.part[:].ext_force = [0.0, 0.0, 1.0] + self.system.part.all().ext_force = [0.0, 0.0, 1.0] self.system.integrator.run(0) obs_data = density_profile.calculate() np.testing.assert_array_equal( @@ -91,7 +91,7 @@ def test_flux_density_profile(self): def test_pid_profile_interface(self): # test setters and getters - params = {'ids': list(self.system.part[:].id), + params = {'ids': list(self.system.part.all().id), 'n_x_bins': 4, 'n_y_bins': 6, 'n_z_bins': 8, diff --git a/testsuite/python/observables.py b/testsuite/python/observables.py index dd52a2c10e6..3ad7974201e 100644 --- a/testsuite/python/observables.py +++ b/testsuite/python/observables.py @@ -24,44 +24,45 @@ def calc_com_x(system, x, id_list): """Mass-weighted average, skipping virtual sites""" - masses = system.part[id_list].mass + partcls = system.part.by_ids(id_list) + masses = partcls.mass # Filter out virtual particles by using mass=0 for them - virtual = system.part[id_list].virtual + virtual = partcls.virtual masses[np.nonzero(virtual)] = 0. - return np.average(getattr(system.part[id_list], x), weights=masses, axis=0) + return np.average(getattr(partcls, x), weights=masses, axis=0) class Observables(ut.TestCase): N_PART = 200 # Handle for espresso system system = espressomd.System(box_l=[10.0, 10.0, 10.0]) - system.part.add( + partcls = 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) + partcls.mass = np.random.random(N_PART) if espressomd.has_features(["DIPOLES"]): - system.part[:].dip = np.random.random((N_PART, 3)) - .3 + partcls.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)) + partcls.omega_body = np.random.random((N_PART, 3)) - .5 + partcls.torque_lab = np.random.random((N_PART, 3)) - .5 + partcls.quat = np.random.random((N_PART, 4)) if espressomd.has_features("DIPOLES"): - system.part[:].dipm = np.random.random(N_PART) + 2 + partcls.dipm = np.random.random(N_PART) + 2 if espressomd.has_features("ELECTROSTATICS"): - system.part[:].q = np.random.random(N_PART) + partcls.q = np.random.random(N_PART) if espressomd.has_features("VIRTUAL_SITES"): - p = system.part[system.part[:].id[8]] + p = system.part.by_id(partcls.id[8]) p.virtual = True def generate_test_for_pid_observable( @@ -80,7 +81,7 @@ def func(self): # Randomly pick a subset of the particles id_list = sorted( np.random.choice( - self.system.part[:].id, + self.system.part.all().id, size=int( self.N_PART * .9), replace=False)) @@ -90,10 +91,10 @@ def func(self): # Get data from particles if pprop_name == "f": for p_id in id_list: - if self.system.part[p_id].virtual: + if self.system.part.by_id(p_id).virtual: id_list.remove(p_id) - part_data = getattr(self.system.part[id_list], pprop_name) + part_data = getattr(self.system.part.by_ids(id_list), pprop_name) # Reshape and aggregate to linear array if len(part_data.shape) > 1: @@ -147,7 +148,7 @@ def func(self): @utx.skipIfMissingFeatures(['ROTATION']) def test_particle_body_velocities(self): obs = espressomd.observables.ParticleBodyVelocities( - ids=self.system.part[:].id) + ids=self.system.part.all().id) obs_data = obs.calculate() part_data = np.array([p.convert_vector_space_to_body(p.v) for p in self.system.part]) @@ -188,9 +189,9 @@ def test_pressure_tensor(self): @utx.skipIfMissingFeatures('ELECTROSTATICS') def test_dipolemoment(self): - obs = espressomd.observables.DipoleMoment(ids=self.system.part[:].id) + obs = espressomd.observables.DipoleMoment(ids=self.partcls.id) obs_data = obs.calculate() - part_data = self.system.part[:].q.dot(self.system.part[:].pos) + part_data = self.partcls.q.dot(self.partcls.pos) self.assertEqual(obs_data.shape, part_data.shape) np.testing.assert_array_almost_equal( obs_data, part_data, err_msg="Data did not agree for observable 'DipoleMoment'", decimal=9) @@ -198,7 +199,7 @@ def test_dipolemoment(self): def test_com_force(self): id_list = sorted( np.random.choice( - self.system.part[:].id, + self.partcls.id, size=int( self.N_PART * .9), replace=False)) diff --git a/testsuite/python/oif_volume_conservation.py b/testsuite/python/oif_volume_conservation.py index 0d7cc332e92..7169b718280 100644 --- a/testsuite/python/oif_volume_conservation.py +++ b/testsuite/python/oif_volume_conservation.py @@ -46,27 +46,26 @@ def test(self): cell0 = oif.OifCell( cell_type=cell_type, particle_type=0, origin=[5.0, 5.0, 5.0]) self.assertEqual(system.max_oif_objects, 1) - # cell0.output_vtk_pos_folded(file_name="cell0_0.vtk") + partcls = system.part.all() # fluid - diameter_init = cell0.diameter() print(f"initial diameter = {diameter_init}") # OIF object is being stretched by factor 1.5 - system.part[:].pos = (system.part[:].pos - 5) * 1.5 + 5 + partcls.pos = (partcls.pos - 5) * 1.5 + 5 diameter_stretched = cell0.diameter() print(f"stretched diameter = {diameter_stretched}") # Apply non-isotropic deformation - system.part[:].pos = system.part[:].pos * np.array((0.96, 1.05, 1.02)) + partcls.pos = partcls.pos * np.array((0.96, 1.05, 1.02)) # Test that restoring forces net to zero and don't produce a torque system.integrator.run(1) np.testing.assert_allclose( np.sum( - system.part[:].f, axis=0), [ + partcls.f, axis=0), [ 0., 0., 0.], atol=1E-12) total_torque = np.zeros(3) diff --git a/testsuite/python/p3m_tuning_exceptions.py b/testsuite/python/p3m_tuning_exceptions.py index e6134d87b1e..52781405816 100644 --- a/testsuite/python/p3m_tuning_exceptions.py +++ b/testsuite/python/p3m_tuning_exceptions.py @@ -242,7 +242,7 @@ def test_04_invalid_params_p3m_elc_cpu(self): with self.assertRaisesRegex(ValueError, "maxPWerror must be > 0"): self.system.actors.add(solver_elc) - self.system.part[0].q = -2 + self.system.part.by_id(0).q = -2 self.system.actors.clear() solver_elc = espressomd.electrostatics.ELC( @@ -258,7 +258,7 @@ def test_04_invalid_params_p3m_elc_cpu(self): with self.assertRaisesRegex(RuntimeError, "ELC does not work for non-neutral systems and non-metallic dielectric contrast"): self.system.actors.add(solver_elc) - self.system.part[0].q = -1 + self.system.part.by_id(0).q = -1 self.system.actors.clear() solver_elc = espressomd.electrostatics.ELC( diff --git a/testsuite/python/pair_criteria.py b/testsuite/python/pair_criteria.py index 1e7c1fa3d88..95790e0d2a4 100644 --- a/testsuite/python/pair_criteria.py +++ b/testsuite/python/pair_criteria.py @@ -45,7 +45,7 @@ def test_distance_crit_periodic(self): # Decisions # Periodic system. Particles in range via minimum image convention - self.system.periodicity = (1, 1, 1) + self.system.periodicity = 3 * [True] self.assertTrue(dc.decide(self.p1, self.p2)) self.assertTrue(dc.decide(self.p1.id, self.p2.id)) @@ -101,13 +101,13 @@ def test_bond_crit(self): self.assertTrue(not bc.decide(self.p1.id, self.p2.id)) # Add bond. Then the criterion should match - self.system.part[self.p1.id].bonds = ((bt, self.p2.id),) + self.p1.bonds = ((bt, self.p2.id),) self.assertTrue(bc.decide(self.p1, self.p2)) self.assertTrue(bc.decide(self.p1.id, self.p2.id)) # Place bond on the 2nd particle. The criterion should still match - self.system.part[self.p1.id].bonds = () - self.system.part[self.p2.id].bonds = ((bt, self.p1.id),) + self.p1.bonds = () + self.p2.bonds = ((bt, self.p1.id),) self.assertTrue(bc.decide(self.p1, self.p2)) self.assertTrue(bc.decide(self.p1.id, self.p2.id)) diff --git a/testsuite/python/particle.py b/testsuite/python/particle.py index ae3ddda75e5..d68e391d3ae 100644 --- a/testsuite/python/particle.py +++ b/testsuite/python/particle.py @@ -43,7 +43,7 @@ class ParticleProperties(ut.TestCase): def setUp(self): if not self.system.part.exists(self.pid): - self.system.part.add(id=self.pid, pos=(0, 0, 0)) + self.partcl = self.system.part.add(id=self.pid, pos=(0, 0, 0)) def tearDown(self): self.system.part.clear() @@ -62,9 +62,9 @@ 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 - setattr(self.system.part[self.pid], propName, value) + setattr(self.partcl, propName, value) np.testing.assert_allclose( - np.array(getattr(self.system.part[self.pid], propName)), value, + np.array(getattr(self.partcl, propName)), value, err_msg=propName + ": value set and value gotten back differ.", atol=self.tol) @@ -84,8 +84,8 @@ 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 - setattr(self.system.part[self.pid], propName, value) - self.assertEqual(getattr(self.system.part[self.pid], propName), + setattr(self.partcl, propName, value) + self.assertEqual(getattr(self.partcl, propName), value, propName + ": value set and value gotten back differ.") return func @@ -129,9 +129,9 @@ def test_director(self): sample_vector_normalized = sample_vector / \ np.linalg.norm(sample_vector) - setattr(self.system.part[self.pid], "director", sample_vector) + setattr(self.partcl, "director", sample_vector) np.testing.assert_allclose( - np.array(getattr(self.system.part[self.pid], "director")), + np.array(getattr(self.partcl, "director")), sample_vector_normalized ) @@ -141,9 +141,9 @@ def test_director(self): "gamma", np.array([2., 9., 0.23])) def test_gamma_single(self): - self.system.part[self.pid].gamma = 17.4 + self.partcl.gamma = 17.4 np.testing.assert_array_equal( - np.copy(self.system.part[self.pid].gamma), + np.copy(self.partcl.gamma), np.array([17.4, 17.4, 17.4]), "gamma: value set and value gotten back differ.") else: @@ -154,9 +154,9 @@ def test_gamma_single(self): "gamma_rot", np.array([5., 10., 0.33])) def test_gamma_rot_single(self): - self.system.part[self.pid].gamma_rot = 15.4 + self.partcl.gamma_rot = 15.4 np.testing.assert_array_equal( - np.copy(self.system.part[self.pid].gamma_rot), + np.copy(self.partcl.gamma_rot), np.array([15.4, 15.4, 15.4]), "gamma_rot: value set and value gotten back differ.") else: @@ -185,12 +185,12 @@ def test_gamma_rot_single(self): if espressomd.has_features(["VIRTUAL_SITES_RELATIVE"]): def test_yy_vs_relative(self): self.system.part.add(id=0, pos=(0, 0, 0)) - self.system.part.add(id=1, pos=(0, 0, 0)) - self.system.part[1].vs_relative = (0, 5.0, (0.5, -0.5, -0.5, -0.5)) - self.system.part[1].vs_quat = [1, 2, 3, 4] + p1 = self.system.part.add(id=1, pos=(0, 0, 0)) + p1.vs_relative = (0, 5.0, (0.5, -0.5, -0.5, -0.5)) + p1.vs_quat = [1, 2, 3, 4] np.testing.assert_array_equal( - self.system.part[1].vs_quat, [1, 2, 3, 4]) - res = self.system.part[1].vs_relative + p1.vs_quat, [1, 2, 3, 4]) + res = p1.vs_relative self.assertEqual(res[0], 0, "vs_relative: " + res.__str__()) self.assertEqual(res[1], 5.0, "vs_relative: " + res.__str__()) np.testing.assert_allclose( @@ -198,13 +198,13 @@ def test_yy_vs_relative(self): err_msg="vs_relative: " + res.__str__(), atol=self.tol) # check exceptions with self.assertRaisesRegex(ValueError, "needs input in the form"): - self.system.part[1].vs_relative = (0, 5.0) + p1.vs_relative = (0, 5.0) with self.assertRaisesRegex(ValueError, "particle id has to be given as an int"): - self.system.part[1].vs_relative = ('0', 5.0, (1, 0, 0, 0)) + p1.vs_relative = ('0', 5.0, (1, 0, 0, 0)) with self.assertRaisesRegex(ValueError, "distance has to be given as a float"): - self.system.part[1].vs_relative = (0, '5', (1, 0, 0, 0)) + p1.vs_relative = (0, '5', (1, 0, 0, 0)) with self.assertRaisesRegex(ValueError, "quaternion has to be given as a tuple of 4 floats"): - self.system.part[1].vs_relative = (0, 5.0, (1, 0, 0)) + p1.vs_relative = (0, 5.0, (1, 0, 0)) @utx.skipIfMissingFeatures("DIPOLES") def test_contradicting_properties_dip_dipm(self): @@ -260,20 +260,20 @@ def test_image_box(self): pos = 1.5 * s.box_l - s.part.add(pos=pos) + p = s.part.add(pos=pos) - np.testing.assert_equal(np.copy(s.part[0].image_box), [1, 1, 1]) + np.testing.assert_equal(np.copy(p.image_box), [1, 1, 1]) def test_accessing_invalid_id_raises(self): self.system.part.clear() - handle_to_non_existing_particle = self.system.part[42] + handle_to_non_existing_particle = self.system.part.by_id(42) with self.assertRaises(RuntimeError): handle_to_non_existing_particle.id def test_parallel_property_setters(self): s = self.system s.part.clear() - s.part.add(pos=s.box_l * np.random.random((100, 3))) + partcls = s.part.add(pos=s.box_l * np.random.random((100, 3))) # Copy individual properties of particle 0 print( @@ -282,21 +282,21 @@ def test_parallel_property_setters(self): # Uncomment to identify guilty property # print( p) - assert hasattr(s.part[0], p), \ + assert hasattr(s.part.by_id(0), p), \ "Inconsistency between ParticleHandle and particle_data.particle_attributes" try: - setattr(s.part[:], p, getattr(s.part[0], p)) + setattr(partcls, p, getattr(s.part.by_id(0), p)) except AttributeError: print("Skipping read-only", p) # Cause a different mpi callback to uncover deadlock immediately - _ = getattr(s.part[:], p) + _ = getattr(partcls, p) def test_remove_particle(self): """Tests that if a particle is removed, it no longer exists and bonds to the removed particle are also removed.""" - p1 = self.system.part[self.pid] + p1 = self.system.part.by_id(self.pid) p2 = self.system.part.add(pos=p1.pos, bonds=[(self.f1, p1.id)]) p1.remove() @@ -306,7 +306,7 @@ def test_remove_particle(self): def test_bonds(self): """Tests bond addition and removal.""" - p1 = self.system.part[self.pid] + p1 = self.system.part.by_id(self.pid) p2 = self.system.part.add(pos=p1.pos) inactive_bond = espressomd.interactions.FeneBond(k=1, d_r_max=2) p2.add_bond([self.f1, p1]) @@ -329,17 +329,17 @@ def test_bonds(self): p2.delete_bond((self.f1, 'p1')) def test_zz_remove_all(self): - for pid in self.system.part[:].id: - self.system.part[pid].remove() + for p in self.system.part.all(): + p.remove() self.system.part.add( pos=np.random.random((100, 3)) * self.system.box_l, id=np.arange(100, dtype=int)) - ids = self.system.part[:].id + ids = self.system.part.all().id np.random.shuffle(ids) for pid in ids: - self.system.part[pid].remove() + self.system.part.by_id(pid).remove() with self.assertRaises(Exception): - self.system.part[17].remove() + self.system.part.by_id(17).remove() def test_coord_fold_corner_cases(self): system = self.system @@ -381,11 +381,12 @@ def test_particle_slice(self): # Empty slice system.part.clear() - self.assertEqual(len(system.part[:]), 0) - self.assertEqual(len(system.part[:].pos), 0) - self.assertEqual(len(system.part[:].id), 0) + all_partcls_empty = system.part.all() + self.assertEqual(len(all_partcls_empty), 0) + self.assertEqual(len(all_partcls_empty), 0) + self.assertEqual(len(all_partcls_empty), 0) with self.assertRaises(AttributeError): - system.part[:].pos = ((1, 2, 3,),) + all_partcls_empty.pos = ((1, 2, 3,),) # Slice containing particles ids = [1, 4, 6, 3, 8, 9] @@ -393,34 +394,30 @@ def test_particle_slice(self): system.part.add(id=ids, pos=pos) # All particles - self.assertEqual(len(system.part[:]), len(ids)) - np.testing.assert_equal(system.part[:].id, sorted(ids)) - np.testing.assert_equal(system.part[:].pos, pos[np.argsort(ids)]) + all_partcls = system.part.all() + self.assertEqual(len(all_partcls), len(ids)) + np.testing.assert_equal(all_partcls.id, sorted(ids)) + np.testing.assert_equal(all_partcls.pos, pos[np.argsort(ids)]) # Access via slicing - np.testing.assert_equal(system.part[4:9].id, + np.testing.assert_equal(system.part.by_ids(range(4, 9)).id, [i for i in sorted(ids) if i >= 4 and i < 9]) - np.testing.assert_equal(system.part[9:4:-1].id, + np.testing.assert_equal(system.part.by_ids(range(9, 4, -1)).id, [i for i in sorted(ids, key=lambda i:-i) if i > 4 and i <= 9]) - # Check that negative start and end on slices are not accepted - with self.assertRaises(IndexError): - system.part[-1:] - with self.assertRaises(IndexError): - system.part[:-1] # Setting particle properties on a slice - system.part[:5].pos = 0, 0, 0 - np.testing.assert_equal(system.part[:].pos, + system.part.by_ids(range(5)).pos = 0, 0, 0 + np.testing.assert_equal(system.part.all().pos, [pos[i] if ids[i] >= 5 else [0, 0, 0] for i in np.argsort(ids)]) # Slice access via explicit list of ids - np.testing.assert_equal(system.part[ids[1:4]].id, ids[1:4]) + np.testing.assert_equal(system.part.by_ids(ids[1:4]).id, ids[1:4]) # Check that ids passed in an explicit list must exist with self.assertRaises(IndexError): - system.part[99, 3] + system.part.by_ids([99, 3]) # Check that wrong types are not accepted with self.assertRaises(TypeError): - system.part[[ids[0], 1.2]] + system.part.by_ids([[ids[0], 1.2]]) def test_to_dict(self): self.system.part.clear() diff --git a/testsuite/python/particle_slice.py b/testsuite/python/particle_slice.py index 1c1fe85a870..26a40eb393d 100644 --- a/testsuite/python/particle_slice.py +++ b/testsuite/python/particle_slice.py @@ -22,80 +22,83 @@ class ParticleSliceTest(ut.TestCase): - state = [[0, 0, 0], [0, 0, 1]] - system = espressomd.System(box_l=[10, 10, 10]) + system = espressomd.System(box_l=3 * [10]) - def __init__(self, *args, **kwargs): - super().__init__(*args, **kwargs) - self.system.part.clear() - for i in range(4): - self.system.part.add(pos=[0, 0, i]) + def setUp(self): + self.fix_flags = [[False, False, False], [False, False, True]] + self.p0 = self.system.part.add(pos=[0, 0, 0]) + self.p1 = self.system.part.add(pos=[0, 0, 1]) + self.p2 = self.system.part.add(pos=[0, 0, 2]) + self.p3 = self.system.part.add(pos=[0, 0, 3]) + + self.all_partcls = self.system.part.all() + self.p0p1 = self.system.part.by_ids([0, 1]) + self.p2p3 = self.system.part.by_ids([2, 3]) if espressomd.has_features(["EXTERNAL_FORCES"]): - self.system.part[1].fix = self.state[1] - np.testing.assert_array_equal( - np.copy(self.system.part[0].fix), self.state[0]) + self.p1.fix = self.fix_flags[1] np.testing.assert_array_equal( - np.copy(self.system.part[1].fix), self.state[1]) + np.copy(self.p0.fix), self.fix_flags[0]) np.testing.assert_array_equal( - np.copy(self.system.part[:2].fix), self.state) - xs = self.system.part[:].pos - for i in range(len(xs)): + np.copy(self.p1.fix), self.fix_flags[1]) np.testing.assert_array_equal( - xs[i], np.copy(self.system.part[i].pos)) + np.copy(self.p0p1.fix), self.fix_flags) + + def tearDown(self): + self.system.part.clear() @utx.skipIfMissingFeatures(["EXTERNAL_FORCES"]) def test_1_set_different_values(self): - self.state[0] = [1, 0, 0] - self.state[1] = [1, 0, 0] - self.system.part[:2].fix = self.state + self.fix_flags[0] = [1, 0, 0] + self.fix_flags[1] = [1, 0, 0] + self.p0p1.fix = self.fix_flags np.testing.assert_array_equal( - np.copy(self.system.part[:2].fix), self.state) + np.copy(self.p0p1.fix), self.fix_flags) @utx.skipIfMissingFeatures(["EXTERNAL_FORCES"]) def test_2_set_same_value(self): - self.state[0] = [0, 1, 0] - self.state[1] = [0, 1, 0] - self.system.part[:2].fix = self.state[1] + self.fix_flags[0] = [0, 1, 0] + self.fix_flags[1] = [0, 1, 0] + self.p0p1.fix = self.fix_flags[1] np.testing.assert_array_equal( - np.copy(self.system.part[:2].fix), self.state) + np.copy(self.p0p1.fix), self.fix_flags) @utx.skipIfMissingFeatures(["EXTERNAL_FORCES"]) def test_3_set_one_value(self): - self.state[1] = [0, 0, 1] - self.system.part[1:2].fix = self.state[1] + self.fix_flags[1] = [0, 0, 1] + self.p1.fix = self.fix_flags[1] np.testing.assert_array_equal( - np.copy(self.system.part[:2].fix), self.state) + np.copy(self.p0p1.fix), self.fix_flags) @utx.skipIfMissingFeatures(["EXTERNAL_FORCES"]) def test_4_str(self): - self.assertEqual(repr(self.system.part[0].fix), - repr(np.array([0, 1, 0]))) - self.assertEqual(repr(self.system.part[:2].fix), - repr(np.array([[0, 1, 0], [0, 0, 1]]))) + self.assertEqual(repr(self.p0.fix), + repr(np.array([0, 0, 0]))) + self.assertEqual(repr(self.p0p1.fix), + repr(np.array([[0, 0, 0], [0, 0, 1]]))) def test_pos_str(self): - self.system.part[0].pos = [0, 0, 0] - self.system.part[1].pos = [0, 0, 1] - self.assertEqual(repr(self.system.part[0].pos), + self.p0.pos = [0, 0, 0] + self.p1.pos = [0, 0, 1] + self.assertEqual(repr(self.p0.pos), repr(np.array([0.0, 0.0, 0.0]))) - self.assertEqual(repr(self.system.part[:2].pos), + self.assertEqual(repr(self.p0p1.pos), repr(np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]))) @utx.skipIfMissingFeatures(["ELECTROSTATICS"]) def test_scalar(self): - self.system.part[:1].q = 1.3 - self.assertEqual(self.system.part[0].q, 1.3) - self.system.part[:2].q = 2.0 - self.assertEqual(self.system.part[0].q, 2) - self.assertEqual(self.system.part[1].q, 2) - self.system.part[:2].q = 3 - self.assertEqual(self.system.part[0].q, 3) - self.assertEqual(self.system.part[1].q, 3) - self.system.part[:2].q = [-1, 1.0] - self.assertEqual(self.system.part[0].q, -1) - self.assertEqual(self.system.part[1].q, 1) - qs = self.system.part[:2].q + self.p0.q = 1.3 + self.assertEqual(self.p0.q, 1.3) + self.p0p1.q = 2.0 + self.assertEqual(self.p0.q, 2) + self.assertEqual(self.p1.q, 2) + self.p0p1.q = 3 + self.assertEqual(self.p0.q, 3) + self.assertEqual(self.p1.q, 3) + self.p0p1.q = [-1, 1.0] + self.assertEqual(self.p0.q, -1) + self.assertEqual(self.p1.q, 1) + qs = self.p0p1.q self.assertEqual(qs[0], -1) self.assertEqual(qs[1], 1) @@ -108,246 +111,245 @@ def test_bonds(self): # tuple b = fene, 0 - self.system.part[2].bonds = b - self.assertEqual(self.system.part[:].bonds, [(), (), ((fene, 0),), ()]) + self.p2.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0),), ()]) # list - self.system.part[:].bonds = [] + self.all_partcls.bonds = [] b = [fene, 0] - self.system.part[2].bonds = b - self.assertEqual(self.system.part[:].bonds, [(), (), ((fene, 0),), ()]) + self.p2.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0),), ()]) # nested list single - self.system.part[:].bonds = [] + self.all_partcls.bonds = [] b = [[fene, 0]] - self.system.part[2].bonds = b - self.assertEqual(self.system.part[:].bonds, [(), (), ((fene, 0),), ()]) + self.p2.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0),), ()]) # nested list multi - self.system.part[:].bonds = [] + self.all_partcls.bonds = [] b = [[fene, 0], [fene, 1]] - self.system.part[2].bonds = b - self.assertEqual(self.system.part[:].bonds, + self.p2.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0), (fene, 1)), ()]) # nested tuple single - self.system.part[:].bonds = [] + self.all_partcls.bonds = [] b = ((fene, 0),) - self.system.part[2].bonds = b - self.assertEqual(self.system.part[:].bonds, [(), (), ((fene, 0),), ()]) + self.p2.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0),), ()]) # nested tuple multi - self.system.part[:].bonds = [] + self.all_partcls.bonds = [] b = ((fene, 0), (fene, 1)) - self.system.part[2].bonds = b - self.assertEqual(self.system.part[:].bonds, + self.p2.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0), (fene, 1)), ()]) # Add/Del bonds - self.system.part[:].bonds = [] - self.system.part[2].add_bond((fene, 0)) - self.assertEqual(self.system.part[:].bonds, [(), (), ((fene, 0),), ()]) - self.system.part[2].add_bond((fene, 1)) - self.assertEqual(self.system.part[:].bonds, + self.all_partcls.bonds = [] + self.p2.add_bond((fene, 0)) + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0),), ()]) + self.p2.add_bond((fene, 1)) + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0), (fene, 1)), ()]) - self.system.part[2].delete_bond((fene, 1)) - self.assertEqual(self.system.part[:].bonds, [(), (), ((fene, 0),), ()]) - self.system.part[2].delete_bond((fene, 0)) - self.assertEqual(self.system.part[:].bonds, [(), (), (), ()]) + self.p2.delete_bond((fene, 1)) + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0),), ()]) + self.p2.delete_bond((fene, 0)) + self.assertEqual(self.all_partcls.bonds, [(), (), (), ()]) - self.system.part[:].bonds = [] - self.system.part[2].add_bond((fene, 0)) - self.system.part[2].add_bond((fene, 1)) - self.system.part[2].delete_all_bonds() - self.assertEqual(self.system.part[:].bonds, [(), (), (), ()]) + self.all_partcls.bonds = [] + self.p2.add_bond((fene, 0)) + self.p2.add_bond((fene, 1)) + self.p2.delete_all_bonds() + self.assertEqual(self.all_partcls.bonds, [(), (), (), ()]) # Slices # tuple for all - self.system.part[:].bonds = [] + self.all_partcls.bonds = [] b = fene, 0 - self.system.part[2:].bonds = b - self.assertEqual(self.system.part[:].bonds, + self.p2p3.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0),), ((fene, 0),)]) # list for all - self.system.part[:].bonds = [] + self.all_partcls.bonds = [] b = [fene, 0] - self.system.part[2:].bonds = b - self.assertEqual(self.system.part[:].bonds, + self.p2p3.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0),), ((fene, 0),)]) # nested list single for all - self.system.part[:].bonds = [] + self.all_partcls.bonds = [] b = [[fene, 0]] - self.system.part[2:].bonds = b - self.assertEqual(self.system.part[:].bonds, + self.p2p3.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0),), ((fene, 0),)]) # nested list multi for all - self.system.part[:].bonds = [] + self.all_partcls.bonds = [] b = [[fene, 0], [fene, 1]] - self.system.part[2:].bonds = b - self.assertEqual(self.system.part[:].bonds, + self.p2p3.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0), (fene, 1)), ((fene, 0), (fene, 1))]) # tuples for each - self.system.part[:].bonds = [] + self.all_partcls.bonds = [] b = (((fene, 0),), ((fene, 1),)) - self.system.part[2:].bonds = b - self.assertEqual(self.system.part[:].bonds, + self.p2p3.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0),), ((fene, 1),)]) # lists for each - self.system.part[:].bonds = [] + self.all_partcls.bonds = [] b = [[[fene, 0]], [[fene, 1]]] - self.system.part[2:].bonds = b - self.assertEqual(self.system.part[:].bonds, + self.p2p3.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0),), ((fene, 1),)]) # multi tuples for each - self.system.part[:].bonds = [] + self.all_partcls.bonds = [] b = (((fene, 0), (fene, 1)), ((fene, 0), (fene, 1))) - self.system.part[2:].bonds = b - self.assertEqual(self.system.part[:].bonds, + self.p2p3.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0), (fene, 1)), ((fene, 0), (fene, 1))]) # multi lists for each - self.system.part[:].bonds = [] + self.all_partcls.bonds = [] b = [[[fene, 0], [fene, 1]], [[fene, 0], [fene, 1]]] - self.system.part[2:].bonds = b - self.assertEqual(self.system.part[:].bonds, + self.p2p3.bonds = b + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0), (fene, 1)), ((fene, 0), (fene, 1))]) # Add/Del bonds - self.system.part[:].bonds = [] - self.system.part[2:].add_bond((fene, 0)) - self.assertEqual(self.system.part[:].bonds, + self.all_partcls.bonds = [] + self.p2p3.add_bond((fene, 0)) + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0),), ((fene, 0),)]) - self.system.part[2:].add_bond((fene, 1)) - self.assertEqual(self.system.part[:].bonds, + self.p2p3.add_bond((fene, 1)) + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0), (fene, 1)), ((fene, 0), (fene, 1))]) - self.system.part[2:].delete_bond((fene, 1)) - self.assertEqual(self.system.part[:].bonds, + self.p2p3.delete_bond((fene, 1)) + self.assertEqual(self.all_partcls.bonds, [(), (), ((fene, 0),), ((fene, 0),)]) - self.system.part[2:].delete_bond((fene, 0)) - self.assertEqual(self.system.part[:].bonds, [(), (), (), ()]) + self.p2p3.delete_bond((fene, 0)) + self.assertEqual(self.all_partcls.bonds, [(), (), (), ()]) b = [[[fene, 0], [fene, 1]], [[fene, 0], [fene, 1]]] - self.system.part[2:].bonds = b - self.system.part[:].delete_all_bonds() - self.assertEqual(self.system.part[:].bonds, [(), (), (), ()]) + self.p2p3.bonds = b + self.all_partcls.delete_all_bonds() + self.assertEqual(self.all_partcls.bonds, [(), (), (), ()]) @utx.skipIfMissingFeatures(["EXCLUSIONS"]) def test_exclusions(self): def assert_exclusions_equal(B): - A = self.system.part[:].exclusions + A = self.all_partcls.exclusions self.assertEqual([x.tolist() for x in A], B) # Setter # int - self.system.part[:].exclusions = [] + self.all_partcls.exclusions = [] b = 1 - self.system.part[2].exclusions = b + self.p2.exclusions = b assert_exclusions_equal([[], [2], [1], []]) # single list - self.system.part[:].exclusions = [] + self.all_partcls.exclusions = [] b = [1] - self.system.part[2].exclusions = b + self.p2.exclusions = b assert_exclusions_equal([[], [2], [1], []]) # tuple - self.system.part[:].exclusions = [] + self.all_partcls.exclusions = [] b = (0, 1) - self.system.part[2].exclusions = b + self.p2.exclusions = b assert_exclusions_equal([[2], [2], [0, 1], []]) # list - self.system.part[:].exclusions = [] + self.all_partcls.exclusions = [] b = [0, 1] - self.system.part[2].exclusions = b + self.p2.exclusions = b assert_exclusions_equal([[2], [2], [0, 1], []]) # Add/Del exclusions - self.system.part[:].exclusions = [] - self.system.part[2].add_exclusion(1) + self.all_partcls.exclusions = [] + self.p2.add_exclusion(1) assert_exclusions_equal([[], [2], [1], []]) - self.system.part[2].add_exclusion(0) + self.p2.add_exclusion(0) assert_exclusions_equal([[2], [2], [1, 0], []]) - self.system.part[2].delete_exclusion(0) + self.p2.delete_exclusion(0) assert_exclusions_equal([[], [2], [1], []]) - self.system.part[2].delete_exclusion(1) + self.p2.delete_exclusion(1) assert_exclusions_equal([[], [], [], []]) # Slices # single list for all - self.system.part[:].exclusions = [] + self.all_partcls.exclusions = [] b = [1] - self.system.part[2:].exclusions = b + self.p2p3.exclusions = b assert_exclusions_equal([[], [2, 3], [1], [1]]) # list for all - self.system.part[:].exclusions = [] + self.all_partcls.exclusions = [] b = [0, 1] - self.system.part[2:].exclusions = b + self.p2p3.exclusions = b assert_exclusions_equal([[2, 3], [2, 3], [0, 1], [0, 1]]) # single list for each - self.system.part[:].exclusions = [] + self.all_partcls.exclusions = [] b = [[0], [0]] - self.system.part[2:].exclusions = b + self.p2p3.exclusions = b assert_exclusions_equal([[2, 3], [], [0], [0]]) # multi list for each - self.system.part[:].exclusions = [] + self.all_partcls.exclusions = [] b = [[0, 1], [0, 1]] - self.system.part[2:].exclusions = b + self.p2p3.exclusions = b assert_exclusions_equal([[2, 3], [2, 3], [0, 1], [0, 1]]) # Add/Del exclusions - self.system.part[:].exclusions = [] - self.system.part[2:].add_exclusion(1) + self.all_partcls.exclusions = [] + self.p2p3.add_exclusion(1) assert_exclusions_equal([[], [2, 3], [1], [1]]) - self.system.part[2:].add_exclusion(0) + self.p2p3.add_exclusion(0) assert_exclusions_equal([[2, 3], [2, 3], [1, 0], [1, 0]]) - self.system.part[2:].delete_exclusion(0) + self.p2p3.delete_exclusion(0) assert_exclusions_equal([[], [2, 3], [1], [1]]) - self.system.part[2:].delete_exclusion(1) + self.p2p3.delete_exclusion(1) assert_exclusions_equal([[], [], [], []]) @utx.skipIfMissingFeatures("VIRTUAL_SITES_RELATIVE") def test_vs_relative(self): self.system.part.clear() - self.system.part.add(pos=[0, 0, 0]) - self.system.part.add(pos=[0, 0, 0]) - self.system.part.add(pos=[0, 0, 0]) - self.system.part.add(pos=[0, 0, 0]) - self.system.part[0].vs_relative = [1, 1.0, (1.0, 1.0, 1.0, 1.0)] + p0 = self.system.part.add(pos=[0, 0, 0]) + self.system.part.add(pos=3 * [[0, 0, 0]]) + p0.vs_relative = [1, 1.0, (1.0, 1.0, 1.0, 1.0)] + all_partcls = self.system.part.all() - self.assertEqual(repr(self.system.part[:].vs_relative), + self.assertEqual(repr(all_partcls.vs_relative), repr([(1, 1.0, np.array([1., 1., 1., 1.])), (0, 0.0, np.array([1., 0., 0., 0.])), (0, 0.0, np.array([1., 0., 0., 0.])), (0, 0.0, np.array([1., 0., 0., 0.]))])) - self.system.part[:].vs_relative = [1, 1.0, (1.0, 1.0, 1.0, 1.0)] + all_partcls.vs_relative = [1, 1.0, (1.0, 1.0, 1.0, 1.0)] - self.assertEqual(repr(self.system.part[:].vs_relative), + self.assertEqual(repr(all_partcls.vs_relative), repr([(1, 1.0, np.array([1., 1., 1., 1.])), (1, 1.0, np.array([1., 1., 1., 1.])), (1, 1.0, np.array([1., 1., 1., 1.])), (1, 1.0, np.array([1., 1., 1., 1.]))])) - self.system.part[:].vs_relative = [[1, 1.0, (1.0, 1.0, 1.0, 1.0)], - [1, 2.0, (1.0, 1.0, 1.0, 1.0)], - [1, 3.0, (1.0, 1.0, 1.0, 1.0)], - [1, 4.0, (1.0, 1.0, 1.0, 1.0)]] + all_partcls.vs_relative = [[1, 1.0, (1.0, 1.0, 1.0, 1.0)], + [1, 2.0, (1.0, 1.0, 1.0, 1.0)], + [1, 3.0, (1.0, 1.0, 1.0, 1.0)], + [1, 4.0, (1.0, 1.0, 1.0, 1.0)]] - self.assertEqual(repr(self.system.part[:].vs_relative), + self.assertEqual(repr(all_partcls.vs_relative), repr([(1, 1.0, np.array([1., 1., 1., 1.])), (1, 2.0, np.array([1., 1., 1., 1.])), (1, 3.0, np.array([1., 1., 1., 1.])), @@ -356,11 +358,12 @@ def test_vs_relative(self): def test_multiadd(self): self.system.part.clear() positions = [[1, 2, 3], [4, 5, 6], [7, 8, 9]] - self.system.part.add(pos=positions) + partcls_multiadd = self.system.part.add(pos=positions) for p in positions: self.system.part.add(pos=p) + partcls_singleadd = self.system.part.by_ids(range(3, 6)) np.testing.assert_allclose( - np.copy(self.system.part[:3].pos), np.copy(self.system.part[3:6].pos)) + np.copy(partcls_multiadd.pos), np.copy(partcls_singleadd.pos)) with self.assertRaises(ValueError): self.system.part.add(pos=([1, 1, 1], [2, 2, 2]), type=0) @@ -368,21 +371,23 @@ def test_multiadd(self): self.system.part.add(pos=([1, 1, 1], [2, 2, 2]), type=(0, 1, 2)) self.system.part.clear() - self.system.part.add(pos=([1, 1, 1], [2, 2, 2]), type=(0, 1)) - self.assertEqual(self.system.part[0].type, 0) - self.assertEqual(self.system.part[1].type, 1) + p0, p1 = self.system.part.add(pos=([1, 1, 1], [2, 2, 2]), type=(0, 1)) + self.assertEqual(p0.type, 0) + self.assertEqual(p1.type, 1) def test_empty(self): - np.testing.assert_array_equal(self.system.part[0:0].pos, np.empty(0)) + np.testing.assert_array_equal( + self.system.part.by_ids( + []).pos, np.empty(0)) def test_len(self): - self.assertEqual(len(self.system.part[0:0]), 0) - self.assertEqual(len(self.system.part[0:1]), 1) - self.assertEqual(len(self.system.part[0:2]), 2) + self.assertEqual(len(self.system.part.by_ids([])), 0) + self.assertEqual(len(self.system.part.by_ids([0])), 1) + self.assertEqual(len(self.system.part.by_ids([0, 1, 2])), 3) def test_non_existing_property(self): with self.assertRaises(AttributeError): - self.system.part[:].thispropertydoesnotexist = 1.0 + self.all_partcls.thispropertydoesnotexist = 1.0 if __name__ == "__main__": diff --git a/testsuite/python/polymer_diamond.py b/testsuite/python/polymer_diamond.py index 4a622ba7a34..4692eb6a456 100644 --- a/testsuite/python/polymer_diamond.py +++ b/testsuite/python/polymer_diamond.py @@ -85,7 +85,7 @@ def test_particle_properties(self): np.testing.assert_allclose(self.charged_mono.q, len(self.charged_mono) * [self.diamond_params['val_cM']]) # particle id - self.assertGreaterEqual(min(self.system.part[:].id), + self.assertGreaterEqual(min(self.system.part.all().id), self.diamond_params['start_id']) def test_bonds(self): @@ -109,9 +109,10 @@ def test_connected(self): test that all particles in the polymer are connected by pushing one particle """ - self.system.part[self.diamond_params['start_id']].ext_force = 3 * [2] + self.system.part.by_id( + self.diamond_params['start_id']).ext_force = 3 * [2] self.system.integrator.run(200) - vels = np.linalg.norm(self.system.part[:].v, axis=1) + vels = np.linalg.norm(self.system.part.all().v, axis=1) self.assertGreater(np.min(vels), 1e-6) def test_geometry(self): diff --git a/testsuite/python/pressure.py b/testsuite/python/pressure.py index 1edb2a176f5..5a8bc3c2e35 100644 --- a/testsuite/python/pressure.py +++ b/testsuite/python/pressure.py @@ -126,8 +126,9 @@ def test(self): p2.v = [27.0, 23.0, 17.0] p3.v = [13.0, 11.0, 19.0] - pos = system.part[:].pos - vel = system.part[:].v + partcls = system.part.all() + pos = partcls.pos + vel = partcls.v sim_pressure_tensor = system.analysis.pressure_tensor() sim_pressure_tensor_kinetic = np.copy(sim_pressure_tensor['kinetic']) diff --git a/testsuite/python/random_pairs.py b/testsuite/python/random_pairs.py index 2cbb5a48067..e08b01956b9 100644 --- a/testsuite/python/random_pairs.py +++ b/testsuite/python/random_pairs.py @@ -66,7 +66,8 @@ def pairs_n2(self, dist): pairs = [] parts = self.system.part for p in self.all_pairs: - if self.system.distance(parts[p[0]], parts[p[1]]) <= dist: + if self.system.distance(parts.by_id( + p[0]), parts.by_id(p[1])) <= dist: pairs.append(p) return set(pairs) diff --git a/testsuite/python/rdf.py b/testsuite/python/rdf.py index e5a049d29c4..29771c05dc9 100644 --- a/testsuite/python/rdf.py +++ b/testsuite/python/rdf.py @@ -51,7 +51,7 @@ def test_single_type(self): r_bins = 50 r_min = 0.5 * dx r_max = r_bins * dx - obs = espressomd.observables.RDF(ids1=system.part[:].id, min_r=r_min, + obs = espressomd.observables.RDF(ids1=system.part.all().id, min_r=r_min, max_r=r_max, n_r_bins=r_bins) rdf = obs.calculate() r = obs.bin_centers() @@ -72,12 +72,13 @@ def test_mixed(self): for i in range(n_part): system.part.add( pos=[i * dx, 0.5 * system.box_l[1], 0.5 * system.box_l[2]], type=(i % 2)) + partcls = system.part.all() r_bins = 50 r_min = 0.5 * dx r_max = r_bins * dx - obs = espressomd.observables.RDF(ids1=system.part[:].id[0::2], - ids2=system.part[:].id[1::2], + obs = espressomd.observables.RDF(ids1=partcls.id[0::2], + ids2=partcls.id[1::2], min_r=r_min, max_r=r_max, n_r_bins=r_bins) rdf01 = obs.calculate() @@ -93,8 +94,8 @@ def test_mixed(self): np.testing.assert_allclose(parts_in_bin[1::2], 0.0) # Check symmetry - obs = espressomd.observables.RDF(ids1=system.part[:].id[1::2], - ids2=system.part[:].id[0::2], + obs = espressomd.observables.RDF(ids1=partcls.id[1::2], + ids2=partcls.id[0::2], min_r=r_min, max_r=r_max, n_r_bins=r_bins) rdf10 = obs.calculate() @@ -104,9 +105,9 @@ def test_mixed(self): def test_rdf_interface(self): # test setters and getters system = self.system - system.part.add(pos=4 * [(0, 0, 0)], type=[0, 1, 0, 1]) - pids1 = system.part[:].id[0::2] - pids2 = system.part[:].id[1::2] + partcls = system.part.add(pos=4 * [(0, 0, 0)], type=[0, 1, 0, 1]) + pids1 = partcls.id[0::2] + pids2 = partcls.id[1::2] params = { 'ids1': pids1, 'ids2': pids2, @@ -117,8 +118,8 @@ def test_rdf_interface(self): # check pids np.testing.assert_array_equal(np.copy(observable.ids1), pids1) np.testing.assert_array_equal(np.copy(observable.ids2), pids2) - new_pids1 = [system.part[:].id[0]] - new_pids2 = [system.part[:].id[1]] + new_pids1 = [partcls.id[0]] + new_pids2 = [partcls.id[1]] observable = espressomd.observables.RDF( **{**params, 'ids1': new_pids1, 'ids2': new_pids2}) np.testing.assert_array_equal(np.copy(observable.ids1), new_pids1) diff --git a/testsuite/python/rescale.py b/testsuite/python/rescale.py index 6c89379ac29..eb314beb98b 100644 --- a/testsuite/python/rescale.py +++ b/testsuite/python/rescale.py @@ -33,7 +33,11 @@ class RescaleTest(ut.TestCase): def setUp(self): N = 100 self.system.box_l = 3 * [10] - self.system.part.add(pos=self.system.box_l * np.random.random((N, 3))) + self.partcls = self.system.part.add( + pos=self.system.box_l * np.random.random((N, 3))) + + def tearDown(self): + self.system.part.clear() def test_iso(self): """Test 'isotropic' case (dir="xyz"). @@ -41,9 +45,9 @@ def test_iso(self): scale = 1.3 new_box_l = scale * self.system.box_l[0] - old_pos = self.system.part[:].pos + old_pos = self.partcls.pos self.system.change_volume_and_rescale_particles(new_box_l) - new_pos = self.system.part[:].pos + new_pos = self.partcls.pos max_diff = np.max(np.abs(new_pos / old_pos - scale)) self.assertAlmostEqual(0., max_diff, places=10) @@ -54,9 +58,9 @@ def dir_test(self, dir): scale = 0.7 new_box_l = scale * self.system.box_l[dir] - old_pos = self.system.part[:].pos + old_pos = self.partcls.pos self.system.change_volume_and_rescale_particles(new_box_l, dir=dir) - new_pos = self.system.part[:].pos + new_pos = self.partcls.pos for i in range(3): if i == dir: diff --git a/testsuite/python/rotational-diffusion-aniso.py b/testsuite/python/rotational-diffusion-aniso.py index 3afb473c75f..5959a900efc 100644 --- a/testsuite/python/rotational-diffusion-aniso.py +++ b/testsuite/python/rotational-diffusion-aniso.py @@ -56,10 +56,10 @@ def tearDown(self): def add_particles(self): positions = np.random.random((self.n_part, 3)) * self.system.box_l[0] - self.system.part.add(pos=positions, rotation=self.n_part * [(1, 1, 1)], - rinertia=self.n_part * [self.J]) + self.partcls = self.system.part.add(pos=positions, rotation=self.n_part * [(1, 1, 1)], + rinertia=self.n_part * [self.J]) if espressomd.has_features("ROTATION"): - self.system.part[:].omega_body = [0.0, 0.0, 0.0] + self.partcls.omega_body = [0.0, 0.0, 0.0] def set_anisotropic_param(self): """ @@ -142,7 +142,7 @@ def check_rot_diffusion(self): # The body angular velocity is rotated now, but there is only the # thermal velocity, hence, this does not impact the test and its # physical context. - self.system.part[:].quat = [1.0, 0.0, 0.0, 0.0] + self.partcls.quat = [1.0, 0.0, 0.0, 0.0] # Average direction cosines # Diagonal ones: dcosjj_validate = np.zeros((3)) diff --git a/testsuite/python/scafacos_dipoles_1d_2d.py b/testsuite/python/scafacos_dipoles_1d_2d.py index 78723e9fb02..9d1cde69274 100644 --- a/testsuite/python/scafacos_dipoles_1d_2d.py +++ b/testsuite/python/scafacos_dipoles_1d_2d.py @@ -76,7 +76,7 @@ def test_scafacos(self): data = np.genfromtxt(tests_common.abspath( file_prefix + "_reference_data_forces_torques.dat")) s.part.add(pos=data[:, 1:4], dip=data[:, 4:7]) - s.part[:].rotation = (1, 1, 1) + s.part.all().rotation = 3 * [True] if dim == 2: scafacos = espressomd.magnetostatics.Scafacos( @@ -122,8 +122,8 @@ def test_scafacos(self): # Calculate errors - err_f = self.vector_error(s.part[:].f, data[:, 7:10]) - err_t = self.vector_error(s.part[:].torque_lab, data[:, 10:13]) + err_f = self.vector_error(s.part.all().f, data[:, 7:10]) + err_t = self.vector_error(s.part.all().torque_lab, data[:, 10:13]) err_e = s.analysis.energy()["dipolar"] - ref_E tol_f = 2E-3 diff --git a/testsuite/python/scafacos_interface.py b/testsuite/python/scafacos_interface.py index 82aa9e0c7f7..59069e141e5 100644 --- a/testsuite/python/scafacos_interface.py +++ b/testsuite/python/scafacos_interface.py @@ -32,7 +32,7 @@ class ScafacosInterface(ut.TestCase): system = espressomd.System(box_l=3 * [5]) system.time_step = 0.01 system.cell_system.skin = 0.5 - system.periodicity = [1, 1, 1] + system.periodicity = 3 * [True] def tearDown(self): self.system.part.clear() @@ -153,8 +153,8 @@ def p3m_data(self): system.integrator.run(0, recalc_forces=True) ref_E_coulomb = system.analysis.energy()["coulomb"] ref_E_dipoles = system.analysis.energy()["dipolar"] - ref_forces = np.copy(system.part[:].f) - ref_torques = np.copy(system.part[:].torque_lab) + ref_forces = np.copy(system.part.all().f) + ref_torques = np.copy(system.part.all().torque_lab) system.actors.clear() @@ -197,8 +197,8 @@ def fcs_data(self): system.integrator.run(0, recalc_forces=True) ref_E_coulomb = system.analysis.energy()["coulomb"] ref_E_dipoles = system.analysis.energy()["dipolar"] - ref_forces = np.copy(system.part[:].f) - ref_torques = np.copy(system.part[:].torque_lab) + ref_forces = np.copy(system.part.all().f) + ref_torques = np.copy(system.part.all().torque_lab) system.actors.clear() diff --git a/testsuite/python/tabulated.py b/testsuite/python/tabulated.py index 921ca306b7b..ce157e9dec3 100644 --- a/testsuite/python/tabulated.py +++ b/testsuite/python/tabulated.py @@ -46,9 +46,9 @@ def tearDown(self): self.system.part.clear() def check(self): - p0, p1 = self.system.part[:] + p0, p1 = self.system.part.all() # Below cutoff - np.testing.assert_allclose(np.copy(self.system.part[:].f), 0.0) + np.testing.assert_allclose(np.copy(self.system.part.all().f), 0.0) for z in np.linspace(0, self.max_ - self.min_, 200, endpoint=False): p1.pos = [5., 5., 6. + z] @@ -86,7 +86,7 @@ def test_bonded(self): self.assertAlmostEqual(tb.params['min'], self.min_) self.assertAlmostEqual(tb.params['max'], self.max_) - p0, p1 = self.system.part[:] + p0, p1 = self.system.part.all() p0.add_bond((tb, p1)) self.check() diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index 83a353498f8..914e26052fc 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -127,7 +127,7 @@ def test_system_variables(self): np.copy(system.periodicity), self.ref_periodicity) def test_part(self): - p1, p2 = system.part[0:2] + p1, p2 = system.part.by_ids([0, 1]) np.testing.assert_allclose(np.copy(p1.pos), np.array([1.0, 2.0, 3.0])) np.testing.assert_allclose(np.copy(p2.pos), np.array([1.0, 1.0, 2.0])) np.testing.assert_allclose(np.copy(p1.f), particle_force0) @@ -139,7 +139,7 @@ def test_bonded_interactions_serialization(self): experience the force from a harmonic bond. The thermostat friction is negligible compared to the harmonic bond strength. ''' - p3, p4 = system.part[3:5] + p3, p4 = system.part.by_ids([3, 4]) np.testing.assert_allclose(np.copy(p3.pos), system.box_l / 2. - 1.) np.testing.assert_allclose(np.copy(p4.pos), system.box_l / 2. + 1.) np.testing.assert_allclose(np.copy(p3.f), -np.copy(p4.f), rtol=1e-4) @@ -153,7 +153,7 @@ def test_shape_based_constraints_serialization(self): experience the force from a shape-based constraint. The thermostat friction is negligible compared to the LJ force. ''' - p3, p4 = system.part[3:5] + p3, p4 = system.part.by_ids([3, 4]) old_force = np.copy(p3.f) system.constraints.remove(system.constraints[0]) system.integrator.run(0, recalc_forces=True) @@ -312,17 +312,18 @@ def test_bonded_inter(self): bond_ids = system.bonded_inter.call_method('get_bond_ids') self.assertEqual(len(bond_ids), len(system.bonded_inter)) # check bonded interactions - state = system.part[1].bonds[0][0].params + partcl_1 = system.part.by_id(1) + state = partcl_1.bonds[0][0].params reference = {'r_0': 0.0, 'k': 1.0, 'r_cut': 0.0} self.assertEqual(state, reference) - state = system.part[1].bonds[0][0].params + state = partcl_1.bonds[0][0].params self.assertEqual(state, reference) if 'THERM.LB' not in modes: - state = system.part[1].bonds[1][0].params + state = partcl_1.bonds[1][0].params reference = {'temp_com': 0., 'gamma_com': 0., 'temp_distance': 0.2, 'gamma_distance': 0.5, 'r_cut': 2.0, 'seed': 51} self.assertEqual(state, reference) - state = system.part[1].bonds[1][0].params + state = partcl_1.bonds[1][0].params self.assertEqual(state, reference) # immersed boundary bonds self.assertEqual( @@ -359,7 +360,7 @@ def test_drude_helpers(self): @utx.skipIfMissingFeatures(['VIRTUAL_SITES', 'VIRTUAL_SITES_RELATIVE']) def test_virtual_sites(self): - self.assertTrue(system.part[1].virtual) + self.assertTrue(system.part.by_id(1).virtual) self.assertIsInstance( system.virtual_sites, espressomd.virtual_sites.VirtualSitesRelative) @@ -492,9 +493,9 @@ def test_collision_detection(self): @utx.skipIfMissingFeatures('EXCLUSIONS') def test_exclusions(self): - self.assertEqual(list(system.part[0].exclusions), [2]) - self.assertEqual(list(system.part[1].exclusions), [2]) - self.assertEqual(list(system.part[2].exclusions), [0, 1]) + self.assertEqual(list(system.part.by_id(0).exclusions), [2]) + self.assertEqual(list(system.part.by_id(1).exclusions), [2]) + self.assertEqual(list(system.part.by_id(2).exclusions), [0, 1]) @ut.skipIf(not LB or not (espressomd.has_features("LB_BOUNDARIES") or espressomd.has_features("LB_BOUNDARIES_GPU")), "Missing features") diff --git a/testsuite/python/tests_common.py b/testsuite/python/tests_common.py index eeb21d89fd8..6285c8b5221 100644 --- a/testsuite/python/tests_common.py +++ b/testsuite/python/tests_common.py @@ -89,7 +89,7 @@ def verify_lj_forces(system, tolerance, ids_to_skip=()): # Initialize dict with expected forces f_expected = {} - for pid in system.part[:].id: + for pid in system.part.all().id: f_expected[pid] = np.zeros(3) # Cache some stuff to speed up pair loop @@ -98,7 +98,7 @@ def verify_lj_forces(system, tolerance, ids_to_skip=()): non_bonded_inter = system.non_bonded_inter # LJ parameters lj_params = {} - all_types = np.unique(system.part[:].type) + all_types = np.unique(system.part.all().type) for i in all_types: for j in all_types: lj_params[i, j] = non_bonded_inter[i, j].lennard_jones.get_params() diff --git a/testsuite/python/thermalized_bond.py b/testsuite/python/thermalized_bond.py index f5d6c449db4..1af3ad7df92 100644 --- a/testsuite/python/thermalized_bond.py +++ b/testsuite/python/thermalized_bond.py @@ -45,7 +45,7 @@ def test_com_langevin(self): """Test for COM thermalization.""" N = 200 - N2 = int(N / 2) + N_half = int(N / 2) self.system.part.clear() self.system.time_step = 0.02 self.system.periodicity = [0, 0, 0] @@ -53,9 +53,13 @@ def test_com_langevin(self): m1 = 1.0 m2 = 10.0 # Place particles - for i in range(0, N, 2): - self.system.part.add(pos=np.random.random(3), mass=m1) - self.system.part.add(pos=np.random.random(3), mass=m2) + + partcls_m1 = self.system.part.add( + pos=np.random.random( + (N_half, 3)), mass=N_half * [m1]) + partcls_m2 = self.system.part.add( + pos=np.random.random( + (N_half, 3)), mass=N_half * [m2]) t_dist = 0 g_dist = 0 @@ -67,7 +71,7 @@ def test_com_langevin(self): gamma_distance=g_dist, r_cut=2.0, seed=55) self.system.bonded_inter.add(thermalized_dist_bond) - for p1, p2 in zip(self.system.part[::2], self.system.part[1::2]): + for p1, p2 in zip(partcls_m1, partcls_m2): p1.add_bond((thermalized_dist_bond, p2)) # Warmup @@ -75,13 +79,13 @@ def test_com_langevin(self): # Sampling loops = 150 - v_stored = np.zeros((N2 * loops, 3)) + v_stored = np.zeros((N_half * loops, 3)) for i in range(loops): self.system.integrator.run(5) v_com = 1.0 / \ (m1 + m2) * \ - (m1 * self.system.part[::2].v + m2 * self.system.part[1::2].v) - v_stored[i * N2:(i + 1) * N2, :] = v_com + (m1 * partcls_m1.v + m2 * partcls_m2.v) + v_stored[i * N_half:(i + 1) * N_half, :] = v_com v_minmax = 5 bins = 50 @@ -93,7 +97,7 @@ def test_dist_langevin(self): """Test for dist thermalization.""" N = 100 - N2 = int(N / 2) + N_half = int(N / 2) self.system.part.clear() self.system.time_step = 0.02 self.system.periodicity = [1, 1, 1] @@ -101,9 +105,12 @@ def test_dist_langevin(self): m1 = 1.0 m2 = 10.0 # Place particles - for i in range(0, N, 2): - self.system.part.add(pos=np.random.random(3), mass=m1) - self.system.part.add(pos=np.random.random(3), mass=m2) + partcls_m1 = self.system.part.add( + pos=np.random.random( + (N_half, 3)), mass=N_half * [m1]) + partcls_m2 = self.system.part.add( + pos=np.random.random( + (N_half, 3)), mass=N_half * [m2]) t_dist = 2.0 g_dist = 4.0 @@ -115,7 +122,7 @@ def test_dist_langevin(self): gamma_distance=g_dist, r_cut=9, seed=51) self.system.bonded_inter.add(thermalized_dist_bond) - for p1, p2 in zip(self.system.part[::2], self.system.part[1::2]): + for p1, p2 in zip(partcls_m1, partcls_m2): p1.add_bond((thermalized_dist_bond, p2)) # Warmup @@ -123,11 +130,11 @@ def test_dist_langevin(self): # Sampling loops = 150 - v_stored = np.zeros((N2 * loops, 3)) + v_stored = np.zeros((N_half * loops, 3)) for i in range(loops): self.system.integrator.run(5) - v_dist = self.system.part[1::2].v - self.system.part[::2].v - v_stored[i * N2:(i + 1) * N2, :] = v_dist + v_dist = partcls_m2.v - partcls_m1.v + v_stored[i * N_half:(i + 1) * N_half, :] = v_dist v_minmax = 5 bins = 50 diff --git a/testsuite/python/thermostats_common.py b/testsuite/python/thermostats_common.py index ae876f4f827..205a1724ab3 100644 --- a/testsuite/python/thermostats_common.py +++ b/testsuite/python/thermostats_common.py @@ -79,23 +79,23 @@ def check_global(self, N, kT, steps, v_minmax, system = self.system # Place particles - system.part.add(pos=np.random.random((N, 3))) + partcls = system.part.add(pos=np.random.random((N, 3))) # Enable rotation if compiled in if espressomd.has_features("ROTATION"): - system.part[:].rotation = [1, 1, 1] + partcls.rotation = 3 * [True] # Warmup system.integrator.run(20) vel_obs = espressomd.observables.ParticleVelocities( - ids=system.part[:].id) + ids=partcls.id) vel_acc = espressomd.accumulators.TimeSeries(obs=vel_obs) system.auto_update_accumulators.add(vel_acc) if espressomd.has_features("ROTATION"): omega_obs = espressomd.observables.ParticleBodyAngularVelocities( - ids=system.part[:].id) + ids=partcls.id) omega_acc = espressomd.accumulators.TimeSeries(obs=omega_obs) system.auto_update_accumulators.add(omega_acc) @@ -104,9 +104,9 @@ def check_global(self, N, kT, steps, v_minmax, omega_stored = np.zeros((steps, N, 3)) for i in range(steps): system.integrator.run(1, recalc_forces=recalc_forces) - v_stored[i] = system.part[:].v + v_stored[i] = partcls.v if espressomd.has_features("ROTATION"): - omega_stored[i] = system.part[:].omega_body + omega_stored[i] = partcls.omega_body vel = vel_acc.time_series().reshape((-1, 3)) self.check_velocity_distribution(vel, v_minmax, n_bins, error_tol, kT) @@ -140,26 +140,26 @@ def check_per_particle(self, N, kT, gamma_local, steps, v_minmax, n_bins, Error tolerance. """ system = self.system - system.part.add(pos=np.random.random((N, 3))) + partcls = system.part.add(pos=np.random.random((N, 3))) if espressomd.has_features("ROTATION"): - system.part[:].rotation = [1, 1, 1] + partcls.rotation = [1, 1, 1] if espressomd.has_features("PARTICLE_ANISOTROPY"): gamma_local = 3 * [gamma_local] # Set different gamma on 2nd half of particles - system.part[N // 2:].gamma = gamma_local + system.part.by_ids(range(N // 2, N)).gamma = gamma_local system.integrator.run(50) vel_obs = espressomd.observables.ParticleVelocities( - ids=system.part[:].id) + ids=partcls.id) vel_acc = espressomd.accumulators.TimeSeries(obs=vel_obs) system.auto_update_accumulators.add(vel_acc) if espressomd.has_features("ROTATION"): omega_obs = espressomd.observables.ParticleBodyAngularVelocities( - ids=system.part[:].id) + ids=partcls.id) omega_acc = espressomd.accumulators.TimeSeries(obs=omega_obs) system.auto_update_accumulators.add(omega_acc) @@ -186,14 +186,15 @@ def check_noise_correlation(self, kT, steps, delta): """ system = self.system + partcls = system.part.all() vel_obs = espressomd.observables.ParticleVelocities( - ids=system.part[:].id) + ids=partcls.id) vel_series = espressomd.accumulators.TimeSeries(obs=vel_obs) system.auto_update_accumulators.add(vel_series) if espressomd.has_features("ROTATION"): - system.part[:].rotation = (1, 1, 1) + partcls.rotation = 3 * [True] omega_obs = espressomd.observables.ParticleBodyAngularVelocities( - ids=system.part[:].id) + ids=partcls.id) omega_series = espressomd.accumulators.TimeSeries(obs=omega_obs) system.auto_update_accumulators.add(omega_series) diff --git a/testsuite/python/thole.py b/testsuite/python/thole.py index 663126bb37d..8abb54be275 100644 --- a/testsuite/python/thole.py +++ b/testsuite/python/thole.py @@ -43,8 +43,14 @@ class TestThole(ut.TestCase): def setUp(self): self.system.time_step = 0.01 self.system.cell_system.skin = 0.4 - self.system.part.add(pos=[0, 0, 0], type=0, fix=[1, 1, 1], q=self.q1) - self.system.part.add(pos=[2, 0, 0], type=0, fix=[1, 1, 1], q=self.q2) + self.p0 = self.system.part.add( + pos=[ + 0, 0, 0], type=0, fix=[ + 1, 1, 1], q=self.q1) + self.p1 = self.system.part.add( + pos=[ + 2, 0, 0], type=0, fix=[ + 1, 1, 1], q=self.q2) p3m = espressomd.electrostatics.P3M( prefactor=COULOMB_PREFACTOR, accuracy=1e-6, mesh=3 * [52], cao=4) @@ -59,7 +65,7 @@ def test(self): ns = 100 for i in range(1, ns): x = 20.0 * i / ns - self.system.part[1].pos = [x, 0, 0] + self.p1.pos = [x, 0, 0] self.system.integrator.run(0) sd = x * self.thole_s @@ -74,7 +80,7 @@ def test(self): E = self.system.analysis.energy() E_tot = E["total"] - res_dForce.append(self.system.part[1].f[0] - F_calc) + res_dForce.append(self.p1.f[0] - F_calc) res_dEnergy.append(E_tot - E_calc) for f in res_dForce: diff --git a/testsuite/python/virtual_sites_relative.py b/testsuite/python/virtual_sites_relative.py index 3be08927621..9d0f2fc87a8 100644 --- a/testsuite/python/virtual_sites_relative.py +++ b/testsuite/python/virtual_sites_relative.py @@ -54,7 +54,7 @@ def verify_vs(self, vs, verify_velocity=True): vs_r = vs.vs_relative # Get related particle - rel = self.system.part[vs_r[0]] + rel = self.system.part.by_id(vs_r[0]) # Distance d = self.system.distance(rel, vs) @@ -190,8 +190,7 @@ def test_pos_vel_forces(self): f2 = np.array((3, 4, 5)) f3 = np.array((-4, 5, 6)) # Add forces to vs - p2 = system.part[2] - p3 = system.part[3] + p2, p3 = system.part.by_ids([2, 3]) p2.ext_force = f2 p3.ext_force = f3 system.integrator.run(0) @@ -240,22 +239,24 @@ def run_test_lj(self): # For n spheres, n/2 dumbells. for i in range(n // 2): # Type=1, i.e., no lj ia for the center of mass particles - system.part.add( - rotation=(1, 1, 1), id=3 * i, pos=np.random.random(3) * l, type=1, + p3i = system.part.add( + rotation=3 * [True], id=3 * i, pos=np.random.random(3) * l, type=1, omega_lab=0.3 * np.random.random(3), v=np.random.random(3)) # lj spheres - system.part.add(rotation=(1, 1, 1), id=3 * i + 1, - pos=system.part[3 * i].pos + - system.part[3 * i].director / 2., - type=0) - system.part.add(rotation=(1, 1, 1), id=3 * i + 2, - pos=system.part[3 * i].pos - - system.part[3 * i].director / 2., - type=0) - system.part[3 * i + 1].vs_auto_relate_to(3 * i) - self.verify_vs(system.part[3 * i + 1], verify_velocity=False) - system.part[3 * i + 2].vs_auto_relate_to(3 * i) - self.verify_vs(system.part[3 * i + 2], verify_velocity=False) + p3ip1 = system.part.add(rotation=3 * [True], + id=3 * i + 1, + pos=p3i.pos + + p3i.director / 2., + type=0) + p3ip2 = system.part.add(rotation=3 * [True], + id=3 * i + 2, + pos=p3i.pos - + p3i.director / 2., + type=0) + p3ip1.vs_auto_relate_to(p3i.id) + self.verify_vs(p3ip1, verify_velocity=False) + p3ip2.vs_auto_relate_to(p3i.id) + self.verify_vs(p3ip2, verify_velocity=False) system.integrator.run(0, recalc_forces=True) # interactions system.non_bonded_inter[0, 0].lennard_jones.set_params( @@ -277,8 +278,8 @@ def run_test_lj(self): system.integrator.run(2) # Check the virtual sites config,pos and vel of the lj spheres for j in range(int(n / 2)): - self.verify_vs(system.part[3 * j + 1]) - self.verify_vs(system.part[3 * j + 2]) + self.verify_vs(system.part.by_id(3 * j + 1)) + self.verify_vs(system.part.by_id(3 * j + 2)) # Verify lj forces on the particles. The non-virtual particles are # skipped because the forces on them originate from the vss and not @@ -288,7 +289,8 @@ def run_test_lj(self): # Test applying changes enegry_pre_change = system.analysis.energy()['total'] pressure_pre_change = system.analysis.pressure()['total'] - system.part[0].pos = system.part[0].pos + (2.2, -1.4, 4.2) + p0 = system.part.by_id(0) + p0.pos = p0.pos + (2.2, -1.4, 4.2) enegry_post_change = system.analysis.energy()['total'] pressure_post_change = system.analysis.pressure()['total'] self.assertNotAlmostEqual(enegry_pre_change, enegry_post_change) diff --git a/testsuite/python/writevtf.py b/testsuite/python/writevtf.py index a9e842ae6ad..96c0aba43b5 100644 --- a/testsuite/python/writevtf.py +++ b/testsuite/python/writevtf.py @@ -55,9 +55,10 @@ class CommonTests(ut.TestCase): system.bonded_inter.add( espressomd.interactions.FeneBond(k=1., d_r_max=10.0)) - system.part[0].add_bond((0, 1)) - system.part[0].add_bond((0, 2)) - system.part[0].add_bond((0, 3)) + p0 = system.part.by_id(0) + p0.add_bond((0, 1)) + p0.add_bond((0, 2)) + p0.add_bond((0, 3)) system.integrator.run(steps=0) @@ -101,10 +102,9 @@ def test_atoms(self): def test_vtf_pid_map(self): system = self.system types = self.types_to_write - if types == 'all': - ids = system.part[:].id - else: - ids = system.part[:].id[np.isin(system.part[:].type, types)] + ids = system.part.all().id + if types != 'all': + ids = ids[np.isin(system.part.all().type, types)] mapping_ref = dict(zip(ids, range(len(ids)))) mapping = espressomd.io.writer.vtf.vtf_pid_map(system, types=types) self.assertEqual(mapping, mapping_ref) diff --git a/testsuite/scripts/samples/test_MDAnalysisIntegration.py b/testsuite/scripts/samples/test_MDAnalysisIntegration.py index 69746436958..8703ee58344 100644 --- a/testsuite/scripts/samples/test_MDAnalysisIntegration.py +++ b/testsuite/scripts/samples/test_MDAnalysisIntegration.py @@ -36,18 +36,19 @@ class Sample(ut.TestCase): def test_universe(self): system = self.system + partcls = system.part.all() u = sample.u self.assertEqual(len(u.atoms), 10) np.testing.assert_equal(u.atoms.ids, np.arange(10) + 1) np.testing.assert_equal(u.atoms.types, sorted(5 * ['T0', 'T1'])) np.testing.assert_almost_equal( - u.atoms.charges, system.part[:].q, decimal=6) + u.atoms.charges, partcls.q, decimal=6) np.testing.assert_almost_equal( - u.atoms.positions, system.part[:].pos, decimal=6) + u.atoms.positions, partcls.pos, decimal=6) np.testing.assert_almost_equal( - u.atoms.velocities, system.part[:].v, decimal=6) + u.atoms.velocities, partcls.v, decimal=6) np.testing.assert_almost_equal( - u.atoms.forces, system.part[:].f, decimal=6) + u.atoms.forces, partcls.f, decimal=6) if __name__ == "__main__": diff --git a/testsuite/scripts/samples/test_visualization_elc.py b/testsuite/scripts/samples/test_visualization_elc.py index 1d04c264065..ac6a19ed073 100644 --- a/testsuite/scripts/samples/test_visualization_elc.py +++ b/testsuite/scripts/samples/test_visualization_elc.py @@ -38,7 +38,8 @@ class Sample(ut.TestCase): def test_dipole_moment(self): import espressomd.observables - obs = espressomd.observables.DipoleMoment(ids=self.system.part[:].id) + obs = espressomd.observables.DipoleMoment( + ids=self.system.part.all().id) dipm = obs.calculate() self.assertLess(dipm[2], 0, msg="charges moved in the wrong direction") # the dipole moment should be the strongest along the z-axis From dfcd36f95081d61aa8c444d818eaf3ae9e5edf59 Mon Sep 17 00:00:00 2001 From: Christoph Lohrmann Date: Wed, 8 Dec 2021 13:40:40 +0100 Subject: [PATCH 2/2] sphinxdoc --- doc/sphinx/inter_bonded.rst | 4 ++-- doc/sphinx/introduction.rst | 12 +++++----- doc/sphinx/io.rst | 2 +- doc/sphinx/particles.rst | 45 ++++++++++++++----------------------- doc/sphinx/running.rst | 4 ++-- doc/sphinx/system_setup.rst | 2 +- 6 files changed, 29 insertions(+), 40 deletions(-) diff --git a/doc/sphinx/inter_bonded.rst b/doc/sphinx/inter_bonded.rst index 88981ae259e..84e6458e32d 100644 --- a/doc/sphinx/inter_bonded.rst +++ b/doc/sphinx/inter_bonded.rst @@ -41,7 +41,7 @@ Note that alternatively to particle handles, the particle's ids can be used to setup bonded interactions. For example, to create a bond between the particles with the ids 12 and 43:: - system.part[12].add_bond((fene, 43)) + system.part.by_id(12).add_bond((fene, 43)) .. _Distance-dependent bonds: @@ -492,7 +492,7 @@ where ``softID`` identifies the soft particle and ``kappaV`` is a volumetric spring constant. Note that this ``volCons`` bond does not have a bond partner. It is added to a particle as follows:: - system.part[0].add_bond((volCons,)) + system.part.by_id(0).add_bond((volCons,)) The comma is needed to create a tuple containing a single item. diff --git a/doc/sphinx/introduction.rst b/doc/sphinx/introduction.rst index e10fae6acdb..60ba0ca4687 100644 --- a/doc/sphinx/introduction.rst +++ b/doc/sphinx/introduction.rst @@ -182,21 +182,21 @@ particles>` over all particles:: Internally, each particle is automatically assigned a unique numerical id by |es|. Note that in principle it is possible to explicitly set this particle id (if not in use already) on particle creation. -Using the id in square brackets, the respective particle can be accessed from the particle list:: +Using the id, the respective particle can be accessed from the particle list:: >>> system.part.add(id=3, pos=[2.1, 1.2, 3.3], type=0) - >>> system.part[3].pos = [1.0, 1.0, 2.0] - >>> print(system.part[3]) + >>> system.part.by_id(3).pos = [1.0, 1.0, 2.0] + >>> print(system.part.by_id(3).pos) [1.0, 1.0, 2.0] For larger simulation setups, explicit handling of numerical ids can quickly become confusing and is thus error-prone. We therefore highly recommend using :class:`~espressomd.particle_data.ParticleHandle` instead wherever possible. -:ref:`Properties of several particles` -can be accessed by using Python slices: :: +:ref:`Properties of all particles` +can be accessed via: :: - positions = system.part[:].pos + positions = system.part.all().pos .. rubric:: Interactions diff --git a/doc/sphinx/io.rst b/doc/sphinx/io.rst index d8d17b29f37..1c07a34362f 100644 --- a/doc/sphinx/io.rst +++ b/doc/sphinx/io.rst @@ -442,7 +442,7 @@ To read a PDB file containing a single frame:: # create particles and add bonds between them system.part.add(pos=np.array(chainA.positions, dtype=float)) for i in range(0, len(chainA) - 1): - system.part[i].add_bond((system.bonded_inter[0], system.part[i + 1].id)) + system.part.by_id(i).add_bond((system.bonded_inter[0], i + 1)) # visualize protein in 3D import espressomd.visualization visualizer = espressomd.visualization.openGLLive(system, bond_type_radius=[0.2]) diff --git a/doc/sphinx/particles.rst b/doc/sphinx/particles.rst index bee2d7789e8..f46a041d2b9 100644 --- a/doc/sphinx/particles.rst +++ b/doc/sphinx/particles.rst @@ -55,7 +55,7 @@ automatically. Alternatively, you can assign an id manually when adding them to The id provides an alternative way to access particles in the system. To retrieve the handle of the particle with id ``INDEX``, call:: - p = system.part[] + p = system.part.by_id() .. _Accessing particle properties: @@ -134,11 +134,11 @@ Exclusions do not apply to the short range part of electrostatics and magnetosta To create exclusions for particles pairs 0 and 1:: - system.part[0].add_exclusion(1) + system.part.by_id(0).add_exclusion(1) To delete the exclusion:: - system.part[0].delete_exclusion(1) + system.part.by_id(0).delete_exclusion(1) See :attr:`espressomd.particle_data.ParticleHandle.exclusions` @@ -346,30 +346,19 @@ to retrieve a particle slice: When adding several particles at once, a particle slice is returned instead of a particle handle. -- By slicing :py:attr:`espressomd.system.System.part` +- By calling :meth:`espressomd.particle_data.ParticleList.by_ids` - The :class:`~espressomd.particle_data.ParticleList` supports slicing - similarly to lists and NumPy arrays, however with the distinction that - particle slices can have gaps. + It is also possible to get a slice containing particles of specific ids, e.g.:: - Using a colon returns a slice containing all particles:: + system.part.by_ids([1, 4, 3]) - print(system.part[:]) - - To access particles with ids ranging from 0 to 9, use:: - - system.part[0:10] - - Note that, like in other cases in Python, the lower bound is inclusive and - the upper bound is non-inclusive. The length of the slice does not have to - be 10, it can be for example 2 if there are only 2 particles in the system - with an id between 0 and 9. - - It is also possible to get a slice containing particles of specific ids:: + would contain the particles with ids 1, 4, and 3 in that specific order. + +- By calling :meth:`espressomd.particle_data.ParticleList.all` - system.part[[1, 4, 3]] + You can get a slice containing all particles using:: - would contain the particles with ids 1, 4, and 3 in that specific order. + system.part.all() - By calling :meth:`espressomd.particle_data.ParticleList.select` @@ -381,7 +370,7 @@ to retrieve a particle slice: Properties of particle slices can be accessed just like with single particles. A list of all values is returned:: - print(system.part[:].q) + print(system.part.all().q) A particle slice can be iterated over, see :ref:`Iterating over particles and pairs of particles`. @@ -389,27 +378,27 @@ Setting properties of slices can be done by - supplying a *single value* that is assigned to each entry of the slice, e.g.:: - system.part[0:10].ext_force = [1, 0, 0] + system.part.by_ids(range(10)).ext_force = [1, 0, 0] - supplying an *array of values* that matches the length of the slice which sets each entry individually, e.g.:: - system.part[0:3].ext_force = [[1, 0, 0], [2, 0, 0], [3, 0, 0]] + system.partby_ids(range(3)).ext_force = [[1, 0, 0], [2, 0, 0], [3, 0, 0]] For list properties that have no fixed length like ``exclusions`` or ``bonds``, some care has to be taken. There, *single value* assignment also accepts lists/tuples just like setting the property of an individual particle. For example:: - system.part[0].exclusions = [1, 2] + system.part.by_id(0).exclusions = [1, 2] would both exclude short-range interactions of the particle pairs ``0 <-> 1`` and ``0 <-> 2``. Similarly, a list can also be assigned to each entry of the slice:: - system.part[2:4].exclusions = [0, 1] + system.part.by_ids(range(2,4)).exclusions = [0, 1] This would exclude interactions between ``2 <-> 0``, ``2 <-> 1``, ``3 <-> 0`` and ``3 <-> 1``. Now when it is desired to supply an *array of values* with individual values for each slice entry, the distinction can no longer be done by the length of the input, as slice length and input length can be equal. Here, the nesting level of the input is the distinctive criterion:: - system.part[2:4].exclusions = [[0, 1], [0, 1]] + system.part.by_ids(range(2,4)).exclusions = [[0, 1], [0, 1]] The above code snippet would lead to the same exclusions as the one before. The same accounts for the ``bonds`` property by interchanging the integer entries of the exclusion list with diff --git a/doc/sphinx/running.rst b/doc/sphinx/running.rst index 84846f80039..546b7520158 100644 --- a/doc/sphinx/running.rst +++ b/doc/sphinx/running.rst @@ -258,7 +258,7 @@ force/torque observable used in the custom convergence criterion. Since these two properties can be cast to boolean values, they can be used as masks to remove forces/torques that are ignored by the integrator:: - particles = system.part[:] + particles = system.part.all() max_force = np.max(np.linalg.norm(particles.f * np.logical_not(particles.fix), axis=1)) max_torque = np.max(np.linalg.norm(particles.torque_lab * np.logical_not(particles.rotation), axis=1)) @@ -273,7 +273,7 @@ The correct forces need to be re-calculated after running the integration:: p2 = system.part.add(pos=[0, 0, 0.1], type=1) p2.vs_auto_relate_to(p1) system.integrator.set_steepest_descent(f_max=800, gamma=1.0, max_displacement=0.01) - while convergence_criterion(system.part[:].f): + while convergence_criterion(system.part.all().f): system.integrator.run(10) system.integrator.run(0, recalc_forces=True) # re-calculate forces from virtual sites system.integrator.set_vv() diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index 61958be23b4..84f9569606d 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -73,7 +73,7 @@ Some variables like or are no longer directly available as attributes. In these cases they can be easily derived from the corresponding Python objects like:: - n_part = len(system.part[:].pos) + n_part = len(system.part) or by calling the corresponding ``get_state()`` methods like::