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:
espressomd.particle_data.ParticleHandle
provides access to a single particle in the simulation.espressomd.particle_data.ParticleList
provides access to all particles in the simulationespressomd.particle_data.ParticleSlice
provides access to a subset of particles in the simulation identified by a list of ids.
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)
4.8. Create particular particle configurations¶
4.8.1. Setting up polymer chains¶
from espressomd.polymer import create_polymer
A function that allows to create a number of polymers and polyelectrolytes.
See espressomd.polymer.create_polymer()
for a detailed list of
arguments.
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.
The polymer can be created using several different random walk modes (via the parameter modes
):
- (Random walk)
mode = 1
The monomers are randomly placed by a random walk with a step size ofbond_length
.- (Pruned self-avoiding walk)
mode = 2
The position of a monomer is randomly chosen within a distance to the previous monomer. If the position is closer to another particle thanshield
, the attempt is repeated up tomax_tries
times. Note, that this is not a real self-avoiding random walk, as the particle distribution is not the same. If you want a real self-avoiding walk, use the mode 0. However, this mode is several orders of magnitude faster than a true self-avoiding random walk, especially for long chains.- (Self-avoiding random walk)
mode = 0
The positions of the monomers are chosen as in the plain random walk. However, if this results in a chain that has a monomer that is closer to another particle thanshield
, a new attempt of setting up the whole chain is done, up tomax_tries
times.
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\). 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.
See espressomd.diamond.Diamond()
for more details.
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
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,
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.
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.Repeat the previous step with more virtual sites, if desired.
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
andespressomd.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 in the grand canonical ensemble, or other purposes, when particles of certain types are created and deleted frequently is often of interest. Particle ids can be stored in a map for each individual type and so random ids of particles of a certain type can be drawn.
import espressomd
system = espressomd.System()
system.setup_type_map([_type])
system.find_particle(_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.
When using the
keyword find_particle
and a particle type, the command will return a randomly
chosen particle id, for a particle of the given type. 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. Therefore calling select()
is slow compared to using espressomd.system.System.number_of_particles()
which 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].
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 Open MPI
node, i.e. n_nodes
= 1.