<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Ferrofluid - Part 1
Table of Contents
Introduction
The Model
Structure of this tutorial
Compiling ESPResSo for this Tutorial
A Monolayer-Ferrofluid System in ESPResSo
Setup
Sampling
Sampling with animation
Sampling without animation
Cluster distribution
Introduction
Ferrofluids are colloidal suspensions of ferromagnetic single-domain particles in a liquid carrier. As the single particles contain only one magnetic domain, they can be seen as small permanent magnets. To prevent agglomeration of the particles, due to van-der-Waals or magnetic attraction, they are usually sterically or electrostatically stabilized (see <a href='#fig_1'>figure 1</a>). The former is achieved by adsorption of long chain molecules onto the particle surface, the latter by adsorption of charged coating particles. The size of the ferromagnetic particles are in the region of 10 nm. With the surfactant layer added they can reach a size of a few hundred nanometers. Have in mind that if we refer to the particle diameter $\sigma$ we mean the diameter of the magnetic core plus two times the thickness of the surfactant layer.
Some of the liquid properties, like the viscosity, the phase behavior or the optical birefringence can be altered via an external magnetic field or simply the fluid can be guided by such an
field. Thus ferrofluids possess a wide range of biomedical applications like  magnetic drug
targeting or magnetic thermoablation and technical applications like fine positioning systems or adaptive bearings and
dampers.
In <a href='#fig_2'>figure 2</a> the picture of a ferrofluid exposed to the magnetic field of a permanent magnet is shown. The famous energy minimizing thorn-like surface is clearly visible.
<a id='fig_1'></a><figure>
<img src="figures/Electro-Steric_Stabilization.jpg" style="float
Step1:  Now we set up all simulation parameters.
Step2:  Note that we declared a <tt>lj_cut</tt>. This will be used as the cut-off radius of the Lennard-Jones potential to obtain a purely repulsive WCA potential.
Now we set up the system. The length of the simulation box is calculated using the desired area fraction and the area all particles occupy. Then we create the ESPResSo system and pass the simulation step. For the Verlet list skin parameter we use the built-in tuning algorithm of ESPResSo.
Exercise
Step3:  Now we set up the interaction between the particles as a non-bonded interaction and use the Lennard-Jones potential as the interaction potential. Here we use the above mentioned cut-off radius to get a purely repulsive interaction.
Step4:  Now we generate random positions and orientations of the particles and their dipole moments. 
Hint
Step5:  Now we add the particles with their positions and orientations to our system. Thereby we activate all degrees of freedom for the orientation of the dipole moments. As we want a two dimensional system we only allow the particles to translate in $x$- and $y$-direction and not in $z$-direction by using the <tt>fix</tt> argument.
Step6:  Be aware that we do not set the magnitude of the magnetic dipole moments to the particles. As in our case all particles have the same dipole moment it is possible to rewrite the dipole-dipole interaction potential to
\begin{equation}
    u_{\text{DD}}(\vec{r}{ij}, \vec{\mu}_i, \vec{\mu}_j) =  \gamma \mu^2 \left(\frac{\vec{\hat{\mu}}_i \cdot \vec{\hat{\mu}}_j}{r{ij}^3} - 3\frac{(\vec{\hat{\mu}}i \cdot \vec{r}{ij}) \cdot (\vec{\hat{\mu}}j \cdot \vec{r}{ij})}{r_{ij}^5}\right)
\end{equation}
where $\vec{\hat{\mu}}_i$ is the unit vector of the dipole moment $i$ and $\mu$ is the magnitude of the dipole moments.
Thus we can only prescribe the initial orientation of the dipole moment to the particles and take the magnitude of the moments into account when calculating the dipole-dipole interaction with Dipolar P3M, by modifying the original Dipolar P3M prefactor $\gamma$ such that
\begin{equation}
    \tilde{\gamma} = \gamma \mu^2 = \frac{\mu_0}{4\pi}\mu^2 = \lambda \sigma^3 k_{\text{B}}T
