Ferrofluid - Part 1¶
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 figure 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.
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 figure 4) 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}
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
#define DIPOLES
#define DP3M
#define LENNARD_JONES
Please uncomment them in the myconfig.hpp and compile ESPResSo using this myconfig.hpp.
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.
import espressomd
import espressomd.magnetostatics
import espressomd.cluster_analysis
import espressomd.pair_criteria
espressomd.assert_features(['DIPOLES', 'DP3M', 'LENNARD_JONES'])
import numpy as np
import tqdm
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
plt.rcParams.update({'font.size': 18})
Now we set up all simulation parameters.
# 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
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.
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}} $$
BOX_SIZE = (N_PART * np.pi * (LJ_SIGMA / 2.)**2 / PHI)**0.5
# 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
Box size 97.08129562778495
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.
# 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.
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.
# 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))
# 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))))
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.
# Add particles
system.part.add(pos=pos, rotation=N_PART * [(1, 1, 1)], dip=dip, fix=N_PART * [(0, 0, 1)])
<espressomd.particle_data.ParticleSlice at 0x764e5b00880>
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.
# Set integrator to steepest descent method
system.integrator.set_steepest_descent(
f_max=0, gamma=0.1, max_displacement=0.05)
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.
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:.4f}')
energy = energy_new
Minimization, relative change in energy: 1.0000 Minimization, relative change in energy: 0.9983 Minimization, relative change in energy: 1.0000
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.
# Switch to velocity Verlet integrator
system.integrator.set_vv()
system.thermostat.set_langevin(kT=KT, gamma=GAMMA, seed=1)
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.
CI_DP3M_PARAMS = {'cao':3,'r_cut':8.34,'mesh':[8,8,8],'alpha':0.2115,'tune':False}
# Setup dipolar P3M and dipolar layer correction
dp3m = espressomd.magnetostatics.DipolarP3M(accuracy=5E-4, prefactor=DIP_LAMBDA * LJ_SIGMA**3 * KT, **CI_DP3M_PARAMS)
mdlc = espressomd.magnetostatics.DLC(actor=dp3m, maxPWerror=1E-4, gap_size=BOX_SIZE - LJ_SIGMA)
system.actors.add(mdlc)
# 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(f'tuned skin = {system.cell_system.skin:.2f}')
tuned skin = 0.45
Now we equilibrate the dipole-dipole interaction for some time
# equilibrate
EQUIL_ROUNDS = 10
EQUIL_STEPS = 100
for i in tqdm.trange(EQUIL_ROUNDS):
system.integrator.run(EQUIL_STEPS)
100%|██████████████████████████████████████████████████████████████████████| 10/10 [00:03<00:00, 2.53it/s]
Sampling¶
The system will be sampled over 100 loops.
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
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.
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,
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.)
def run(i):
# < excercise >
# Save current system state as a plot
x_data, y_data = # < excercise >
ax.figure.canvas.draw()
part.set_data(x_data, y_data)
print(f'progress: {(i + 1) * 100. / LOOPS:3.0f}%', end='\r')
return part,
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(f'progress: {(i + 1) * 100. / LOOPS:3.0f}%', 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, as well as some isolated monomers.
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)
progress: 100%
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 nearest neighbors are closer than $1.3\sigma_{\text{LJ}}$.
# 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.
n_clusters = []
cluster_sizes = []
Sampling without animation¶
The following code just samples the system and does a cluster analysis every loops (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:
for i in tqdm.trange(LOOPS):
# Run cluster analysis
cluster_structure.run_for_all_pairs()
# Gather statistics:
n_clusters.append(# < excercise >)
for c in cluster_structure.clusters:
cluster_sizes.append(# < excercise >)
system.integrator.run(100)
for i in tqdm.trange(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)
100%|████████████████████████████████████████████████████████████████████| 100/100 [00:42<00:00, 2.34it/s]
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.
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()
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.
size_dist = np.histogram(cluster_sizes, range=(2, 21), bins=19)
Now we can plot this histogram and should see the number of clusters decreasing roughly exponentially with the cluster size, i.e. the number of monomers in the cluster.
import scipy.optimize
def kernel(x, a, b):
return a * np.exp(-np.abs(b) * x)
xdata = size_dist[1][:-2]
ydata = size_dist[0][:-1] / float(LOOPS)
popt, _ = scipy.optimize.curve_fit(kernel, xdata[2:], ydata[2:])
xdata_fit = np.linspace(np.min(xdata), np.max(xdata), 100)
ydata_fit = kernel(xdata_fit, *popt)
plt.figure(figsize=(10, 6))
plt.plot(xdata_fit, ydata_fit, label='Exponential fit')
plt.plot(xdata, ydata, 'o', label='Simulation results')
plt.xlabel('Cluster size', fontsize=20)
plt.ylabel('Probability density function', fontsize=20)
plt.legend()
plt.gca().get_xaxis().set_major_locator(ticker.MaxNLocator(integer=True))
plt.show()
References¶
[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. Chemical Physics Letters 400(1–3): 62–67, 2004. DOI:10.1016/j.cplett.2004.10.086
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
