4. Setting up particles

4.1. Overview of the relevant Python classes

For understanding this chapter, it is helpful to be aware of the Python classes provided by ESPResSo to interact with particles:

in almost no case have these classes to be instantiated explicitly by the user. Rather, access is provided via the espressomd.system.System.part attribute. The details are explained in the following sections.

4.2. Adding particles

In order to add particles to the system, call espressomd.particle_data.ParticleList.add():

import espressomd
system = espressomd.System()
system.part.add(pos=[1.0, 1.0, 1.0], id=0, type=0)

This command adds a single particle to the system with properties given as arguments. The pos property is required, the id can be omitted, in which case it is chosen automatically. All available particle properties are members of espressomd.particle_data.ParticleHandle.

It is also possible to add several particles at once:

import numpy as np
system.part.add(pos=np.random.random((10, 3) * box_length))

Furthermore, the espressomd.particle_data.ParticleList.add() method returns the added particle(s):

tracer = system.part.add(pos=(0, 0, 0))
print(tracer.pos)

Note that the instance of espressomd.particle_data.ParticleHandle returned by espressomd.particle_data.ParticleList.add() are handles for the live particles in the simulation, rather than offline copies. Changing their properties will affect the simulation.

4.3. Accessing particle properties

Particle properties can be accessed like a class member.

To access property PROPERTY of the particle with index INDEX:

system.part[<INDEX>].<PROPERTY>

For example, to print the current position of the particle with index 0 in the system, call:

print(system.part[0].pos)

Similarly, the position can be set:

system.part[0].pos = (1, 2.5, 3)

4.3.1. Vectorial properties

For vectorial particle properties, component-wise manipulation like system.part[0].pos[0] = 1 or in-place operators like += or *= are not allowed and result in an error. This behavior is inherited, so the same applies to a after a = system.part[0].pos. If you want to use a vectorial property for further calculations, you should explicitly make a copy e.g. via a = numpy.copy(system.part[0].pos).

4.4. Interacting with groups of particles

The espressomd.particle_data.ParticleList support slicing similarly to lists and NumPy arrays. To access all existing particles, use a colon:

print(sysstem.part[:].pos)
system.part[:].q = 0

To access particles with indices ranging from 0 to 9, use:

system.part[0:10].pos

Note that, like in other cases in Python, the lower bound is inclusive and the upper bound is non-inclusive. Setting 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]

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

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]

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]

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

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 the tuple (bond, partners).

You can select a subset of particles via using the select method. For example you can obtain a list of particles with charge -1 via using

system.part.select(q=-1)

For further information on how to use selections see espressomd.particle_data.ParticleList.select().

4.5. Deleting particles

Particles can be easily deleted in Python using particle ids or ranges of particle ids. For example, to delete all particles with particle index greater than 10, run:

system.part[10:].remove()

To delete all particles, use:

system.part.clear()

4.6. Iterating over particles and pairs of particles

You can iterate over all particles or over a subset of particles as follows:

for p in system.part:
    print(p.pos)

for p in system.part[0:10]:
    print(p.pos)

You can iterate over all pairs of particles using:

for pair in system.part.pairs():
    print(pair[0].id, pair[1].id)

4.7. Exclusions

Particles can have an exclusion list of all other particles where non-bonded interactions are ignored. This is typically used in atomistic simulations, where nearest and next nearest neighbor interactions along the chain have to be omitted since they are included in the bonding potentials. Exclusions do not apply to the short range part of electrostatics and magnetostatics methods, e.g. to P3M.

system.part[0].add_exclusion(1)

Create exclusions for particles pairs 0 and 1.

To delete the exclusion, simply use

system.part[0].delete_exclusion(1)

See espressomd.particle_data.ParticleHandle.exclusions

4.8. Create particular particle configurations

4.8.1. Setting up polymer chains

If you want to have polymers in your system, you can use the function espressomd.polymer.positions() to determine suitable positions.

