3. Setting up the system¶
3.1. Setting global variables in Python¶
The global variables in Python are controlled via the
espressomd.system.System
class.
Global system variables can be read and set in Python simply by accessing the
attribute of the corresponding Python object. Those variables that are already
available in the Python interface are listed in the following. Note that for the
vectorial properties box_l
and periodicity
, componentwise manipulation
like system.box_l[0] = 1
or inplace operators like +=
or *=
are not
allowed and result in an error. This behavior is inherited, so the same applies
to a
after a = system.box_l
. If you want to use a vectorial property
for further calculations, you should explicitly make a copy e.g. via
a = numpy.copy(system.box_l)
.

(float[3]) Simulation box lengths of the cuboid box used by ESPResSo. Note that if you change the box length during the simulation, the folded particle coordinates will remain the same, i.e., the particle stay in the same image box, but at the same relative position in their image box. If you want to scale the positions, use the command
change_volume_and_rescale_particles()
. 
(int[3]) Specifies periodicity for the three directions. ESPResSo can be instructed to treat some dimensions as nonperiodic. By default ESPResSo assumes periodicity in all directions which equals setting this variable to
[True, True, True]
. A dimension is specified as nonperiodic via setting the periodicity variable for this dimension toFalse
. E.g. Periodicity only in zdirection is obtained by[False, False, True]
. Caveat: Be aware of the fact that making a dimension nonperiodic does not hinder particles from leaving the box in this direction. In this case for keeping particles in the simulation box a constraint has to be set. 
(float) Time step for MD integration.

(float) The simulation time.

(float) Minimal total cutoff for real space. Effectively, this plus the
skin
is the minimally possible cell size. ESPResSo typically determines this value automatically, but some algorithms, virtual sites, require you to specify it manually. 
readonly Maximal cutoff of bonded real space interactions.

