New particle access syntax (#4402)

Fixes #4389 

Description of changes:
- The `[]`-operator on `system.part` was removed
- An individual particle can be accessed by id using `system.part.by_id(id)`, where `id` is the numerical particle id. This returns a `ParticleHandle`
- A slice of particles can be accessed via an iterable of ids with `system.part.by_ids(id_list)`, which returns a `ParticleSlice`
- `system.part.all()` returns a slice containing all particles
- The order in which particles will be provided when iterating over a slice will stay fixed, once the slice has been instanced. I.e, `for pos, vel in zip(system.part.all().pos, system.part.all().v)` will produce the expected result.
- Selection by criterion (either using `key=value` or a lambda) can be done with the existing `system.part.select()`

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.
 

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<Interacting with groups of particles>`
-can be accessed by using Python slices: ::
+:ref:`Properties of all particles<Interacting with groups of particles>`
+can be accessed via: ::
 
-    positions = system.part[:].pos
+    positions = system.part.all().pos
 
 .. rubric:: Interactions
 

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])

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[<INDEX>]
+    p = system.part.by_id(<INDEX>)
 
 .. _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,35 +370,35 @@ 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`.
 
 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

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()

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::
 

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,15 +303,15 @@
     "\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",
     "    while i < sd_params['max_steps'] // sd_params['emstep']:\n",
     "        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,

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,

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,

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,