Required arguments are the desired number of polymers n_polymers, the number of monomers per polymer chain beads_per_chain, and the parameter bond_length, which determines the distance between adjacent monomers within the polymer chains. Determining suitable particle positions pseudo-randomly requires the use of a pseudo-random number generator, which has to be seeded. This seed is therefore also a mandatory parameter.

The function espressomd.polymer.positions() returns a three-dimensional numpy array, namely a list of polymers containing the positions of monomers (x, y, z). A quick example of how to set up polymers:

import espressomd
from espressomd import polymer

system = espressomd.System([50, 50, 50])
polymers = polymer.positions(n_polymers=10,
                             beads_per_chain=25,
                             bond_length=0.9, seed=23)
for p in polymers:
    for i, m in enumerate(p):
       id = len(system.part)
       system.part.add(id=id, pos=m)
       if i > 0:
           system.part[id].add_bond((<BOND_TYPE>, id - 1))

If there are constraints present in your system which you want to be taken into account when creating the polymer positions, you can set the optional boolean parameter respect_constraint=True. To simulate excluded volume while drawing the polymer positions, a minimum distance between all particles can be set via min_distance. This will also respect already existing particles in the system. Both when setting respect_constraints and choosing a min_distance trial positions are pseudo-randomly chosen and only accepted if the requested requirement is fulfilled. Otherwise, a new attempt will be made, up to max_tries times per monomer and if this fails max_tries per polymer. The default value is max_tries=1000. Depending on the total number of beads and constraints, this value may need to be adapted. If determining suitable polymer positions within this limit fails, a runtime error is thrown.

Note that the distance between adjacent monomers during the course of the simulation depends on the applied potentials. For fixed bond length please refer to the Rattle Shake algorithm[And83]. The algorithm is based on Verlet algorithm and satisfy internal constraints for molecular models with internal constraints, using Lagrange multipliers.

4.8.2. Setting up diamond polymer networks

from espressomd import diamond

Creates a diamond-structured polymer network with 8 tetra-functional nodes connected by \(2*8\) polymer chains of length (MPC) in a unit cell of length \(a\). The diamond command creates 16*MPC+8 many particles which are connected via the provided bond type (the term plus 8 stems from adding 8 nodes which are connecting the chains). Chain monomers are placed at a mutual distance along the vector connecting network nodes. The polymer is created starting from particle ID 0. Nodes are assigned type 0, monomers (both charged and uncharged) are type 1 and counterions type 2. For inter-particle bonds interaction \(0\) is taken which must be a two-particle bond.

Diamond-like polymer network with MPC=15.

Diamond-like polymer network with MPC=15.

See espressomd.diamond.Diamond for more details. For simulating compressed or stretched gels the function espressomd.system.System.change_volume_and_rescale_particles() may be used.

4.9. Virtual sites

Virtual sites are particles, the positions and velocities of which are not obtained by integrating an equation of motion. Rather, their coordinates are obtained from the position (and orientation) of one or more other particles. In this way, rigid arrangements of particles can be constructed and a particle can be placed in the center of mass of a set of other particles. Virtual sites can interact with other particles in the system by means of interactions. Forces are added to them according to their respective particle type. Before the next integration step, the forces accumulated on a virtual site are distributed back to those particles, from which the virtual site was derived.

There are different schemes for virtual sites, described in the following sections. To switch the active scheme, the attribute espressomd.system.System.virtual_sites of the system class can be used:

import espressomd
from espressomd.virtual_sites import VirtualSitesOff, VirtualSitesRelative

s = espressomd.System()
s.virtual_sites = VirtualSitesRelative(have_velocity=True, have_quaternion=False)
# or
s.virtual_sites = VirtualSitesOff()