\end{equation}
Of course it would also be possible to prescribe the whole dipole moment vectors to every particle and leave the prefactor of Dipolar P3M unchanged ($\gamma$). In fact we have to do this if we want to simulate polydisperse systems.
Now we choose the steepest descent integrator to remove possible overlaps of the particles.
Step7:  Exercise
Step8:  To calculate the dipole-dipole interaction we use the Dipolar P3M method (see Ref. <a href='#[1]'>[1]</a>) which is based on the Ewald summation. By default the boundary conditions of the system are set to conducting which means the dielectric constant is set to infinity for the surrounding medium. As we want to simulate a two dimensional system we additionally use the dipolar layer correction (DLC) (see Ref. <a href='#[2]'>[2]</a>). As we add <tt>DipolarP3M</tt> to our system as an actor, a tuning function is started automatically which tries to find the optimal parameters for Dipolar P3M and prints them to the screen. The last line of the output is the value of the tuned skin.
Step9:  Now we equilibrate the dipole-dipole interaction for some time
Step10:  Sampling
The system will be sampled over 100 loops.
Step12:  As the system is two dimensional, we can simply do a scatter plot to get a visual representation of a system state. To get a better insight of how a ferrofluid system develops during time we will create a video of the development of our system during the sampling. If you only want to sample the system simply go to Sampling without animation
Sampling with animation
To get an animation of the system development we have to create a function which will save the video and embed it in an html string.
Step13:  Exercise
Step14:  Cluster analysis
To quantify the number of clusters and their respective sizes, we now want to perform a cluster analysis.
For that we can use ESPREsSo's cluster analysis class.
Exercise
Step15:  Sampling without animation
The following code just samples the system and does a cluster analysis every <tt>loops</tt> (100 by default) simulation steps.
Exercise
Step16:  In the plot chain-like and ring-like clusters should be visible. Some of them are connected via Y- or X-links to each other. Also some monomers should be present.
Cluster distribution
After having sampled our system we now can calculate estimators for the expectation value of the cluster sizes and their distribution.
Exercise
<ASSISTANT_TASK:>
Python Code:
import espressomd
import espressomd.magnetostatics
import espressomd.magnetostatic_extensions
import espressomd.cluster_analysis
import espressomd.pair_criteria
espressomd.assert_features('DIPOLES', 'LENNARD_JONES')
import numpy as np
Explanation: Ferrofluid - Part 1
Table of Contents
Introduction
The Model
Structure of this tutorial
Compiling ESPResSo for this Tutorial
A Monolayer-Ferrofluid System in ESPResSo
Setup
Sampling
Sampling with animation
Sampling without animation
Cluster distribution
Introduction
Ferrofluids are colloidal suspensions of ferromagnetic single-domain particles in a liquid carrier. As the single particles contain only one magnetic domain, they can be seen as small permanent magnets. To prevent agglomeration of the particles, due to van-der-Waals or magnetic attraction, they are usually sterically or electrostatically stabilized (see <a href='#fig_1'>figure 1</a>). The former is achieved by adsorption of long chain molecules onto the particle surface, the latter by adsorption of charged coating particles. The size of the ferromagnetic particles are in the region of 10 nm. With the surfactant layer added they can reach a size of a few hundred nanometers. Have in mind that if we refer to the particle diameter $\sigma$ we mean the diameter of the magnetic core plus two times the thickness of the surfactant layer.
Some of the liquid properties, like the viscosity, the phase behavior or the optical birefringence can be altered via an external magnetic field or simply the fluid can be guided by such an
field. Thus ferrofluids possess a wide range of biomedical applications like  magnetic drug
targeting or magnetic thermoablation and technical applications like fine positioning systems or adaptive bearings and
dampers.
In <a href='#fig_2'>figure 2</a> the picture of a ferrofluid exposed to the magnetic field of a permanent magnet is shown. The famous energy minimizing thorn-like surface is clearly visible.
<a id='fig_1'></a><figure>
<img src="figures/Electro-Steric_Stabilization.jpg" style="float: center; width: 49%">
<center>
<figcaption>Figure 1: Schematic representation of electrostatically stabilization (picture top) and steric stabilization (picture bottom) <a href='#[3]'>[3]</a></figcaption>
</center>
</figure>
<a id='fig_2'></a><figure>
<img src='figures/Ferrofluid_Magnet_under_glass_edit.jpg' alt='ferrofluid on glass plate under which a strong magnet is placed' style='width: 600px;'/>
<center>
<figcaption>Figure 2: Real Ferrofluid exposed to an external magnetic field (neodymium magnet) <a href='#[4]'>[4]</a></figcaption>
</center>
</figure>
The Model
For simplicity in this tutorial we simulate spherical particles in a monodisperse ferrofluid system which means all particles have the same diameter $\sigma$ and dipole moment $\mu$. The point dipole moment is placed at the center of the particles and is constant both in magnitude and direction (in the coordinate system of the particle). This can be justified as the Néel relaxation times are usually negligible for the usual sizes of ferrofluid particles.
Thus the magnetic interaction potential between two single particles is the dipole-dipole interaction potential which reads
\begin{equation}
        u_{\text{DD}}(\vec{r}{ij}, \vec{\mu}_i, \vec{\mu}_j) = \gamma \left(\frac{\vec{\mu}_i \cdot \vec{\mu}_j}{r{ij}^3} - 3\frac{(\vec{\mu}i \cdot \vec{r}{ij}) \cdot (\vec{\mu}j \cdot \vec{r}{ij})}{r_{ij}^5}\right)
