In this notebook we simulate the thermal equilibrium (Boltzmann distrubion) of two interacting magnetic nanoparticles (dimer), coupled with dipolar interactions. The equilibrium distribution is computed with three different approaches
Two particles are aligned along their anisotropy axes at a distance of $R$. They have identical volume $V$, anisotropy constant $K$, saturation magnetisation $M_s$. The temperature of the environment is $T$. The angle of the magnetic moment to the anisotropy axis is $\theta_1,\theta_2$ for particle 1 and 2 respectively.
The solid angle of the magnetisation angles follow a Boltzmann distribution, such that the angle $\theta$ is distributed:
$$p\left(\theta_1,\theta_2\right) = \frac{\sin(\theta_1)\sin(\theta_2)e^{-E\left(\theta_1,\theta_2\right)/\left(K_BT\right)}}{Z}$$where
$$\frac{E\left(\theta_1,\theta_2\right)}{K_BT}=\sigma\left(\cos^2\theta_1+\cos^2\theta_2\right) -\nu\left(3\cos\theta_1\cos\theta_2 - \cos\left(\theta_1-\theta_2\right)\right)$$$$\sigma=\frac{KV}{K_BT}$$$$\nu=\frac{\mu_0V^2M_s^2}{2\pi R^3K_BT}$$$\sigma,\nu$ are the normalised anisotropy and interaction strength respectively.
In [16]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import ImageGrid
%matplotlib inline
Individual energy terms
In [17]:
def dd(t1, t2, nu):
return -nu*(3*np.cos(t1)*np.cos(t2) - np.cos(t1-t2))
def anis(t1, t2, sigma):
return sigma*(np.sin(t1)**2 + np.sin(t2)**2)
def tot(t1, t2, nu, sigma):
return dd(t1, t2, nu) + anis(t1, t2, sigma)
The unnormalised probability of state $\theta_1\theta_2$
In [18]:
def p_unorm(t1,t2,nu,sigma):
return np.sin(t1)*np.sin(t2)*np.exp(-tot(t1,t2,nu,sigma))
2-dimensional Boltzmann distribution
In [19]:
from scipy.integrate import dblquad
def boltz_2d(ts, nu, sigma):
e = np.array([[p_unorm(t1,t2,nu,sigma) for t1 in ts] for t2 in ts])
Z = dblquad(lambda t1,t2: p_unorm(t1,t2,nu,sigma),
0, ts[-1], lambda x: 0, lambda x: ts[-1])[0]
return e/Z
In [20]:
nus = [0, 0.3, 0.3]
sigmas = [2.0, 2.0, 0.5]
In [21]:
ts = np.linspace(0, np.pi, 100)
fg = plt.figure(figsize=(10,4))
axs = ImageGrid(
fg, 111, nrows_ncols=(1,3), axes_pad=0.15,
share_all=True,cbar_location="right",
cbar_mode="single",cbar_size="7%",
cbar_pad=0.15,
)
for nu, sigma, ax in zip(nus, sigmas, axs):
e = [[tot(t1, t2, nu, sigma) for t1 in ts] for t2 in ts]
cf=ax.contourf(ts, ts, e)
ax.set_xlabel('$\\theta_1$'); ax.set_ylabel('$\\theta_2$')
ax.set_aspect('equal')
axs[0].set_title('No interactions')
axs[1].set_title('Strong interactions')
axs[2].set_title('Weak anisotropy')
ax.cax.colorbar(cf) # fix color bar
Out[21]:
The probability of each point in the phase space is directily related to the energy through the Boltzmann distribution (above). The three different cases are described:
In [22]:
ts = np.linspace(0, np.pi, 100)
fg = plt.figure(figsize=(10,4))
axs = ImageGrid(
fg, 111, nrows_ncols=(1,3), axes_pad=0.15,
share_all=True,cbar_location="right",
cbar_mode="single",cbar_size="7%",
cbar_pad=0.15,
)
for nu, sigma, ax in zip(nus, sigmas, axs):
b = boltz_2d(ts, nu, sigma)
cf=ax.contourf(ts, ts, b)
ax.set_xlabel('$\\theta_1$'); ax.set_ylabel('$\\theta_2$')
ax.set_aspect('equal')
axs[0].set_title('No interactions')
axs[1].set_title('Strong interactions')
axs[2].set_title('Weak anisotropy')
ax.cax.colorbar(cf) # fix color bar
Out[22]:
Computing the analytic solutions required the partition function $Z$ to be computed, which was achieved with numerical integration. For more than 2 particles this becomes extremely difficult. A more efficient approach is to use MCMC. In MCMC we randomly choose a system state and add it to a histogram depending on how energetically favourable the state is. The more states we try, the closer our histogram becomes to the true distribution.
We use PyMC3 for MCMC sampling. The model is simple:
PyMC3 uses the NUTS algorithm to efficicently sample trial states with a high acceptance rate.
In [23]:
import pymc3 as pm
Setting up the PyMC3 model is simple! Specify the priors and the energy function using pm.Potential
In [24]:
nu = 0.3
sigma = 1.5
with pm.Model() as model:
z1 = pm.Uniform('z1', -1, 1)
theta1 = pm.Deterministic('theta1', np.arccos(z1))
z2 = pm.Uniform('z2', -1, 1)
theta2 = pm.Deterministic('theta2', np.arccos(z2))
energy = tot(theta1, theta2, nu, sigma)
like = pm.Potential('energy', -energy)
Chose the NUTS algorithm as our MCMC step and request a large number of random samples. These are returned in the trace variable.
In [25]:
with model:
step = pm.NUTS()
trace = pm.sample(100000, step=step)
In [26]:
b = boltz_2d(ts, nu, sigma)
plt.hist2d(trace['theta1'], trace['theta2'], bins=70, normed=True)
plt.contour(ts, ts, b, cmap='Greys')
plt.gca().set_aspect('equal')
plt.xlabel('$\\theta_1$'); plt.ylabel('$\\theta_2$');
The results look good! We can also check that the marginal distributions of the two angles are correct. The marginal distributions are also easier to check by eye.
$$p(\theta_1) = \int p(\theta_1, \theta_2) \mathrm{d}\theta_2$$$$p(\theta_2) = \int p(\theta_1, \theta_2) \mathrm{d}\theta_1$$
In [27]:
from scipy.integrate import trapz
b_marginal = -b.sum(axis=0) / trapz(ts, b.sum(axis=0))
In [28]:
fg, axs = plt.subplots(ncols=2, figsize=(9,3))
for i in range(2):
axs[i].hist(trace['theta{}'.format(i+1)], bins=100, normed=True)
axs[i].plot(ts, b_marginal, lw=2)
axs[i].set_xlabel('$\\theta_{}$'.format(i+1))
axs[i].set_ylabel('$p(\\theta_{})$'.format(i+1))
plt.suptitle('Marginal distributions', fontsize=18)
Out[28]:
In [29]:
import magpy as mp
In [30]:
K = 1e5
r = 8e-9
T = 330
Ms=400e3
R=1.5e-3
kdir = [0, 0, 1]
location1 = np.array([0, 0, 0], dtype=np.float)
location2 = np.array([0, 0, R], dtype=np.float)
direction = np.array([0, 0, 1], dtype=np.float)
alpha = 1.0
In [31]:
V = 4./3*np.pi*r**3
sigma = K*V/mp.get_KB()/T
nu = mp.get_mu0() * V**2 * Ms**2 / 2.0 / np.pi / R**3 / mp.get_KB() / T
print('Sigma: {:.3f}'.format(sigma))
print(' Nu: {:.3f}'.format(nu))
In [32]:
model = mp.Model(
anisotropy=np.array([K, K], dtype=np.float),
anisotropy_axis=np.array([kdir, kdir], dtype=np.float),
damping=alpha,
location=np.array([location1, location2], dtype=np.float),
magnetisation=Ms,
magnetisation_direction=np.array([direction, direction], dtype=np.float),
radius=np.array([r, r], dtype=np.float),
temperature=T
)
In [33]:
res = model.simulate_ensemble(end_time=1e-9, time_step=1e-12,
max_samples=500, seeds=range(5000),
n_jobs=8, implicit_solve=True,
interactions=False)
m_z = [state['z'][0] for state in res.final_state()]
theta = np.arccos(m_z)
In [19]:
ts = np.linspace(0, np.pi/2, 100)
b = boltz_2d(ts, nu, sigma)
b_marginal = -b.sum(axis=0) / trapz(ts, b.sum(axis=0))
b_noint = b_marginal
plt.hist(theta, bins=50, normed=True)
plt.plot(ts[ts<theta.max()], b_marginal[ts<theta.max()])
Out[19]:
In [715]:
res = model.simulate_ensemble(end_time=1e-9, time_step=1e-12,
max_samples=500, seeds=range(5000),
n_jobs=8, implicit_solve=True,
interactions=True)
m_z = [state['z'][0] for state in res.final_state()]
theta = np.arccos(m_z)
In [716]:
ts = np.linspace(0, np.pi/2, 100)
b = boltz_2d(ts, nu, sigma)
b_marginal = -b.sum(axis=0) / trapz(ts, b.sum(axis=0))
plt.hist(theta, bins=50, normed=True)
plt.plot(ts[ts<theta.max()], b_marginal[ts<theta.max()])
plt.plot(ts[ts<theta.max()], b_noint[ts<theta.max()])
Out[716]:
In [717]:
R=1e-7
location1 = np.array([0, 0, 0], dtype=np.float)
location2 = np.array([0, 0, R], dtype=np.float)
In [718]:
nu = mp.get_mu0() * V**2 * Ms**2 / 2.0 / np.pi / R**3 / mp.get_KB() / T
print('Sigma: {:.3f}'.format(sigma))
print(' Nu: {:.3f}'.format(nu))
In [719]:
model = mp.Model(
anisotropy=np.array([K, K], dtype=np.float),
anisotropy_axis=np.array([kdir, kdir], dtype=np.float),
damping=alpha,
location=np.array([location1, location2], dtype=np.float),
magnetisation=Ms,
magnetisation_direction=np.array([direction, direction], dtype=np.float),
radius=np.array([r, r], dtype=np.float),
temperature=T
)
In [720]:
res = model.simulate_ensemble(end_time=1e-9, time_step=1e-12,
max_samples=500, seeds=range(5000),
n_jobs=8, implicit_solve=True,
interactions=True)
m_z = [state['z'][0] for state in res.final_state()]
theta = np.arccos(m_z)
In [721]:
ts = np.linspace(0, np.pi/2, 100)
b = boltz_2d(ts, nu, sigma)
b_marginal = -b.sum(axis=0) / trapz(ts, b.sum(axis=0))
plt.hist(theta, bins=50, normed=True)
plt.plot(ts[ts<theta.max()], b_marginal[ts<theta.max()])
plt.plot(ts[ts<theta.max()], b_noint[ts<theta.max()])
Out[721]:
In [36]:
R=0.6e-8
location1 = np.array([0, 0, 0], dtype=np.float)
location2 = np.array([0, 0, R], dtype=np.float)
In [37]:
nu = mp.get_mu0() * V**2 * Ms**2 / 2.0 / np.pi / R**3 / mp.get_KB() / T
print('Sigma: {:.3f}'.format(sigma))
print(' Nu: {:.3f}'.format(nu))
In [38]:
model = mp.Model(
anisotropy=np.array([K, K], dtype=np.float),
anisotropy_axis=np.array([kdir, kdir], dtype=np.float),
damping=alpha,
location=np.array([location1, location2], dtype=np.float),
magnetisation=Ms,
magnetisation_direction=np.array([direction, direction], dtype=np.float),
radius=np.array([r, r], dtype=np.float),
temperature=T
)
In [39]:
res = model.simulate_ensemble(end_time=10e-9, time_step=1e-12,
max_samples=500, seeds=range(5000),
n_jobs=8, implicit_solve=True,
interactions=True)
m_z = [state['z'][0] for state in res.final_state()]
theta = np.arccos(m_z)
In [40]:
ts = np.linspace(0, np.pi/2, 100)
b = boltz_2d(ts, nu, sigma)
b_marginal = -b.sum(axis=0) / trapz(ts, b.sum(axis=0))
plt.hist(theta, bins=50, normed=True, label='simulation')
plt.plot(ts[ts<theta.max()], b_marginal[ts<theta.max()], label='analytic')
plt.plot(ts[ts<theta.max()], b_noint[ts<theta.max()], label='analytic ($\\nu=0$)')
plt.legend()
Out[40]:
In [783]:
r = 1e-9
In [784]:
V = 4./3*np.pi*r**3
sigma = K*V/mp.get_KB()/T
nu = mp.get_mu0() * V**2 * Ms**2 / 2.0 / np.pi / R**3 / mp.get_KB() / T
print('Sigma: {:.3f}'.format(sigma))
print(' Nu: {:.3f}'.format(nu))
In [788]:
model = mp.Model(
anisotropy=np.array([K, K], dtype=np.float),
anisotropy_axis=np.array([kdir, kdir], dtype=np.float),
damping=alpha,
location=np.array([location1, location2], dtype=np.float),
magnetisation=Ms,
magnetisation_direction=np.array([direction, direction], dtype=np.float),
radius=np.array([r, r], dtype=np.float),
temperature=T
)
In [795]:
res = model.simulate(end_time=1e-9, time_step=1e-13, max_samples=1000, seed=1001,
implicit_solve=True, interactions=True)
res.plot();
In [798]:
res = model.simulate_ensemble(end_time=1e-9, time_step=1e-13,
max_samples=500, seeds=range(5000),
n_jobs=8, implicit_solve=True,
interactions=True)
m_z = [state['z'][0] for state in res.final_state()]
theta = np.arccos(m_z)
In [799]:
ts = np.linspace(0, np.pi, 100)
b = boltz_2d(ts, nu, sigma)
b_marginal = -b.sum(axis=0) / trapz(ts, b.sum(axis=0))
plt.hist(theta, bins=50, normed=True)
plt.plot(ts, b_marginal)
Out[799]:
In [ ]: