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 fig. 1). 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 figure 2 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.
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_{cut}$ at the minimum of the potential $r_{cut} = r_{min} = 2^{\frac{1}{6}}\cdot \sigma$ and shifted by $\epsilon_{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_{\text{lj}}\left[ \left( \frac{\sigma}{r_{ij}} \right)^{12} - \left( \frac{\sigma}{r_{ij}} \right)^6 \right] + \varepsilon_{\text{lj}} & r_{ij} < r_{cut} \\ 0 & r_{ij} \geq r_{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}_{DD}}{u_T} = \gamma \frac{\mu^2}{k_BT\sigma^3} \end{equation}where $u_{T} = k_BT$ 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 fig. 4) per particle, i.e. in formulas it reads
\begin{equation} \tilde{u}_{DD} = \frac{ \left| u_{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_B T}B \end{equation}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.
For this tutorial the following features of ESPResSo are needed
#define EXTERNAL_FORCES
#define ROTATION
#define DIPOLES
#define LENNARD_JONES
Please uncomment them in the myconfig.hpp and compile ESPResSo using this myconfig.hpp.
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.
We start with checking for the presence of Espresso features and importing all necessary packages.
In [ ]:
import espressomd
espressomd.assert_features('DIPOLES', 'LENNARD_JONES')
from espressomd.magnetostatics import DipolarP3M
from espressomd.magnetostatic_extensions import DLC
from espressomd.cluster_analysis import ClusterStructure
from espressomd.pair_criteria import DistanceCriterion
import numpy as np
Now we set up all simulation parameters.
In [ ]:
# Lennard-Jones parameters
lj_sigma = 1
lj_epsilon = 1
lj_cut = 2**(1./6.) * lj_sigma
# Particles
N = 1000
# 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
dt = 0.01
Note that we declared a lj_cut. 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. 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.
In [ ]:
# System setup
box_size = (N * np.pi * (lj_sigma/2.)**2 /phi)**0.5
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 = dt
# tune verlet list skin
system.cell_system.tune_skin(min_skin=0.4, max_skin=2., tol=0.2, int_steps=100)
system.thermostat.set_langevin(kT=kT,gamma=gamma,seed=1)
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.
In [ ]:
# Lennard-Jones interaction
system.non_bonded_inter[0,0].lennard_jones.set_params(epsilon=lj_epsilon,sigma=lj_sigma,cutoff=lj_cut, shift="auto")
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.
In [ ]:
# Random dipole moments
np.random.seed(seed = 1)
dip_phi = np.random.random((N,1)) *2. * np.pi
dip_cos_theta = 2*np.random.random((N,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))
# Random positions in the monolayer
pos = box_size* np.hstack((np.random.random((N,2)), np.zeros((N,1))))
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 fix argument.
In [ ]:
# Add particles
system.part.add(pos=pos,rotation=N*[(1,1,1)],dip=dip,fix=N*[(0,0,1)])
Be aware that we do not commit 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 commit the orientation of the dipole moment to the particles and take the magnitude into account when calculating the dipole-dipole interaction with Dipolar P3M if we commit the magnitude to Dipolar P3M together with the prefactor $\gamma$ as a new prefactor
\begin{equation} \tilde{\gamma} = \gamma \mu^2 = \frac{\mu_0}{4\pi}\mu^2 = \lambda \sigma^3 k_BT \end{equation}Of course it would also be possible to commit the whole dipole moment vectors to every particle and let 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. Therefore we integrate until the energy of the whole system is lower than $5 k_BT$.
In [ ]:
# Remove overlap between particles by means of the steepest descent method
system.integrator.set_steepest_descent(
f_max=0,gamma=0.1,max_displacement=0.05)
while system.analysis.energy()["total"] > 5*kT*N:
system.integrator.run(20)
For the simulation of our system we choose the velocity Verlet integrator
In [ ]:
# Switch to velocity Verlet integrator
system.integrator.set_vv()
To calculate the dipole-dipole interaction we use the Dipolar P3M method (see Ref. [1]) 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. [2]). As we add DipolarP3M 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.
In [ ]:
# Setup dipolar P3M and dipolar layer correction
dp3m = DipolarP3M(accuracy=5E-4,prefactor=dip_lambda*lj_sigma**3*kT)
dlc = DLC(maxPWerror=1E-4, gap_size=box_size-lj_sigma)
system.actors.add(dp3m)
system.actors.add(dlc)
# tune verlet list skin again
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))
Now we equilibrate the dipole-dipole interaction for some time
In [ ]:
# Equilibrate
print("Equilibration...")
equil_rounds = 20
equil_steps = 1000
for i in range(equil_rounds):
system.integrator.run(equil_steps)
print("progress: {:3.0f}%, dipolar energy: {:9.2f}".format(
(i + 1) * 100. / equil_rounds, system.analysis.energy()["dipolar"]), end="\r")
print("\nEquilibration done")
As we are interested in the cluster sizes and the number of clusters we now set up the cluster analysis with the distance criterion that a particle is added to a cluster if the nearest neighbors are closer than $1.3\cdot\sigma_{lj}$.
In [ ]:
# Setup cluster analysis
cs=ClusterStructure(pair_criterion=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 sizes.
In [ ]:
n_clusters = []
cluster_sizes = []
We sample over 100 loops
In [ ]:
loops = 100
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
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.
In [ ]:
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from tempfile import NamedTemporaryFile
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 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)
xdata, ydata = [], []
part.set_data(xdata, ydata)
return part,
def run(i):
# Run cluster analysis
cs.run_for_all_pairs()
# Gather statistics:
n_clusters.append(len(cs.clusters))
for c in cs.clusters:
cluster_sizes.append(c[1].size())
system.integrator.run(100)
# Save current system state as a plot
xdata, ydata = system.part[:].pos_folded[:,0], system.part[:].pos_folded[:,1]
ax.figure.canvas.draw()
part.set_data(xdata, ydata)
print("progress: {:3.0f}%".format((i + 1) * 100. / loops), end="\r")
return part,
Now we use the animation class of matplotlib 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 which are possibly connected via Y- and X-links to each other. Some monomers should be also present.
In [ ]:
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)
The following code just samples the system and do a cluster analysis every loops (100 by default) simulation steps.
In [ ]:
for i in range(loops):
# Run cluster analysis
cs.run_for_all_pairs()
# Gather statistics:
n_clusters.append(len(cs.clusters))
for c in cs.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 matplotlib.
In [ ]:
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[:].pos_folded[:,0], system.part[:].pos_folded[:,1], 'o')
plt.show()
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.
After sampling our system for a while we now can calculate estimators for the expectation value of the cluster sizes and their distribution. Thus we calculate a histogram over all cluster sizes using the numpy function histogram.
Note: We want to get a histogram up to a cluster size of 19 particles. To not count clusters with size of 20 in our last bin (our bin for cluster size 19) we add an additional bin which we will not use in the following. This has to do with the way numpy defines the intervals between two bin edges. All but the last bin is half-open with the open border at the higher bin edge.
In [ ]:
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
In [ ]:
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()
[1] Juan J. Cerdà, V. Ballenegger, O. Lenz, and Ch. Holm. P3M algorithm for dipolar interactions. Journal of Chemical Physics, 129:234104, 2008.
[2] A. Bródka. “Ewald summation method with electrostatic layer correction for
interactions of point dipoles in slab geometry”. In: Chemical Physics Letters 400.1
(2004), pp. 62–67. issn: 0009-2614. doi: https://doi.org/10.1016/j.cplett.
2004.10.086. url: http://www.sciencedirect.com/science/article/pii/S0009261404016847.
Image sources:
[3] Ayouril, Electro-Steric Stabilization, CC BY-SA 3.0
[4] Gregory F. Maxwell <gmaxwell@gmail.com> PGP:0xB0413BFA, Ferrofluid Magnet under glass edit, CC BY-SA 3.0