\end{equation}
with $\gamma = \frac{\mu_0}{4 \pi}$ and $\mu_0$ the vacuum permeability.
For the steric interaction in this tutorial we use the purely repulsive Weeks-Chandler-Andersen (WCA) potential which is a Lennard-Jones potential with cut-off radius $r_{\text{cut}}$ at the minimum of the potential $r_{\text{cut}} = r_{\text{min}} = 2^{\frac{1}{6}}\cdot \sigma$ and shifted by $\varepsilon_{ij}$ such that the potential is continuous at the cut-off radius. Thus the potential has the shape
\begin{equation}
        u_{\text{sr}}^{\text{WCA}}(r_{ij}) = \left{
        \begin{array}{ll}
        4\varepsilon_{ij}\left[ \left( \frac{\sigma}{r_{ij}} \right)^{12} - \left( \frac{\sigma}{r_{ij}} \right)^6 \right] + \varepsilon_{ij} & r_{ij} < r_{\text{cut}} \
        0 & r_{ij} \geq r_{\text{cut}} \
        \end{array}
        \right.
\end{equation}
where $r_{ij}$ are the distances between two particles.
The purely repulsive character of the potential can be justified by the fact that the particles in real ferrofluids are sterically or electrostatically stabilized against agglomeration.
The whole interaction potential reads
\begin{equation}
        u(\vec{r}{ij}, \vec{\mu}_i, \vec{\mu}_j) = u{\text{sr}}(\vec{r}{ij}) + u{\text{DD}}(\vec{r}_{ij}, \vec{\mu}_i, \vec{\mu}_j)
\end{equation}
The liquid carrier of the system is simulated through a Langevin thermostat.
For ferrofluid systems there are three important parameters. The first is the volume fraction in three dimensions or the area fraction in two dimensions or quasi two dimensions. The second is the dipolar interaction parameter $\lambda$
\begin{equation}
    \lambda = \frac{\tilde{u}{\text{DD}}}{u_T} = \gamma \frac{\mu^2}{k{\text{B}}T\sigma^3}
\end{equation}
where $u_\mathrm{T} = k_{\text{B}}T$ is the thermal energy and $\tilde{u}_{DD}$ is the absolute value of the dipole-dipole interaction energy at close contact (cc) and head-to-tail configuration (htt) (see <a href='#fig_4'>figure 4</a>) per particle, i.e. in formulas it reads
\begin{equation}
    \tilde{u}{\text{DD}} = \frac{ \left| u{\text{DD}}^{\text{htt, cc}} \right| }{2}
\end{equation}
The third parameter takes a possible external magnetic field into account and is called Langevin parameter $\alpha$. It is the ratio between the energy of a dipole moment in the external magnetic field $B$ and the thermal energy
\begin{equation}
    \alpha = \frac{\mu_0 \mu}{k_{\text{B}} T}B
\end{equation}
<a id='fig_4'></a><figure>
<img src='figures/headtotailconf.png' alt='schematic representation of head to tail configuration' style='width: 200px;'/>
<center>
<figcaption>Figure 4: Schematic representation of the head-to-tail configuration of two magnetic particles at close contact.</figcaption>
</center>
</figure>
Structure of this tutorial
The aim of this tutorial is to introduce the basic features of ESPResSo for ferrofluids or dipolar fluids in general. In part I and part II we will do this for a monolayer-ferrofluid, in part III for a three dimensional system. In part I we will examine the clusters which are present in all interesting ferrofluid systems. In part II we will examine the influence of the dipole-dipole-interaction on the magnetization curve of a ferrofluid. In part III we calculate estimators for the initial susceptibility using fluctuation formulas and sample the magnetization curve.
We assume the reader is familiar with the basic concepts of Python and MD simulations.
Remark: The equilibration and sampling times used in this tutorial would be not sufficient for scientific  purposes, but they are long enough to get at least a qualitative insight of the behaviour of ferrofluids. They have been shortened so we achieve reasonable computation times for the purpose of a tutorial.
Compiling ESPResSo for this Tutorial
For this tutorial the following features of ESPResSo are needed
```c++
define EXTERNAL_FORCES
define ROTATION
define DIPOLES
define LENNARD_JONES
```
Please uncomment them in the <tt>myconfig.hpp</tt> and compile ESPResSo using this <tt>myconfig.hpp</tt>.
A Monolayer-Ferrofluid System in ESPResSo
For interesting ferrofluid systems, where the fraction of ferromagnetic particles in the liquid carrier and their dipole moment are not vanishingly small, the ferromagnetic particles form clusters of different shapes and sizes. If the fraction and/or dipole moments are big enough the clusters can interconnect with each other and form a whole space occupying network.
In this part we want to investigate the number of clusters as well as their shape and size in our simulated monolayer ferrofluid system. It should be noted that a monolayer is a quasi three dimensional system (q2D), i.e. two dimensional for the positions and three dimensional for the orientation of the dipole moments.
Setup
We start with checking for the presence of ESPResSo features and importing all necessary packages.
End of explanation
# Lennard-Jones parameters
LJ_SIGMA = 1
LJ_EPSILON = 1
LJ_CUT = 2**(1. / 6.) * LJ_SIGMA
# Particles
N_PART = 1200
# Area fraction of the mono-layer
PHI = 0.1
# Dipolar interaction parameter lambda = mu_0 m^2 /(4 pi sigma^3 KT)
DIP_LAMBDA = 4
# Temperature
KT = 1.0
# Friction coefficient
GAMMA = 1.0
# Time step
TIME_STEP = 0.01
Explanation: Now we set up all simulation parameters.
End of explanation
# System setup
# BOX_SIZE = ...
print("Box size", BOX_SIZE)
# Note that the dipolar P3M and dipolar layer correction need a cubic
# simulation box for technical reasons.
system = espressomd.System(box_l=(BOX_SIZE, BOX_SIZE, BOX_SIZE))
system.time_step = TIME_STEP
Explanation: Note that we declared a <tt>lj_cut</tt>. This will be used as the cut-off radius of the Lennard-Jones potential to obtain a purely repulsive WCA potential.
Now we set up the system. The length of the simulation box is calculated using the desired area fraction and the area all particles occupy. Then we create the ESPResSo system and pass the simulation step. For the Verlet list skin parameter we use the built-in tuning algorithm of ESPResSo.
Exercise:
How large does BOX_SIZE have to be for a system of N_PART particles with a volume (area) fraction PHI?
Define BOX_SIZE.
$$
L_{\text{box}} = \sqrt{\frac{N A_{\text{sphere}}}{\varphi}}
$$
python
BOX_SIZE = (N_PART * np.pi * (LJ_SIGMA / 2.)**2 / PHI)**0.5
End of explanation
# Lennard-Jones interaction
system.non_bonded_inter[0, 0].lennard_jones.set_params(epsilon=LJ_EPSILON, sigma=LJ_SIGMA, cutoff=LJ_CUT, shift="auto")
Explanation: Now we set up the interaction between the particles as a non-bonded interaction and use the Lennard-Jones potential as the interaction potential. Here we use the above mentioned cut-off radius to get a purely repulsive interaction.
End of explanation
# Random dipole moments
# ...
# dip = ...
# Random positions in the monolayer
pos = BOX_SIZE * np.hstack((np.random.random((N_PART, 2)), np.zeros((N_PART, 1))))
Explanation: Now we generate random positions and orientations of the particles and their dipole moments. 
Hint:
It should be noted that we seed the random number generator of numpy. Thus the initial configuration of our system is the same every time this script will be executed. You can change it to another one to simulate with a different initial configuration.
Exercise:
How does one set up randomly oriented dipole moments?
Hint: Think of the way that different methods could introduce a bias in the distribution of the orientations.
Create a variable dip as a N_PART x 3 numpy array, which contains the randomly distributed dipole moments.
```python
Random dipole moments
np.random.seed(seed=1)
dip_phi = 2. * np.pi * np.random.random((N_PART, 1))
dip_cos_theta = 2 * np.random.random((N_PART, 1)) - 1
dip_sin_theta = np.sin(np.arccos(dip_cos_theta))
dip = np.hstack((
    dip_sin_theta * np.sin(dip_phi),
    dip_sin_theta * np.cos(dip_phi),
    dip_cos_theta))
```
End of explanation
# Add particles
system.part.add(pos=pos, rotation=N_PART * [(1, 1, 1)], dip=dip, fix=N_PART * [(0, 0, 1)])
Explanation: Now we add the particles with their positions and orientations to our system. Thereby we activate all degrees of freedom for the orientation of the dipole moments. As we want a two dimensional system we only allow the particles to translate in $x$- and $y$-direction and not in $z$-direction by using the <tt>fix</tt> argument.
End of explanation
# Set integrator to steepest descent method
system.integrator.set_steepest_descent(
    f_max=0, gamma=0.1, max_displacement=0.05)
Explanation: Be aware that we do not set the magnitude of the magnetic dipole moments to the particles. As in our case all particles have the same dipole moment it is possible to rewrite the dipole-dipole interaction potential to
\begin{equation}
    u_{\text{DD}}(\vec{r}{ij}, \vec{\mu}_i, \vec{\mu}_j) =  \gamma \mu^2 \left(\frac{\vec{\hat{\mu}}_i \cdot \vec{\hat{\mu}}_j}{r{ij}^3} - 3\frac{(\vec{\hat{\mu}}i \cdot \vec{r}{ij}) \cdot (\vec{\hat{\mu}}j \cdot \vec{r}{ij})}{r_{ij}^5}\right)
\end{equation}
where $\vec{\hat{\mu}}_i$ is the unit vector of the dipole moment $i$ and $\mu$ is the magnitude of the dipole moments.
Thus we can only prescribe the initial orientation of the dipole moment to the particles and take the magnitude of the moments into account when calculating the dipole-dipole interaction with Dipolar P3M, by modifying the original Dipolar P3M prefactor $\gamma$ such that
\begin{equation}
    \tilde{\gamma} = \gamma \mu^2 = \frac{\mu_0}{4\pi}\mu^2 = \lambda \sigma^3 k_{\text{B}}T
\end{equation}
Of course it would also be possible to prescribe the whole dipole moment vectors to every particle and leave the prefactor of Dipolar P3M unchanged ($\gamma$). In fact we have to do this if we want to simulate polydisperse systems.
Now we choose the steepest descent integrator to remove possible overlaps of the particles.
End of explanation
# Switch to velocity Verlet integrator
system.integrator.set_vv()
system.thermostat.set_langevin(kT=KT, gamma=GAMMA, seed=1)
Explanation: Exercise:
Perform a steepest descent energy minimization.
Track the relative energy change $E_{\text{rel}}$ per minimization loop (where the integrator is run for 10 steps) and terminate once $E_{\text{rel}} \le 0.05$, i.e. when there is less than a 5% difference in the relative energy change in between iterations.
```python
import sys
energy = system.analysis.energy()['total']
relative_energy_change = 1.0
while relative_energy_change > 0.05:
    system.integrator.run(10)
    energy_new = system.analysis.energy()['total']
    # Prevent division by zero errors:
    if energy < sys.float_info.epsilon:
        break
    relative_energy_change = (energy - energy_new) / energy
    print(f"Minimization, relative change in energy: {relative_energy_change}")
    energy = energy_new
```
For the simulation of our system we choose the velocity Verlet integrator.
After that we set up the thermostat which is, in our case, a Langevin thermostat to simulate in an NVT ensemble.
Hint:
It should be noted that we seed the Langevin thermostat, thus the time evolution of the system is partly predefined. Partly because of the numeric accuracy and the automatic tuning algorithms of Dipolar P3M and DLC where the resulting parameters are slightly different every time. You can change the seed to get a guaranteed different time evolution.
End of explanation
CI_DP3M_PARAMS = {} # debug variable for continuous integration, can be left empty
# Setup dipolar P3M and dipolar layer correction
dp3m = espressomd.magnetostatics.DipolarP3M(accuracy=5E-4, prefactor=DIP_LAMBDA * LJ_SIGMA**3 * KT, **CI_DP3M_PARAMS)
dlc = espressomd.magnetostatic_extensions.DLC(maxPWerror=1E-4, gap_size=BOX_SIZE - LJ_SIGMA)
system.actors.add(dp3m)
system.actors.add(dlc)
# tune verlet list skin
system.cell_system.tune_skin(min_skin=0.4, max_skin=2., tol=0.2, int_steps=100)
# print skin value
print('tuned skin = {}'.format(system.cell_system.skin))
Explanation: To calculate the dipole-dipole interaction we use the Dipolar P3M method (see Ref. <a href='#[1]'>[1]</a>) which is based on the Ewald summation. By default the boundary conditions of the system are set to conducting which means the dielectric constant is set to infinity for the surrounding medium. As we want to simulate a two dimensional system we additionally use the dipolar layer correction (DLC) (see Ref. <a href='#[2]'>[2]</a>). As we add <tt>DipolarP3M</tt> to our system as an actor, a tuning function is started automatically which tries to find the optimal parameters for Dipolar P3M and prints them to the screen. The last line of the output is the value of the tuned skin.
End of explanation
# Equilibrate
print("Equilibration...")
EQUIL_ROUNDS = 20
EQUIL_STEPS = 1000
for i in range(EQUIL_ROUNDS):
    system.integrator.run(EQUIL_STEPS)
    print(
        f"progress: {(i + 1) * 100. / EQUIL_ROUNDS}%, dipolar energy: {system.analysis.energy()['dipolar']}",
        end="\r")
print("\nEquilibration done")
Explanation: Now we equilibrate the dipole-dipole interaction for some time
End of explanation
LOOPS = 100
Explanation: Sampling
The system will be sampled over 100 loops.
End of explanation
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import tempfile
import base64
VIDEO_TAG = <video controls>
 <source src="data:video/x-m4v;base64,{0}" type="video/mp4">
 Your browser does not support the video tag.
</video>
def anim_to_html(anim):
    if not hasattr(anim, '_encoded_video'):
        with tempfile.NamedTemporaryFile(suffix='.mp4') as f:
            anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264'])
            with open(f.name, "rb") as g:
                video = g.read()
        anim._encoded_video = base64.b64encode(video).decode('ascii')
        plt.close(anim._fig)
    return VIDEO_TAG.format(anim._encoded_video)
animation.Animation._repr_html_ = anim_to_html
def init():
    # Set x and y range
    ax.set_ylim(0, BOX_SIZE)
    ax.set_xlim(0, BOX_SIZE)
    x_data, y_data = [], []
    part.set_data(x_data, y_data)
    return part,
Explanation: As the system is two dimensional, we can simply do a scatter plot to get a visual representation of a system state. To get a better insight of how a ferrofluid system develops during time we will create a video of the development of our system during the sampling. If you only want to sample the system simply go to Sampling without animation
Sampling with animation
To get an animation of the system development we have to create a function which will save the video and embed it in an html string.
End of explanation
fig, ax = plt.subplots(figsize=(10, 10))
part, = ax.plot([], [], 'o')
animation.FuncAnimation(fig, run, frames=LOOPS, blit=True, interval=0, repeat=False, init_func=init)
Explanation: Exercise:
In the following an animation loop is defined, however it is incomplete.
Extend the code such that in every loop the system is integrated for 100 steps.
Afterwards x_data and y_data have to be populated by the folded $x$- and $y$- positions of the particles.
(You may copy and paste the incomplete code template to the empty cell below.)
```python
def run(i):
    # < excercise >
# Save current system state as a plot
x_data, y_data = # &lt; excercise &gt;
ax.figure.canvas.draw()
part.set_data(x_data, y_data)
print("progress: {:3.0f}%".format((i + 1) * 100. / LOOPS), end="\r")
return part,
```
```python
def run(i):
    system.integrator.run(100)
# Save current system state as a plot
x_data, y_data = system.part.all().pos_folded[:, 0], system.part.all().pos_folded[:, 1]
ax.figure.canvas.draw()
part.set_data(x_data, y_data)
print("progress: {:3.0f}%".format((i + 1) * 100. / LOOPS), end="\r")
return part,
```
Now we use the <tt>animation</tt> class of <tt>matplotlib</tt> to save snapshots of the system as frames of a video which is then displayed after the sampling is finished. Between two frames are 100 integration steps.
In the video chain-like and ring-like clusters should be visible, as well as some isolated monomers.
End of explanation
n_clusters = []
cluster_sizes = []
Explanation: Cluster analysis
To quantify the number of clusters and their respective sizes, we now want to perform a cluster analysis.
For that we can use ESPREsSo's cluster analysis class.
Exercise:
Setup a cluster analysis object (ClusterStructure class) and assign its instance to the variable cluster_structure.
As criterion for the cluster analysis use a distance criterion where particles are assumed to be part of a cluster if the neaest neighbors are closer than $1.3\sigma_{\text{LJ}}$.
```python
Setup cluster analysis
cluster_structure = espressomd.cluster_analysis.ClusterStructure(pair_criterion=espressomd.pair_criteria.DistanceCriterion(cut_off=1.3 * LJ_SIGMA))
```
Now we sample our system for some time and do a cluster analysis in order to get an estimator of the cluster observables.
For the cluster analysis we create two empty lists. The first for the number of clusters and the second for their respective sizes.
End of explanation
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 10))
plt.xlim(0, BOX_SIZE)
plt.ylim(0, BOX_SIZE)
plt.xlabel('x-position', fontsize=20)
plt.ylabel('y-position', fontsize=20)
plt.plot(system.part.all().pos_folded[:, 0], system.part.all().pos_folded[:, 1], 'o')
plt.show()
Explanation: Sampling without animation
The following code just samples the system and does a cluster analysis every <tt>loops</tt> (100 by default) simulation steps.
Exercise:
Write an integration loop which runs a cluster analysis on the system, saving the number of clusters n_clusters and the size distribution cluster_sizes.
Take the following as a starting point:
```python
for i in range(LOOPS):
    # Run cluster analysis
    cluster_structure.run_for_all_pairs()
# Gather statistics:
n_clusters.append(# &lt; excercise &gt;)
for c in cluster_structure.clusters:
    cluster_sizes.append(# &lt; excercise &gt;)
system.integrator.run(100)
print("progress: {:3.0f}%".format((float(i)+1)/LOOPS * 100), end="\r")
```
```python
for i in range(LOOPS):
    # Run cluster analysis
    cluster_structure.run_for_all_pairs()
# Gather statistics:
n_clusters.append(len(cluster_structure.clusters))
for c in cluster_structure.clusters:
    cluster_sizes.append(c[1].size())
system.integrator.run(100)
print("progress: {:3.0f}%".format((float(i) + 1) / LOOPS * 100), end="\r")
```
You may want to get a visualization of the current state of the system. For that we plot the particle positions folded to the simulation box using <tt>matplotlib</tt>.
End of explanation
plt.figure(figsize=(10, 10))
plt.grid()
plt.xticks(range(0, 20))
plt.plot(size_dist[1][:-2], size_dist[0][:-1] / float(LOOPS))
plt.xlabel('size of clusters', fontsize=20)
plt.ylabel('distribution', fontsize=20)
plt.show()
Explanation: In the plot chain-like and ring-like clusters should be visible. Some of them are connected via Y- or X-links to each other. Also some monomers should be present.
Cluster distribution
After having sampled our system we now can calculate estimators for the expectation value of the cluster sizes and their distribution.
Exercise:
Use numpy to calculate a histogram of the cluster sizes and assign it to the variable size_dist.
Take only clusters up to a size of 19 particles into account.
Hint: In order not to count clusters with size 20 or more, one may include an additional bin containing these.
The reason for that is that numpy defines the histogram bins as half-open intervals with the open border at the upper bin edge.
Consequently clusters of larger sizes are attributed to the last bin.
By not using the last bin in the plot below, these clusters can effectively be neglected.
python
size_dist = np.histogram(cluster_sizes, range=(2, 21), bins=19)
Now we can plot this histogram and should see an exponential decrease in the number of particles in a cluster along the size of a cluster, i.e. the number of monomers in it
End of explanation
<END_TASK>