Ferrofluid - Part II

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.

Applying an external magnetic field

In this part we want to investigate the influence of a homogeneous external magnetic field exposed to a ferrofluid system.

We import all necessary packages and check for the required ESPResSo features

In [1]:
import espressomd
espressomd.assert_features('DIPOLES', 'LENNARD_JONES')

from espressomd.magnetostatics import DipolarP3M
from espressomd.magnetostatic_extensions import DLC

import numpy as np

and set up the simulation parameters where we introduce a new dimensionless parameter

\begin{equation} \alpha = \frac{\mu B}{k_BT} = \frac{\mu \mu_0 H}{k_BT} \end{equation}

which is called Langevin parameter. We intentionally choose a relatively high volume fraction $\phi$ and dipolar interaction parameter $\lambda$ to clearly see the influence of the dipole-dipole interaction

In [2]:
# Lennard-Jones parameters
lj_sigma = 1.
lj_epsilon = 1.
lj_cut = 2**(1./6.) * lj_sigma

# Particles
N = 700

# Area fraction of the mono-layer 
phi = 0.06

# 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

# Langevin parameter alpha = mu_0 m H / kT
alpha = 10.

# vacuum permeability
mu_0 = 1.

Now we set up the system, where we, as we did in part I, only commit the orientation of the dipole moment to the particles and take the magnitude into account in the prefactor of Dipolar P3M (for more details see part I).

Hint: It should be noted that we seed both the Langevin thermostat and the random number generator of numpy. The latter means that the initial configuration of our system is the same every time this script will be executed. As the time evolution of the system depends not solely on the Langevin thermostat but also on the numeric accuracy and DP3M as well as DLC (the tuned parameters are slightly different every time) it is only partly predefined. You can change the seeds to simulate with a different initial configuration and a guaranteed different time evolution.

In [3]:
# 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
system.thermostat.set_langevin(kT=kT,gamma=gamma, seed=1)

# Lennard-Jones interaction
system.non_bonded_inter[0,0].lennard_jones.set_params(epsilon=lj_epsilon,sigma=lj_sigma,cutoff=lj_cut, shift="auto")

# 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))))

# Add particles
system.part.add(pos=pos,rotation=N*[(1,1,1)],dip=dip,fix=N*[(0,0,1)])

# 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)

# Switch to velocity Verlet integrator
system.integrator.set_vv()

# tune verlet list skin
system.cell_system.tune_skin(min_skin=0.4, max_skin=2., tol=0.2, int_steps=100)