By default, espressomd.virtual_sites.VirtualSitesOff is selected. This means that virtual particles are not touched during integration. The have_velocity parameter determines whether or not the velocity of virtual sites is calculated, which carries a performance cost. The have_quaternion parameter determines whether the quaternion of the virtual particle is updated (useful in combination with the espressomd.particle_data.ParticleHandle.vs_quat property of the virtual particle which defines the orientation of the virtual particle in the body fixed frame of the related real particle.

4.9.1. Rigid arrangements of particles

The relative implementation of virtual sites allows for the simulation of rigid arrangements of particles. It can be used, for extended dipoles and raspberry-particles, but also for more complex configurations. Position and velocity of a virtual site are obtained from the position and orientation of exactly one non-virtual particle, which has to be placed in the center of mass of the rigid body. Several virtual sites can be related to one and the same non-virtual particle. The position of the virtual site is given by

\[\vec{x_v} =\vec{x_n} +O_n (O_v \vec{E_z}) d,\]

where \(\vec{x_n}\) is the position of the non-virtual particle, \(O_n\) is the orientation of the non-virtual particle, \(O_v\) denotes the orientation of the vector \(\vec{x_v}-\vec{x_n}\) with respect to the non-virtual particles body fixed frame and \(d\) the distance between virtual and non-virtual particle. In words: The virtual site is placed at a fixed distance from the non-virtual particle. When the non-virtual particle rotates, the virtual sites rotates on an orbit around the non-virtual particles center.

To use this implementation of virtual sites, activate the feature VIRTUAL_SITES_RELATIVE. Furthermore, an instance of espressomd.virtual_sites.VirtualSitesRelative has to be set as the active virtual sites scheme (see above). To set up a virtual site,

  1. Place the particle to which the virtual site should be related. It needs to be in the center of mass of the rigid arrangement of particles you create. Let its particle id be n.

  2. Place a particle at the desired relative position, make it virtual and relate it to the first particle:

    p = system.part.add(pos=(1, 2, 3))
p.vs_auto_relate_to(<ID>)

    where <ID> is the id of the central particle. This will also set the espressomd.particle_data.ParticleHandle.virtual attribute on the particle to 1.

  3. Repeat the previous step with more virtual sites, if desired.

  4. To update the positions of all virtual sites, call:

    system.integrator.run(0, recalc_forces=True)

Please note:

  • The relative position of the virtual site is defined by its distance from the non-virtual particle, the id of the non-virtual particle and a quaternion which defines the vector from non-virtual particle to virtual site in the non-virtual particles body-fixed frame. This information is saved in the virtual site’s espressomd.particle_data.ParticleHandle.vs_relative attribute. Take care, not to overwrite it after using vs_auto_relate.

  • Virtual sites can not be placed relative to other virtual sites, as the order in which the positions of virtual sites are updated is not guaranteed. Always relate a virtual site to a non-virtual particle placed in the center of mass of the rigid arrangement of particles.

  • In case you know the correct quaternions, you can also setup a virtual site using its espressomd.particle_data.ParticleHandle.vs_relative and espressomd.particle_data.ParticleHandle.virtual attributes.

  • In a simulation on more than one CPU, the effective cell size needs to be larger than the largest distance between a non-virtual particle and its associated virtual sites. To this aim, when running on more than one core, you need to set the system’s espressomd.system.System.min_global_cut attribute to this largest distance. An error is generated when this requirement is not met.

  • If the virtual sites represent actual particles carrying a mass, the inertia tensor of the non-virtual particle in the center of mass needs to be adapted.

  • The presence of rigid bodies constructed by means of virtual sites adds a contribution to the pressure and stress tensor.

4.9.2. Inertialess lattice Boltzmann tracers

espressomd.virtual_sites.VirtualSitesInertialessTracers

When this implementation is selected, the virtual sites follow the motion of a lattice Boltzmann fluid (both, CPU and GPU). This is achieved by integrating their position using the fluid velocity at the virtual sites’ position. Forces acting on the virtual sites are directly transferred as force density onto the lattice Boltzmann fluid, making the coupling free of inertia. The feature stems from the implementation of the Immersed Boundary Method for soft elastic objects, but can be used independently.

For correct results, the LB thermostat has to be deactivated for virtual sites:

system.thermostat.set_lb(kT=0, act_on_virtual=False)

Please note that the velocity attribute of the virtual particles does not carry valid information for this virtual sites scheme.

4.10. Particle number counting feature

Note

Do not use these methods with the espressomd.collision_detection module since the collision detection may create or delete particles without the particle number counting feature being aware of this. Therefore also the espressomd.reaction_ensemble module may not be used with the collision detection.

Knowing the number of particles of a certain type in simulations where particle numbers can fluctuate is of interest. Particle ids can be stored in a map for each individual type:

import espressomd
system = espressomd.System()
system.setup_type_map([_type])
system.number_of_particles(_type)

If you want to keep track of particle ids of a certain type you have to initialize the method by calling

system.setup_type_map([_type])

After that will keep track of particle ids of that type. Keeping track of particles of a given type is not enabled by default since it requires memory. The keyword number_of_particles as argument will return the number of particles which have the given type. For counting the number of particles of a given type you could also use espressomd.particle_data.ParticleList.select()

import espressomd
system = espressomd.System()
...
number_of_particles = len(system.part.select(type=type))

However calling select(type=type) results in looping over all particles which is slow. In contrast, espressomd.system.System.number_of_particles() directly can return the number of particles with that type.

4.11. Self-propelled swimmers

Note

If you are using this feature, please cite [dGMM+16].

See also

espressomd.particle_data.ParticleHandle.swimming

4.11.1. Langevin swimmers

import espressomd

system = espressomd.System()

system.part.add(id=0, pos=[1, 0, 0], swimming={'f_swim': 0.03})

This enables the particle to be self-propelled in the direction determined by its quaternion. For setting the particle’s quaternion see espressomd.particle_data.ParticleHandle.quat. The self-propulsion speed will relax to a constant velocity, that is specified by v_swim. Alternatively it is possible to achieve a constant velocity by imposing a constant force term f_swim that is balanced by friction of a (Langevin) thermostat. The way the velocity of the particle decays to the constant terminal velocity in either of these methods is completely determined by the friction coefficient. You may only set one of the possibilities v_swim or f_swim as you cannot relax to constant force and constant velocity at the same time. Note that there is no real difference between v_swim and f_swim, since the latter may always be chosen such that the same terminal velocity is achieved for a given friction coefficient.

4.11.2. Lattice Boltzmann (LB) swimmers

import espressomd

system = espressomd.System()

system.part.add(id=1, pos=[2, 0, 0], rotation=[1, 1, 1], swimming={
    'f_swim': 0.01, 'mode': 'pusher', 'dipole_length': 2.0, 'rotational_friction': 20})

For an explanation of the parameters v_swim and f_swim see the previous item. In lattice Boltzmann self-propulsion is less trivial than for regular MD, because the self-propulsion is achieved by a force-free mechanism, which has strong implications for the far-field hydrodynamic flow field induced by the self-propelled particle. In ESPResSo only the dipolar component of the flow field of an active particle is taken into account. This flow field can be generated by a pushing or a pulling mechanism, leading to change in the sign of the dipolar flow field with respect to the direction of motion. You can specify the nature of the particle’s flow field by using one of the modes: pusher or puller. You will also need to specify a dipole_length which determines the distance of the source of propulsion from the particle’s center. Note that you should not put this distance to zero; ESPResSo (currently) does not support mathematical dipole flow fields. The key rotational_friction can be used to set the friction that causes the orientation of the particle to change in shear flow. The torque on the particle is determined by taking the cross product of the difference between the fluid velocity at the center of the particle and at the source point and the vector connecting the center and source.

You may ask: “Why are there two methods v_swim and f_swim for the self-propulsion using the lattice Boltzmann algorithm?” The answer is straightforward. When a particle is accelerating, it has a monopolar flow-field contribution which vanishes when it reaches its terminal velocity (for which there will only be a dipolar flow field). The major difference between the above two methods is that with v_swim the flow field only has a monopolar moment and only while the particle is accelerating. As soon as the particle reaches a constant speed (given by v_swim) this monopolar moment is gone and the flow field is zero! In contrast, f_swim always, i.e., while accelerating and while swimming at constant force possesses a dipolar flow field.

Warning

Please note that even though swimming is interoperable with the CPU version of LB it is only supported on one MPI rank, i.e. n_nodes = 1.