readonly Maximal cutoff of bonded real space interactions.
3.1.1. Accessing module states¶
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)
or by calling the corresponding get_state()
methods like:
temperature = system.thermostat.get_state()[0]['kT']
gamma = system.thermostat.get_state()[0]['gamma']
gamma_rot = system.thermostat.get_state()[0]['gamma_rotation']
3.2. Cellsystems¶
This section deals with the flexible particle data organization of ESPResSo. Due to different needs of different algorithms, ESPResSo is able to change the organization of the particles in the computer memory, according to the needs of the used algorithms. For details on the internal organization, refer to section Internal particle organization.
3.2.1. Global properties¶
The properties of the cell system can be accessed by
espressomd.system.System.cell_system
:
(int[3]) 3D node grid for real space domain decomposition (optional, if unset an optimal set is chosen automatically). The domain decomposition can be visualized with
samples/visualization_cellsystem.py
.(float) Skin for the Verlet list. This value has to be set, otherwise the simulation will not start.
Details about the cell system can be obtained by espressomd.system.System.cell_system.get_state()
:
cell_grid
Dimension of the inner cell grid.
cell_size
Boxlength of a cell.
local_box_l
Local simulation box length of the nodes.
max_cut
Maximal cutoff of real space interactions.
n_nodes
Number of nodes.
type
The current type of the cell system.
verlet_reuse
Average number of integration steps the Verlet list is reused.
3.2.2. Domain decomposition¶
Invoking set_domain_decomposition()
selects the domain decomposition cell scheme, using Verlet lists
for the calculation of the interactions. If you specify use_verlet_lists=False
, only the
domain decomposition is used, but not the Verlet lists.
system = espressomd.System(box_l=[1, 1, 1])
system.cell_system.set_domain_decomposition(use_verlet_lists=True)
The domain decomposition cellsystem is the default system and suits most applications with short ranged interactions. The particles are divided up spatially into small compartments, the cells, such that the cell size is larger than the maximal interaction range. In this case interactions only occur between particles in adjacent cells. Since the interaction range should be much smaller than the total system size, leaving out all interactions between nonadjacent cells can mean a tremendous speedup. Moreover, since for constant interaction range, the number of particles in a cell depends only on the density. The number of interactions is therefore of the order N instead of order \(N^2\) if one has to calculate all pair interactions.
3.2.3. Nsquared¶
Invoking set_n_square()
selects the very primitive nsquared cellsystem, which calculates
the interactions for all particle pairs. Therefore it loops over all
particles, giving an unfavorable computation time scaling of
\(N^2\). However, algorithms like MMM1D or the plain Coulomb
interaction in the cell model require the calculation of all pair
interactions.
system = espressomd.System(box_l=[1, 1, 1])
system.cell_system.set_n_square()
In a multiple processor environment, the nsquared cellsystem uses a simple particle balancing scheme to have a nearly equal number of particles per CPU, \(n\) nodes have \(m\) particles, and \(pn\) nodes have \(m+1\) particles, such that \(n \cdot m + (p  n) \cdot (m + 1) = N\), the total number of particles. Therefore the computational load should be balanced fairly equal among the nodes, with one exception: This code always uses one CPU for the interaction between two different nodes. For an odd number of nodes, this is fine, because the total number of interactions to calculate is a multiple of the number of nodes, but for an even number of nodes, for each of the \(p1\) communication rounds, one processor is idle.
E.g. for 2 processors, there are 3 interactions: 00, 11, 01. Naturally, 00 and 11 are treated by processor 0 and 1, respectively. But the 01 interaction is treated by node 1 alone, so the workload for this node is twice as high. For 3 processors, the interactions are 00, 11, 22, 01, 12, 02. Of these interactions, node 0 treats 00 and 02, node 1 treats 11 and 01, and node 2 treats 22 and 12.
Therefore it is highly recommended that you use nsquared only with an odd number of nodes, if with multiple processors at all.
3.3. Thermostats¶
The thermostat can be controlled by the class espressomd.thermostat.Thermostat
.
The different thermostats available in ESPResSo will be described in the following
subsections.
You may combine different thermostats at your own risk by turning them on
one by one. The list of active thermostats can be cleared at any time with
system.thermostat.turn_off()
.
Not all combinations of thermostats are allowed, though (see
espressomd.thermostat.AssertThermostatType()
for details).
Some integrators only work with a specific thermostat and throw an
error otherwise. Note that there is only one temperature for all
thermostats, although for some thermostats like the Langevin thermostat,
particles can be assigned individual temperatures.
Since ESPResSo does not enforce a particular unit system, it cannot know about the current value of the Boltzmann constant. Therefore, when specifying the temperature of a thermostat, you actually do not define the temperature, but the value of the thermal energy \(k_B T\) in the current unit system (see the discussion on units, Section On units).
All thermostats have a seed
argument that controls the state of the random
number generator (Philox Counterbased RNG). This seed is required on first
activation of a thermostat, unless stated otherwise. It can be omitted in
subsequent calls of the method that activates the same thermostat. The random
sequence also depends on the thermostats counters that are
incremented after each integration step.
3.3.1. Langevin thermostat¶
In order to activate the Langevin thermostat the member function
set_langevin()
of the thermostat
class espressomd.thermostat.Thermostat
has to be invoked.
Best explained in an example:
import espressomd
system = espressomd.System(box_l=[1, 1, 1])
system.thermostat.set_langevin(kT=1.0, gamma=1.0, seed=41)
As explained before the temperature is set as thermal energy \(k_\mathrm{B} T\).
The Langevin thermostat is based on an extension of Newton’s equation of motion to
Here, \(f_i\) are all deterministic forces from interactions, \(\gamma\) the bare friction coefficient and \(\eta\) a random, “thermal” force. The friction term accounts for dissipation in a surrounding fluid whereas the random force mimics collisions of the particle with solvent molecules at temperature \(T\) and satisfies
(\(<\cdot>\) denotes the ensemble average and \(\alpha,\beta\) are spatial coordinates).
In the ESPResSo implementation of the Langevin thermostat, the additional terms only enter in the force calculation. This reduces the accuracy of the Velocity Verlet integrator by one order in \(dt\) because forces are now velocitydependent.
The random process \(\eta(t)\) is discretized by drawing an uncorrelated random number \(\overline{\eta}\) for each component of all the particle forces. The distribution of \(\overline{\eta}\) is uniform and satisfies
If the feature ROTATION
is compiled in, the rotational degrees of freedom are
also coupled to the thermostat. If only the first two arguments are
specified then the friction coefficient for the rotation is set to the
same value as that for the translation.
A separate rotational friction coefficient can be set by inputting
gamma_rotate
. The two options allow one to switch the translational and rotational
thermalization on or off separately, maintaining the frictional behavior. This
can be useful, for instance, in high Péclet number active matter systems, where
one only wants to thermalize only the rotational degrees of freedom and
translational motion is affected by the selfpropulsion.
The keywords gamma
and gamma_rotate
can be specified as a scalar,
or, with feature PARTICLE_ANISOTROPY
compiled in, as the three eigenvalues
of the respective friction coefficient tensor. This is enables the simulation of
the anisotropic diffusion of anisotropic colloids (rods, etc.).
Using the Langevin thermostat, it is possible to set a temperature and a
friction coefficient for every particle individually via the feature
LANGEVIN_PER_PARTICLE
. Consult the reference of the part
command
(chapter Setting up particles) for information on how to achieve this.
3.3.2. LatticeBoltzmann thermostat¶
The LatticeBoltzmann thermostat acts similar to the Langevin thermostat in that the governing equation for particles is
where \(u(x,t)\) is the fluid velocity at position \(x\) and time \(t\). To preserve momentum, an equal and opposite friction force and random force act on the fluid.
Numerically the fluid velocity is determined from the latticeBoltzmann node velocities by interpolating as described in Interpolating velocities. The backcoupling of friction forces and noise to the fluid is also done by distributing those forces amongst the nearest LB nodes. Details for both the interpolation and the force distribution can be found in [ADunweg99] and [DunwegL08].
The LB fluid can be used to thermalize particles, while also including their hydrodynamic interactions.
The LB thermostat expects an instance of either espressomd.lb.LBFluid
or espressomd.lb.LBFluidGPU
.
Temperature is set via the kT
argument of the LB fluid.
The magnitude of the frictional coupling can be adjusted by the
parameter gamma
. To enable the LB thermostat, use:
import espressomd
import espressomd.lb
system = espressomd.System(box_l=[1, 1, 1])
lbf = espressomd.lb.LBFluid(agrid=1, dens=1, visc=1, tau=0.01)
system.actors.add(lbf)
system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1.5)
No other thermostatting mechanism is necessary then. Please switch off any other thermostat before starting the LB thermostatting mechanism.
The LBM implementation provides a fully thermalized LB fluid, all nonconserved modes, including the pressure tensor, fluctuate correctly according to the given temperature and the relaxation parameters. All fluctuations can be switched off by setting the temperature to 0.
Note
Coupling between LB and MD only happens if the LB thermostat is set with a \(\gamma \ge 0.0\).
3.3.3. Dissipative Particle Dynamics (DPD)¶
The DPD thermostat adds friction and noise to the particle dynamics like the Langevin thermostat, but these are not applied to every particle individually but instead encoded in a dissipative interaction between particles [SDunwegK03].
To realize a complete DPD fluid model in ESPResSo, three parts are needed: the DPD thermostat, which controls the temperate, a dissipative interaction between the particles that make up the fluid, see DPD interaction, and a repulsive conservative force, see Hat interaction.
The temperature is set via
espressomd.thermostat.Thermostat.set_dpd()
which takes kT
and seed
as arguments.
The friction coefficients and cutoff are controlled via the DPD interaction on a per typepair basis. For details see there.
The friction (dissipative) and noise (random) term are coupled via the fluctuationdissipation theorem. The friction term is a function of the relative velocity of particle pairs. The DPD thermostat is better for dynamics than the Langevin thermostat, since it mimics hydrodynamics in the system.
As a conservative force any interaction potential can be used, see Isotropic nonbonded interactions. A common choice is a force ramp which is implemented as Hat interaction.
A complete example of setting up a DPD fluid and running it
to sample the equation of state can be found in /samples/dpd.py
.
When using a LennardJones interaction, \({r_\mathrm{cut}} = 2^{\frac{1}{6}} \sigma\) is a good value to choose, so that the thermostat acts on the relative velocities between nearest neighbor particles. Larger cutoffs including next nearest neighbors or even more are unphysical.
Boundary conditions for DPD can be introduced by adding the boundary
as a particle constraint, and setting a velocity and a type on it, see
espressomd.constraints.Constraint
. Then a
DPD interaction with the type can be defined, which acts as a
boundary condition.
3.3.4. Isotropic NPT thermostat¶
This feature allows to simulate an (on average) homogeneous and isotropic system in the NPT ensemble.
In order to use this feature, NPT
has to be defined in the myconfig.hpp
.
Activate the NPT thermostat with the command set_npt()
and setup the integrator for the NPT ensemble with set_isotropic_npt()
.
For example:
import espressomd
system = espressomd.System(box_l=[1, 1, 1])
system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0, seed=41)
system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0)
For an explanation of the algorithm involved, see Isotropic NPT integrator.
Be aware that this feature is neither properly examined for all systems nor is it maintained regularly. If you use it and notice strange behavior, please contribute to solving the problem.
3.3.5. Brownian thermostat¶
Brownian thermostat is a formal name of a thermostat enabling the Brownian Dynamics feature (see [Sch10]) which implies a propagation scheme involving systematic and thermal parts of the classical ErmakMcCammom’s (see [EM78]) Brownian Dynamics. Currently it is implemented without hydrodynamic interactions, i.e. with a diagonal diffusion tensor. The hydrodynamic interactions feature will be available later as a part of the present Brownian Dynamics or implemented separately within the Stokesian Dynamics.
In order to activate the Brownian thermostat, the member function
set_brownian
of the thermostat
class espressomd.thermostat.Thermostat
has to be invoked.
The system integrator should be also changed.
Best explained in an example:
import espressomd
system = espressomd.System(box_l=[1, 1, 1])
system.thermostat.set_brownian(kT=1.0, gamma=1.0, seed=41)
system.integrator.set_brownian_dynamics()
where gamma
(hereinafter \(\gamma\)) is a viscous friction coefficient.
In terms of the Python interface and setup, the Brownian thermostat is very
similar to the Langevin thermostat. The feature
BROWNIAN_PER_PARTICLE
is used to control the perparticle
temperature and the friction coefficient setup. The major differences are
its internal integrator implementation and other temporal constraints.
The integrator is still a symplectic Velocity Verletlike one.
It is implemented via a viscous drag part and a random walk of both the position and
velocity. Due to a nature of the Brownian Dynamics method, its time step \(\Delta t\)
should be large enough compared to the relaxation time
\(m/\gamma\) where \(m\) is the particle mass.
This requirement is just a conceptual one
without specific implementation technical restrictions.
Note that with all similarities of
Langevin and Brownian Dynamics, the Langevin thermostat temporal constraint
is opposite. A velocity is restarting from zero at every step.
Formally, the previous step velocity at the beginning of the the \(\Delta t\) interval
is dissipated further
and does not contribute to the end one as well as to the positional random walk.
Another temporal constraint
which is valid for both Langevin and Brownian Dynamics: conservative forces
should not change significantly over the \(\Delta t\) interval.
The viscous terminal velocity \(\Delta v\) and corresponding positional step \(\Delta r\) are fully driven by conservative forces \(F\):
A positional random walk variance of each coordinate \(\sigma_p^2\) corresponds to a diffusion within the Wiener process:
Each velocity component random walk variance \(\sigma_v^2\) is defined by the heat component:
Note that the velocity random walk is propagated from zero at each step.
A rotational motion is implemented similarly. Note: the rotational Brownian dynamics implementation is compatible with particles which have the isotropic moment of inertia tensor only. Otherwise, the viscous terminal angular velocity is not defined, i.e. it has no constant direction over the time.
3.3.6. Stokesian thermostat¶
Note
Requires STOKESIAN_DYNAMICS
external feature, enabled with
DWITH_STOKESIAN_DYNAMICS=ON
.
In order to thermalize a Stokesian Dynamics simulation, the SD thermostat needs to be activated via:
import espressomd
system = espressomd.System(box_l=[1.0, 1.0, 1.0])
system.periodicity = [False, False, False]
system.time_step = 0.01
system.cell_system.skin = 0.4
system.part.add(pos=[0, 0, 0], rotation=[1, 0, 0], ext_force=[0, 0, 1])
system.thermostat.set_stokesian(kT=1.0, seed=43)
system.integrator.set_stokesian_dynamics(viscosity=1.0, radii={0: 1.0})
system.integrator.run(100)
where kT
denotes the desired temperature of the system, and seed
the
seed for the random number generator.
3.4. CUDA¶
CudaInitHandle()
command can be used to choose the GPU for all subsequent
GPUcomputations. Note that due to driver limitations, the GPU cannot be
changed anymore after the first GPUusing command has been issued, for
example lbfluid
. If you do not choose the GPU manually before that,
CUDA internally chooses one, which is normally the most powerful GPU
available, but loadindependent.
system = espressomd.System(box_l=[1, 1, 1])
dev = system.cuda_init_handle.device
system.cuda_init_handle.device = dev
The first invocation in the sample above returns the id of the set graphics card, the second one sets the device id.
3.4.1. GPU Acceleration with CUDA¶
Note
Feature CUDA
required
ESPResSo is capable of GPU acceleration to speed up simulations. Not every simulation method is parallelizable or profits from GPU acceleration. Refer to Available simulation methods to check whether your desired method can be used on the GPU. In order to use GPU acceleration you need a NVIDIA GPU and it needs to have at least compute capability 2.0.
For more information please check espressomd.cuda_init.CudaInitHandle
.
3.4.2. List available CUDA devices¶
If you want to list available CUDA devices
you should access espressomd.cuda_init.CudaInitHandle.device_list
, e.g.,
system = espressomd.System(box_l=[1, 1, 1])
print(system.cuda_init_handle.device_list)
This attribute is read only and will return a dictionary containing the device id as key and the device name as its value.
3.4.3. Selection of CUDA device¶
When you start pypresso
your first GPU should be selected.
If you wanted to use the second GPU, this can be done
by setting espressomd.cuda_init.CudaInitHandle.device
as follows:
system = espressomd.System(box_l=[1, 1, 1])
system.cuda_init_handle.device = 1
Setting a device id outside the valid range or a device which does not meet the minimum requirements will raise an exception.