<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Thermal equilibrium of interacting dimer
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
 - Analytic solution is computed by hand (easy for two aligned particles)
 - Markov chain Monte-Carlo (MCMC) sampling of the energy functions using PyMC3
 - Simulating the stochastic Landau-Lifshitz-Gilbert equation until equilibrium using magpy
Setup
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.
Analytic solution
The solid angle of the magnetisation angles follow a Boltzmann distribution, such that the angle $\theta$ is distributed
Step1:  Individual energy terms
Step2:  The unnormalised probability of state $\theta_1\theta_2$
Step3:  2-dimensional Boltzmann distribution
Step4:  Analytical results
We show resulting calculations for the following three cases
Step5:  Energy landscape
The plots below show the energy associated with each point of the 2-dimensional phase space (i.e. one angle per particle). Low energy (purple) states are energetically favourable wherea the high energy states (yellow) are not.
Step6:  Probability state-space
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
Step7:  Markov-chain Monte-Carlo (MCMC)
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
Step8:  Setting up the PyMC3 model is simple! Specify the priors and the energy function using pm.Potential
Step9:  Chose the NUTS algorithm as our MCMC step and request a large number of random samples. These are returned in the trace variable.
Step10:  Compare results
The trace contains a large number of samples of the system states. The distribution (histogram) over this large sample of states should be close to the true distribution. We can compare the 2D histogram of the system angles to the analytic solution we computed previously.
Step11:  Marginalisation over one variable
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$$
Step12:  Langevin dynamics simulations (sLLG)
Simulate the stochastic Landau-Lifshitz-Gilbert equation using MagPy.
By generating an ensemble of trajectories until thermal equilibriation we can approximate the equilibrium distribution.
Non-interacting case
Step13:  Negligible interactions case
Step14:  Weakly interacting case
Step15:  Strongly interacting case
Step16:  High noise
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import ImageGrid
%matplotlib inline
Explanation: Thermal equilibrium of interacting dimer
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
 - Analytic solution is computed by hand (easy for two aligned particles)
 - Markov chain Monte-Carlo (MCMC) sampling of the energy functions using PyMC3
 - Simulating the stochastic Landau-Lifshitz-Gilbert equation until equilibrium using magpy
Setup
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.
Analytic solution
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.
End of explanation
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)
Explanation: Individual energy terms
End of explanation
def p_unorm(t1,t2,nu,sigma):
    return np.sin(t1)*np.sin(t2)*np.exp(-tot(t1,t2,nu,sigma))
Explanation: The unnormalised probability of state $\theta_1\theta_2$
End of explanation
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
Explanation: 2-dimensional Boltzmann distribution
End of explanation
nus = [0, 0.3, 0.3]
sigmas = [2.0, 2.0, 0.5]
Explanation: Analytical results
We show resulting calculations for the following three cases:
 - No interactions
 - Strong interactions
 - Strong interactions and weak anisotropy
End of explanation
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
Explanation: Energy landscape
The plots below show the energy associated with each point of the 2-dimensional phase space (i.e. one angle per particle). Low energy (purple) states are energetically favourable wherea the high energy states (yellow) are not.
End of explanation
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
Explanation: Probability state-space
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:
 - With no interactions, the probability space is symetrical since both variables are identically distibuted and independent of one another. Both particle are much more likely to be found close to the anisotropy axis (either parallel or antiparallel).
 - With interactions, the symmetry is broken and the particles prefer to be aligned with the anisotropy axis and one another. This increases the probability of the system being in aligned states around the axes rather.
 - With weak anisotropy, the particles prefer to be aligned with one another but are less likely to be found near the anisotropy axis.
End of explanation
import pymc3 as pm
Explanation: Markov-chain Monte-Carlo (MCMC)
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:
We sample uniformly on the unit sphere (i.e. draw solid angles uniformly), converting into the angle $\theta_1$ and $\theta_2$
We compute the energy for both particles
We pass this energy to PyMC3
PyMC3 uses the NUTS algorithm to efficicently sample trial states with a high acceptance rate.
End of explanation
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)
Explanation: Setting up the PyMC3 model is simple! Specify the priors and the energy function using pm.Potential
End of explanation
with model:
    step = pm.NUTS()
    trace = pm.sample(100000, step=step)
Explanation: Chose the NUTS algorithm as our MCMC step and request a large number of random samples. These are returned in the trace variable.
End of explanation
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$');
Explanation: Compare results
The trace contains a large number of samples of the system states. The distribution (histogram) over this large sample of states should be close to the true distribution. We can compare the 2D histogram of the system angles to the analytic solution we computed previously.
End of explanation
from scipy.integrate import trapz
b_marginal = -b.sum(axis=0) / trapz(ts, b.sum(axis=0))
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)
Explanation: Marginalisation over one variable
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$$
End of explanation
import magpy as mp
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
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))
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
)
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)
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()])
Explanation: Langevin dynamics simulations (sLLG)
Simulate the stochastic Landau-Lifshitz-Gilbert equation using MagPy.
By generating an ensemble of trajectories until thermal equilibriation we can approximate the equilibrium distribution.
Non-interacting case
End of explanation
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)
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()])
Explanation: Negligible interactions case
End of explanation
R=1e-7
location1 = np.array([0, 0, 0], dtype=np.float)
location2 = np.array([0, 0, R], dtype=np.float)
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))
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
)
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)
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()])
Explanation: Weakly interacting case
End of explanation
R=0.6e-8
location1 = np.array([0, 0, 0], dtype=np.float)
location2 = np.array([0, 0, R], dtype=np.float)
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))
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
)
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)
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()
Explanation: Strongly interacting case
End of explanation
r = 1e-9
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))
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
)
res = model.simulate(end_time=1e-9, time_step=1e-13, max_samples=1000, seed=1001,
                     implicit_solve=True, interactions=True)
res.plot();
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)
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)
Explanation: High noise
End of explanation
<END_TASK>