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.
In this part we want to calculate estimators for the initial susceptibility, i.e. the susceptibility at zero external magnetic field. One could carry out several simulations with different external magnetic field strengths and get the initial susceptibility by fitting a line to the results. We want to go a more elegant way by using fluctuation formulas known from statistical mechanics. In three dimensions the initial susceptibility $\chi_{init}$ can be calculated with zero field simulations through
\begin{equation} \chi_{init} = \frac{V \cdot \mu_0}{3 \cdot k_B T} \left( \langle \boldsymbol{M}^2 \rangle - \langle \boldsymbol{M} \rangle^2 \right) = \frac{\mu_0}{3 \cdot k_B T \cdot V} \left( \langle \boldsymbol{\mu}^2 \rangle - \langle \boldsymbol{\mu} \rangle^2 \right) \end{equation}where $\boldsymbol{M}$ is the magnetization vector and $\boldsymbol{\mu}$ is the total magnetic dipole moment of the system. In direction $i$ it reads
\begin{equation} M_i = \frac{1}{V} \Bigg\langle \sum_{j=1}^N \tilde{\mu}_j^i \Bigg\rangle \end{equation}where $\tilde{\mu}_j^i$ is the $j$ th dipole moment in direction $i$.
We want to derive the fluctuation formula. We start with the definition of the magnetic susceptibility. In general this reads
\begin{equation} \chi \equiv \frac{\partial}{\partial H} \langle M_{\boldsymbol{H}} \rangle \end{equation}with $\langle M_{\boldsymbol{H}} \rangle$ the ensemble averaged magnetization in direction of a homogeneous external magnetic field $\boldsymbol{H}$.
In thermal equilibrium the ensemble average of the magnetization reads
\begin{equation} \langle M_{\boldsymbol{H}} \rangle = \frac{1}{V Z_c} \left \lbrack \sum_{\alpha} \mu_{\boldsymbol{H},\alpha} e^{ -\beta E_{\alpha}(H=0) + \beta\mu_0\mu_{\boldsymbol{H},\alpha}H }\right \rbrack \end{equation}with $Z_c$ the canonical partition function, $E_{\alpha}(H=0)$ the energy without an external magnetic field $\boldsymbol{H}$, $\beta$ the inverse thermal energy $\frac{1}{k_BT}$, $\mu_{\boldsymbol{H},\alpha}$ the total magnetic dipole moment of the system in direction of the external magnetic field $\boldsymbol{H}$ in microstate $\alpha$ and $V$ the system volume.
Now we insert the magnetization $\langle M_{\boldsymbol{H}} \rangle$ in the definition of the magnetic susceptibility $\chi$ and let the derivative operate on the ensemble average. We get the fluctuation formula
\begin{equation} \chi = \frac{\beta\mu_0}{V} \left \lbrack \frac{1}{Z_c}\sum_{\alpha} \mu_{\alpha}^2~ e^{ -\beta E_{\alpha}(H=0) + \beta\mu_0\mu_{\boldsymbol{H},\alpha}H } - \frac{1}{Z_c}\sum_{\alpha} \mu_{\alpha}~ e^{ -\beta E_{\alpha}(H=0) + \beta\mu_0\mu_{\boldsymbol{H},\alpha}H }~~ \frac{1}{Z_c}\sum_{\alpha'}\mu_{\alpha'}~ e^{ -\beta E_{\alpha'}(H=0) + \beta\mu_0\mu_{\boldsymbol{H},\alpha}H }\right \rbrack = \frac{\beta\mu_0}{V} \left \lbrack \langle \mu_{\boldsymbol{H}}^2 \rangle - \langle \mu_{\boldsymbol{H}} \rangle^2 \right \rbrack = \frac{\beta\mu_0}{V} \left(\Delta \mu_{\boldsymbol{H}}\right)^2 \end{equation}At zero external magnetic field ($H = 0$) there is no distinct direction for the system, so we can take the fluctuations $\Delta \mu$ in all directions and divide it by the dimension. Thus we can use more data points of our simulation for the average and get a more precise estimator for the susceptibility. Thus finally the fluctuation formula for the initial susceptibility in three dimensions reads
\begin{equation} \chi_{init} = \frac{\beta\mu_0}{3V} \left \lbrack \langle \boldsymbol{\mu}^2 \rangle - \langle \boldsymbol{\mu} \rangle^2 \right \rbrack = \frac{V\beta\mu_0}{3} \left \lbrack \langle \boldsymbol{M}^2 \rangle - \langle \boldsymbol{M} \rangle^2 \right \rbrack \end{equation}where $\boldsymbol{\mu}$ and $\boldsymbol{M}$ are defined above.
In this part we want to consider a three dimensional ferrofluid system and compare our result for the initial susceptibility $\chi_{init}$ with them of Ref. [1].
First we import all necessary packages and check for the required Espresso features
In [ ]:
import espressomd
espressomd.assert_features('DIPOLES', 'LENNARD_JONES')
from espressomd.magnetostatics import DipolarP3M
import numpy as np
Now we set up all necessary simulation parameters
In [ ]:
lj_sigma=1
lj_epsilon=1
lj_cut = 2**(1./6.) * lj_sigma
# magnetic field constant
mu_0 = 1.
# Particles
N = 1000
# Volume fraction
# phi = rho * 4. / 3. * np.pi * ( lj_sigma / 2 )**3.
phi = 0.0262
# Dipolar interaction parameter lambda = mu_0 m^2 /(4 pi sigma^3 kT)
dip_lambda = 3.
# Temperature
kT =1.0
# Friction coefficient
gamma=1.0
# Time step
dt =0.02
# box size 3d
box_size = (N * np.pi * 4./3. * (lj_sigma / 2.)**3. / phi)**(1./3.)
and 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. Latter means that the initial configuration of our system is the same every time this script is executed. As the time evolution of the system depends not solely on the Langevin thermostat but also on the numeric accuracy and DP3M (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 [ ]:
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 system volume
pos = box_size * np.random.random((N,3))
# Add particles
system.part.add(pos=pos, rotation=N*[(1,1,1)], dip=dip)
# 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.skin = 0.8
# Setup dipolar P3M
accuracy = 5E-4
system.actors.add(DipolarP3M(accuracy=accuracy,prefactor=dip_lambda*lj_sigma**3*kT))
Now we equilibrate for a while
In [ ]:
print("Equilibration...")
equil_rounds = 10
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 need the magnetization of our system, we import MagneticDipoleMoment from observables which returns us the total dipole moment of the system which is the magnetization times the volume of the system.
In [ ]:
from espressomd.observables import MagneticDipoleMoment
dipm_tot_calc = MagneticDipoleMoment(ids=system.part[:].id)
Now we set the desired number of loops for the sampling
In [ ]:
# Sampling
loops = 2000
and sample the first and second moment of the magnetization or total dipole moment, by averaging over all total dipole moments occurring during the simulation
In [ ]:
print('Sampling ...')
# calculate initial total dipole moment
dipm_tot_temp = dipm_tot_calc.calculate()
# initiate variables for sum of first and second moment of total dipole moment
dipm_tot_sum = dipm_tot_temp
dipm_tot_2_sum = np.square(dipm_tot_temp)
# sample dipole moment
for i in range(loops):
system.integrator.run(10)
dipm_tot_temp = dipm_tot_calc.calculate()
# add current dipole moment
dipm_tot_sum = np.sum((dipm_tot_sum, dipm_tot_temp), axis=0)
dipm_tot_2_sum = np.sum((dipm_tot_2_sum, np.square(dipm_tot_temp)), axis=0)
# print progress only every 10th cycle
if (i+1)%10 == 0:
print("progress: {:3.0f}%".format((i+1)*100./loops), end="\r")
print("\nSampling done")
# calculate average first and second moment of total dipole moment
dipm_tot = dipm_tot_sum / loops
dipm_tot_2 = dipm_tot_2_sum / loops
For the estimator of the initial susceptibility $\chi_{init}$ we need the magnitude of one single dipole moment
In [ ]:
# dipole moment
dipm = np.sqrt(dip_lambda*4*np.pi*lj_sigma**3.*kT / mu_0)
print("dipm = {}".format(dipm))
Now we can calculate $\chi_{init}$ from our simulation data
In [ ]:
# susceptibility in 3d system
chi = mu_0 /(system.volume() * 3. * kT) * ( np.sum(dipm_tot_2 * dipm**2.) - np.sum(np.square(dipm_tot * dipm)) )
and print the result
In [ ]:
print('chi = %.4f' % chi)
Compared with the value $\chi = 0.822 \pm 0.017$ of Ref. [1] (see table 1) it should be very similar.
Now we want to compare the result with the theoretical expectations. At first with the simple Langevin susceptibility
In [ ]:
chi_L = 8. * dip_lambda * phi
print('chi_L = %.4f' % chi_L)
and at second with the more advanced one (see Ref. [1] eq. (6)) which has a cubic accuracy in $\chi_L$ and reads
\begin{equation} \chi = \chi_L \left( 1 + \frac{\chi_L}{3} + \frac{\chi_L^2}{144} \right) \end{equation}
In [ ]:
chi_I = chi_L * ( 1 + chi_L / 3. + chi_L**2. / 144. )
print('chi_I = %.4f' % chi_I)
Both of them should be smaller than our result, but the second one should be closer to our one. The deviation of the theoretical results to our simulation result can be explained by the fact that in the Langevin model there are no interactions between the particles incorporated at all and the more advanced (mean-field-type) one of Ref. [1] do not take occurring cluster formations into account but assumes a homogeneous distribution of the particles. For higher values of the volume fraction $\phi$ and the dipolar interaction parameter $\lambda$ the deviations will increase as the cluster formation will become more pronounced.
At the end of this tutorial we now want to sample the magnetization curve of a three dimensional system and compare the results with analytical solutions. Again we will compare with the Langevin function but also with the approximation of Ref. [2] (see also Ref. [1] for the right coefficients) which takes the dipole-dipole interaction into account. For this approximation, which is a modified mean-field theory based on the pair correlation function, the Langevin parameter $\alpha$ is replaced by
\begin{equation} \alpha' = \alpha + \chi_L~L(\alpha) + \frac{\chi_L^{2}}{16} L(\alpha) \frac{d L(\alpha)}{d\alpha} \end{equation}where $\chi_L$ is the Langevin susceptibility
\begin{equation} \chi_L = \frac{N}{V}\frac{\mu_0 \mu^2}{3k_BT} = 8 \cdot \lambda \cdot \phi \end{equation}Analogous to part II we start at zero external magnetic field and increase the external field successively. At every value of the external field we sample the total dipole moment which is proportional to the magnetization as we have a fixed volume.
First we create a list of values of the Langevin parameter $\alpha$. As we already sampled the magnetization at zero external field in the last section, we take this value and continue with the sampling of an external field unequal zero
In [ ]:
alphas = [0.25, 0.5, 1, 2, 3, 4, 8]
Now for each value in this list we sample the total dipole moment / magnetization of the system for a while. Have in mind that we have only committed the orientation of the dipole moments to the particles. Thus we have to commit $H\cdot \mu$ as the external magnetic field to ESPResSo, where $\mu$ is the magnitude of a single magnetic dipole moment.
As in part II we use the same system for every value of the Langevin parameter $\alpha$. Thus we use that the system is already pre-equilibrated from the previous run so we save some equilibration time. For scientific purposes one would use a new system for every value for the Langevin parameter to ensure that the systems are independent and no correlation effects are measured. Also one would perform more than just one simulation for each value of $\alpha$ to increase the precision of the results.
In [ ]:
# remove all constraints
system.constraints.clear()
# list of magnetization in field direction
magnetization = []
# append result for alpha=0 from previous chapter
magnetization.append(np.average(dipm_tot))
# number of loops for sampling
loops_m = 500
for alpha in alphas:
print("Sample for alpha = {}".format(alpha))
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)
print("done\n")
# Equilibration
print("Equilibration...")
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\n")
# Sampling
print("Sampling...")
magn_temp = 0
for i in range(loops_m):
system.integrator.run(20)
magn_temp += dipm_tot_calc.calculate()[0]
print("progress: {:3.0f}%".format((i+1)*100./loops_m), end="\r")
print("\n")
# save average magnetization
magnetization.append(magn_temp / loops_m)
print("Sampling for alpha = {} done \n".format(alpha))
print("magnetizations = {}".format(magnetization))
print("total progress: {:5.1f}%\n".format((alphas.index(alpha)+1)*100./len(alphas)))
# remove constraint
system.constraints.clear()
print("Magnetization curve sampling done")
Now we define the Langevin function and the modified mean-field-approximation of the Langevin parameter of Ref. [2]
In [ ]:
# Langevin function
def L(y):
return np.cosh(y)/np.sinh(y)-1/y
In [ ]:
# second order mean-field-model from Ref. [2]
def alpha_mean_field(alpha, dip_lambda, phi):
chi = 8. * dip_lambda * phi
return alpha + chi * L(alpha) + chi**2. / 16. * L(alpha) * ( 1./(alpha**2.) - 1./((np.sinh(alpha))**2.) )
We also want to plot the linear approximation at $\alpha = 0$ to see for which values of $\alpha$ this approximation holds. We use the initial susceptibility calculated in the first chapter of this part as the gradient. As we want the gradient of $M^*$ with respect to $\alpha$ which fulfills the relation
\begin{equation} \frac{\partial M^*}{\partial \alpha} = \frac{1}{M_{sat}}\frac{\partial M}{\partial \left( \frac{\mu_0\mu}{k_BT} H\right)} = \frac{k_BT~V}{\mu_0\mu^2N}\frac{\partial M}{\partial H} = \frac{k_BT~V}{\mu_0\mu^2N}~\chi \end{equation}we have to scale our calculated initial susceptibility $\chi_{init}$ by a factor to get it in our dimensionless units.
Now we plot the resulting curves together with our simulation results and the linear approximation
In [ ]:
import matplotlib.pyplot as plt
In [ ]:
y = np.arange(0.01,10, 0.1, dtype=float).tolist()
alphas.insert(0, 0.0)
L_func = []
L_mean_field = []
init_susceptibility = []
for i in y:
L_func.append(L(i))
L_mean_field.append(L(alpha_mean_field(i, dip_lambda, phi)))
init_susceptibility.append(i*system.volume()*kT/(N*mu_0*dipm**2)*chi)
# divide all entries in the magnetization list by N to get the dimensionless magnetization
magnetization_star = []
for i in range(len(magnetization)):
magnetization_star.append(magnetization[i] / N)
print('magnetization_star = {}'.format(magnetization_star))
plt.figure(figsize=(10,10))
plt.ylim(0, 1.)
plt.xlabel(r'$\alpha$', fontsize=20)
plt.ylabel(r'$M^*$', fontsize=20)
plt.plot(y, L_func, label='Langevin function')
plt.plot(y, L_mean_field, label='modified mean-field-theory')
plt.plot(alphas, magnetization_star, 'o', label='simulation results')
plt.plot(y, init_susceptibility, label=r'linear approximation at $\alpha = 0$')
plt.legend(fontsize=20)
plt.show()
We can see that the magnetization curve where we used the Langevin parameter of the modified mean-field-theory is closer to our simulation results. Also we can clearly see that the linear approximation holds only for very small values of $\alpha$.
At this point is should be mentioned, that the modified mean-field-model assumes a spatially homogeneous system which is not the case at higher volume fractions $\phi$ and dipolar interaction parameters $\lambda$ as the particles form clusters. We can already see this with our simulation results as they visibly deviate from the modified mean-field-model.
At sufficiently high volume fractions $\phi$ and dipolar interaction parameters $\lambda$ these clusters can be thus rigid, that simulation with normal methods are impossible as the relaxation times exceeds normal simulation times by far resulting in strongly correlated configurations and thus measurements.
[1] Zuowei Wang, Christian Holm, and Hanns Walter Müller. “Molecular dynamics study on the equilibrium magnetization properties and structure of ferrofluids”. In: Phys. Rev. E 66 (2 Aug. 2002), p. 021405. doi: 10.1103/PhysRevE.66.021405. url: https://link.aps.org/doi/10.1103/PhysRevE.66.021405.
[2] Alexey O. Ivanov and Olga B. Kuznetsova. “Magnetic properties of dense ferrofluids: An influence of interparticle correlations”. In: Phys. Rev. E 64 (4 Sept. 2001), p. 041405. doi: 10.1103/PhysRevE.64.041405. url: https://link.aps.org/doi/10.1103/PhysRevE.64.041405.