# Setup dipolar P3M and dipolar layer correction (DLC)
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))
Box size 95.72344839677595
Dipolar P3M tune parameters: Accuracy goal = 5.00000e-04 prefactor = 4.00000e+00
System: box_l = 9.57234e+01 # charged part = 700 Sum[q_i^2] = 7.00000e+02
Dmesh cao Dr_cut_iL   Dalpha_L     Derr         Drs_err    Dks_err    time [ms]
8    3   7.31480e-02 2.54050e+01 4.96625e-04 3.536e-04 3.488e-04 1       
8    2   7.97126e-02 2.19599e+01 4.95928e-04 3.536e-04 3.478e-04 1       
8    4   6.93968e-02 2.76533e+01 4.99061e-04 3.536e-04 3.522e-04 1       
8    1   9.19039e-02 1.67627e+01 4.94683e-04 3.536e-04 3.460e-04 1       
10   2   7.97126e-02 2.19599e+01 5.85508e-04 3.536e-04 4.667e-04 accuracy not achieved
10   3   7.97126e-02 2.19599e+01 5.31120e-04 3.536e-04 3.963e-04 accuracy not achieved
10   4   7.97126e-02 2.19599e+01 5.05715e-04 3.536e-04 3.616e-04 accuracy not achieved
10   5   7.90898e-02 2.22637e+01 4.98738e-04 3.536e-04 3.518e-04 2       
10   6   7.84671e-02 2.25719e+01 4.98056e-04 3.536e-04 3.508e-04 3       
10   7   7.84671e-02 2.25719e+01 4.95362e-04 3.536e-04 3.470e-04 5       
12   2   7.90898e-02 2.22637e+01 6.65028e-04 3.536e-04 5.633e-04 accuracy not achieved
12   3   7.90898e-02 2.22637e+01 5.83887e-04 3.536e-04 4.647e-04 accuracy not achieved
12   4   7.90898e-02 2.22637e+01 5.43132e-04 3.536e-04 4.123e-04 accuracy not achieved
12   5   7.90898e-02 2.22637e+01 5.20953e-04 3.536e-04 3.826e-04 accuracy not achieved
12   6   7.90898e-02 2.22637e+01 5.08993e-04 3.536e-04 3.662e-04 accuracy not achieved
12   7   7.90898e-02 2.22637e+01 5.02675e-04 3.536e-04 3.573e-04 accuracy not achieved
14   2   7.90898e-02 2.22637e+01 7.01387e-04 3.536e-04 6.058e-04 accuracy not achieved
14   3   7.90898e-02 2.22637e+01 6.00259e-04 3.536e-04 4.851e-04 accuracy not achieved
14   4   7.90898e-02 2.22637e+01 5.49254e-04 3.536e-04 4.203e-04 accuracy not achieved
14   5   7.90898e-02 2.22637e+01 5.19944e-04 3.536e-04 3.812e-04 accuracy not achieved
14   6   7.90898e-02 2.22637e+01 5.02711e-04 3.536e-04 3.574e-04 accuracy not achieved
14   7   7.90898e-02 2.22637e+01 4.92655e-04 3.536e-04 3.431e-04 5       
14   6   7.90898e-02 2.22637e+01 5.02711e-04 3.536e-04 3.574e-04 accuracy not achieved

resulting parameters: mesh: 8, cao: 2, r_cut_iL: 7.9713e-02,
                      alpha_L: 2.1960e+01, accuracy: 4.9593e-04, time: 1

tuned skin = 1.1

We now apply the external magnetic field which is

\begin{equation} B = \mu_0 H = \frac{\alpha~k_BT}{\mu} \end{equation}

As only the current orientation of the dipole moments, i.e. the unit vector of the dipole moments, is saved in the particle list but not their strength we have to commit $B\cdot \mu$ as the external magnetic field to ESPResSo. We applying the field in x-direction using the class constraints of ESPResSo

In [4]:
# magnetic field times dipole moment
H_dipm = (alpha*kT)
H_field = [H_dipm,0,0]
print("Set magnetic field constraint...")
H_constraint = espressomd.constraints.HomogeneousMagneticField(H=H_field)
system.constraints.add(H_constraint)
Set magnetic field constraint...
Out[4]:
<espressomd.constraints.HomogeneousMagneticField at 0x7f310766ea98>

and equilibrate the system for a while

In [5]:
# Equilibrate
print("Equilibration...")
equil_rounds = 10
equil_steps = 100
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")
Equilibration...
progress: 100%, dipolar energy:  -9647.34
Equilibration done

Now we can visualize the current state and see that the particles mostly create chains oriented in the direction of the external magnetic field. Also some monomers should be present.

In [6]:
import matplotlib.pyplot as plt
In [7]:
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()

Video of the development of the system

You may want to get an insight of how the system develops in time. Thus we now create a function which will save a video and embed it in an html string to create a video of the systems development

In [8]:
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):
    system.integrator.run(50)

    # 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), end="\r")
    return part,

We now can start the sampling over the animation class of matplotlib

In [9]:
fig, ax = plt.subplots(figsize=(10,10))
part, = ax.plot([],[], 'o')

animation.FuncAnimation(fig, run, frames=100, blit=True, interval=0, repeat=False, init_func=init)
progress: 100%
Out[9]: