An ensemble of trajectories obtained from simulating Langevin dynamics will tend to a stable distribution: the Boltzmann distribution.
Two identical magnetic nanoparticles, aligned along their anisotropy axes. The system has 6 degrees of freedom (x,y,z components of magnetisation for each particle) but the energy is defined by the two angles $\theta_1,\theta_2$ alone.
The Boltzmann distribution represents the probabability that the system will be found within certain sets (i.e. the angles of magnetisation). The distribution is parameterised by the temperature of the system and the energy landscape of the problem.
Note the sine terms appear because the distribution is over the solid angles. In other words, the distribution is over the surface of a unit sphere. The sine terms project these solid angles onto a simple elevation angle between the magnetisation and the anisotropy axis ($\theta$).
$$p\left(\theta_1,\theta_2,\phi_1,\phi_2\right) = \frac{\sin(\theta_1)\sin(\theta_2)e^{-E\left(\theta_1,\theta_2,\phi_1,\phi_2\right)/\left(K_BT\right)}}{Z}$$The energy function for a single domain magnetic nanoparticle is given by the Stoner-Wohlfarth equation:
$$\frac{E\left(\theta_1,\theta_2,\phi_1,\phi_2\right)}{K_BT}=\sigma\left(\sin^2\theta_1+\sin^2\theta_2\right) -\nu\left[2\cos\theta_1\cos\theta_2 - \sin\theta_1\sin\theta_2\cos\left(\phi_1-\phi_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 [1]:
import numpy as np
In [2]:
# dipole interaction energy
def dd(t1, t2, p1, p2, nu):
return -nu*(2*np.cos(t1)*np.cos(t2) - np.sin(t1)*np.sin(t2)*np.cos(p1-p2))
# anisotropy energy
def anis(t1, t2, sigma):
return sigma*(np.sin(t1)**2 + np.sin(t2)**2)
# total energy
def tot(t1, t2, p1, p2, nu, sigma):
return dd(t1, t2, p1, p2, nu) + anis(t1, t2, sigma)
In [3]:
# numerator of the Boltzmann distribution (i.e. ignoring the partition function Z)
def p_unorm(t1,t2,p1,p2,nu,sigma):
return np.sin(t1)*np.sin(t2)*np.exp(-tot(t1, t2, p1, p2, nu, sigma))
In [4]:
# non interacting
from scipy.integrate import nquad
sigma, nu = 1.0, 0.0
Z = nquad(
lambda t1, t2, p1, p2: p_unorm(t1, t2, p1, p2, nu, sigma),
ranges=[(0, np.pi), (0, np.pi), (0, 2*np.pi), (0, 2*np.pi)]
)
print(Z[0])
In [5]:
Z = nquad(
lambda t1, t2: p_unorm(t1, t2, 0, 0, nu, sigma),
ranges=[(0, np.pi), (0, np.pi)]
)
print(Z[0] * 4 * np.pi**2)
In [6]:
import magpy as mp
In [7]:
K = 1e5
r = 7e-9
T = 330
Ms=400e3
R=9e-9
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 [64]:
base_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
)
ensemble = mp.EnsembleModel(50000, base_model)
Simulate an ensemble of 10,000 dimers without interactions.
In [65]:
res = ensemble.simulate(end_time=1e-9, time_step=1e-12,
max_samples=500, random_state=1002,
n_jobs=8, implicit_solve=True,
interactions=False)
m_z0 = np.array([state['z'][0] for state in res.final_state()])/Ms
m_z1 = np.array([state['z'][1] for state in res.final_state()])/Ms
theta0 = np.arccos(m_z0)
theta1 = np.arccos(m_z1)
System magnetisation shows that the system has relaxed into the local minima (we could relax the system globally but it would take much longer to run since the energy barrier must be overcome).
In [66]:
import matplotlib.pyplot as plt
%matplotlib inline
In [67]:
plt.plot(res.results[0].time, res.ensemble_magnetisation())
plt.title('Non-interacting dimer ensemble magnetisation');
In [68]:
# Dimensionless parameters
V = 4./3*np.pi*r**3
sigma = K*V/mp.core.get_KB()/T
nu = 0
print('Sigma: {:.3f}'.format(sigma))
print(' Nu: {:.3f}'.format(nu))
The joint distribution for both angles are computed analytically and compared with the numerical result.
The resulting distrubiton is symmetric because both particles are independent and identically distributed.
In [69]:
Z = nquad(
lambda t1, t2, p1, p2: p_unorm(t1, t2, p1, p2, nu, sigma),
ranges=[(0, np.pi/2), (0, np.pi/2), (0, 2*np.pi), (0, 2*np.pi)]
)
print(Z[0])
Z=Z[0]
In [70]:
ts = np.linspace(min(theta0), max(theta0), 100)
In [71]:
bdist = [[
nquad(lambda p1, p2: p_unorm(t1, t2, p1, p2, nu, sigma)/Z, ranges=[(0, 2*np.pi), (0, 2*np.pi)])[0]
for t1 in ts]
for t2 in ts]
In [72]:
plt.hist2d(theta0, theta1, bins=30, normed=True);
plt.contour(ts, ts, bdist, cmap='Greys')
plt.title('Joint distribution')
plt.xlabel('$\\theta_1$'); plt.ylabel('$\\theta_2$');
We can also compute the marginal distribution (i.e. the equilibrium of just 1 particle). It is easier to see the alignment of the two distributions.
In [73]:
b_marginal = [nquad(
lambda t2, p1, p2: p_unorm(t1, t2, p1, p2, nu, sigma)/Z,
ranges=[(0, np.pi/2), (0, 2*np.pi), (0, 2*np.pi)])[0]
for t1 in ts]
In [74]:
plt.hist(theta0, bins=50, normed=True)
plt.plot(ts, np.array(b_marginal))
Out[74]:
In [75]:
# Dimensionless parameters
V = 4./3*np.pi*r**3
sigma = K*V/mp.core.get_KB()/T
nu = mp.core.get_mu0() * V**2 * Ms**2 / 2.0 / np.pi**2 / R**3 / mp.core.get_KB() / T
print('Sigma: {:.3f}'.format(sigma))
print(' Nu: {:.3f}'.format(nu))
The interaction strength is very strong (actually the particles are impossibly close). The following command is identical to above except that interactions=True
In [76]:
res = ensemble.simulate(end_time=1e-9, time_step=1e-13,
max_samples=500, random_state=1001,
n_jobs=8, implicit_solve=False,
interactions=True, renorm=True)
m_z0i = np.array([state['z'][0] for state in res.final_state()])/Ms
m_z1i = np.array([state['z'][1] for state in res.final_state()])/Ms
theta0i = np.arccos(m_z0i)
theta1i = np.arccos(m_z1i)
In [77]:
plt.plot(res.results[0].time, res.ensemble_magnetisation())
plt.title('Interacting dimer ensemble magnetisation');
The stationary distributions align BUT:
introduced fudge factor of $\pi$ into the denominator of the interaction constant $\nu$. In other words we use $$\nu_\textrm{fudge}=\frac{\nu}{\pi}$$
This factor of $1/\pi$ could come from integrating somewhere. I think this is an error with my analytic calculations and not the code. This is because the code certainly uses the correct term for interaction strength, whereas I derived this test myself.
In [78]:
# Dimensionless parameters
V = 4./3*np.pi*r**3
sigma = K*V/mp.core.get_KB()/T
nu = 1.0 * mp.core.get_mu0() * V**2 * Ms**2 / 2.0 / np.pi / np.pi / R**3 / mp.core.get_KB() / T
print('Sigma: {:.3f}'.format(sigma))
print(' Nu: {:.3f}'.format(nu))
In [79]:
Z = nquad(
lambda t1, t2, p1, p2: p_unorm(t1, t2, p1, p2, nu, sigma),
ranges=[(0, np.pi/2), (0, np.pi/2), (0, 2*np.pi), (0, 2*np.pi)]
)[0]
In [80]:
ts = np.linspace(min(theta0), max(theta0), 100)
bdist = [[
nquad(lambda p1, p2: p_unorm(t1, t2, p1, p2, nu, sigma)/Z, ranges=[(0, 2*np.pi), (0, 2*np.pi)])[0]
for t1 in ts]
for t2 in ts]
In [81]:
plt.hist2d(theta0i, theta1i, bins=30, normed=True);
# ts = np.linspace(min(theta0), max(theta0), 100)
# b = boltz_2d(ts, nu, sigma)
plt.contour(ts, ts, bdist, cmap='Greys')
plt.title('Joint distribution')
plt.xlabel('$\\theta_1$'); plt.ylabel('$\\theta_2$');
We use the marginal distribution again to check the convergence. We also compare to the interacting case
In [82]:
b_marginal = [nquad(
lambda t2, p1, p2: p_unorm(t1, t2, p1, p2, nu, sigma)/Z,
ranges=[(0, np.pi/2), (0, 2*np.pi), (0, 2*np.pi)])[0]
for t1 in ts]
In [83]:
plt.hist(theta0i, bins=50, normed=True, alpha=0.6, label='Magpy + inter.')
plt.hist(theta0, bins=50, normed=True, alpha=0.6, label='Magpy + no inter.')
plt.plot(ts, b_marginal, label='Analytic')
plt.legend();
One way to test this could be to simulate a 3D system. If another factor of $\pi$ appears then it is definitely something missing from my analytic calculations.
In [47]:
import pymc3 as pm
with pm.Model() as model:
z1 = pm.Uniform('z1', 0, 1)
theta1 = pm.Deterministic('theta1', np.arccos(z1))
z2 = pm.Uniform('z2', 0, 1)
theta2 = pm.Deterministic('theta2', np.arccos(z2))
phi1 = pm.Uniform('phi1', 0, 2*np.pi)
phi2 = pm.Uniform('phi2', 0, 2*np.pi)
energy = tot(theta1, theta2, phi1, phi2, nu, sigma)
like = pm.Potential('energy', -energy)
In [52]:
with model:
step = pm.NUTS()
trace = pm.sample(500000, step=step)
In [54]:
pm.traceplot(trace)
Out[54]:
In [61]:
plt.hist(trace['theta1'], bins=200, normed=True);
plt.plot(ts, b_marginal, label='Analytic')
Out[61]:
In [63]:
plt.hist(trace['theta1'], bins=200, normed=True, alpha=0.6);
plt.hist(theta0i, bins=50, normed=True, alpha=0.6, label='Magpy + inter.');
In [ ]:
def cart_energy(mx1, my1, mz1, mx2, my2, my3, nu, sigma):
t1 = np.arccos(mz1)
t2 = np.arccos(mz2)
anis = sigma*np.sin(t1)**2 + sigma*np.sin(t2)**2
inter = -nu*(3*np.dot)
In [ ]:
import pymc3 as pm
with pm.Model() as model:
z1 = pm.Uniform('z1', 0, 1)
theta1 = pm.Deterministic('theta1', np.arccos(z1))
z2 = pm.Uniform('z2', 0, 1)
theta2 = pm.Deterministic('theta2', np.arccos(z2))
phi1 = pm.Uniform('phi1', 0, 2*np.pi)
phi2 = pm.Uniform('phi2', 0, 2*np.pi)
energy = tot(theta1, theta2, phi1, phi2, nu, sigma)
like = pm.Potential('energy', -energy)