<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
An ensemble of trajectories obtained from simulating Langevin dynamics will tend to a stable distribution
Step1:  Magpy non-interacting case
Using magpy we initialise the dimers with both magnetisation vectors aligned along their anisotropy axes. We allow the system to relax.
Step2:  System properties
Set up the parameters of the dimers. They are identical and aligned along their anisotropy axes
Step3:  Magpy model and simulation
Build a magpy model of the dimer
Step4:  Simulate an ensemble of 10,000 dimers without interactions.
Step5:  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).
Step6:  Compare to analytic thermal equilibrium
Compute the expected boltzmann distribution
Step7:  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.
Step8:  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.
Step9:  Magpy interacting case
We now simulate the exact same ensemble of dimers but with the interactions enabled. We can compute the dimensionless parameters to understand the strength of the interactions (vs. the anisotropy strength).
Step10:  The interaction strength is very strong (actually the particles are impossibly close). The following command is identical to above except that interactions=True
Step11:  System relaxation
The system quickly relaxes into the first minima again, as before.
Step12:  Thermal equilibrium
The stationary distributions align BUT
Step13:  We use the marginal distribution again to check the convergence. We also compare to the interacting case
Step14:  Possible sources of error
Implementation of interaction strength in magpy (but there are many tests against the true equations)
Analytic calculations for equilibrium
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.
Step15:  Cartesian
<ASSISTANT_TASK:>
Python Code:
import numpy as np
# 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)
# 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))
# 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])
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)
Explanation: Thermal equilibrium
An ensemble of trajectories obtained from simulating Langevin dynamics will tend to a stable distribution: the Boltzmann distribution.
Problem setup
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.
Boltzmann distribution
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}$$
Stoner-Wohlfarth model
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.
Functions for analytic solution
End of explanation
import magpy as mp
Explanation: Magpy non-interacting case
Using magpy we initialise the dimers with both magnetisation vectors aligned along their anisotropy axes. We allow the system to relax.
End of explanation
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
Explanation: System properties
Set up the parameters of the dimers. They are identical and aligned along their anisotropy axes
End of explanation
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)
Explanation: Magpy model and simulation
Build a magpy model of the dimer
End of explanation
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)
Explanation: Simulate an ensemble of 10,000 dimers without interactions.
End of explanation
import matplotlib.pyplot as plt
%matplotlib inline
plt.plot(res.results[0].time, res.ensemble_magnetisation())
plt.title('Non-interacting dimer ensemble magnetisation');
Explanation: 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).
End of explanation
# 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))
Explanation: Compare to analytic thermal equilibrium
Compute the expected boltzmann distribution
End of explanation
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]
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]
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$');
Explanation: 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.
End of explanation
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]
plt.hist(theta0, bins=50, normed=True)
plt.plot(ts, np.array(b_marginal))
Explanation: 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.
End of explanation
# 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))
Explanation: Magpy interacting case
We now simulate the exact same ensemble of dimers but with the interactions enabled. We can compute the dimensionless parameters to understand the strength of the interactions (vs. the anisotropy strength).
End of explanation
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)
Explanation: The interaction strength is very strong (actually the particles are impossibly close). The following command is identical to above except that interactions=True
End of explanation
plt.plot(res.results[0].time, res.ensemble_magnetisation())
plt.title('Interacting dimer ensemble magnetisation');
Explanation: System relaxation
The system quickly relaxes into the first minima again, as before.
End of explanation
# 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))
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]
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]
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$');
Explanation: Thermal equilibrium
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.
End of explanation
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]
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();
Explanation: We use the marginal distribution again to check the convergence. We also compare to the interacting case
End of explanation
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)
with model:
    step = pm.NUTS()
    trace = pm.sample(500000, step=step)
pm.traceplot(trace)
plt.hist(trace['theta1'], bins=200, normed=True);
plt.plot(ts, b_marginal, label='Analytic')
plt.hist(trace['theta1'], bins=200, normed=True, alpha=0.6);
plt.hist(theta0i, bins=50, normed=True, alpha=0.6, label='Magpy + inter.');
Explanation: Possible sources of error
Implementation of interaction strength in magpy (but there are many tests against the true equations)
Analytic calculations for equilibrium
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.
End of explanation
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)
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)
Explanation: Cartesian
End of explanation
<END_TASK>