This goal of this project is to simulate and run experiments on the lattice structure of Europium-Oxide (EuO), a face-center cubic (FCC) crystal. This is version is in the making and shall function as the main module for all experiments. The first task is to fix two bugs; randomisation of lattice site and magnetisation values. The lattice site [shown below] are not randomised as it should, and the magnetisation values are much higher than excpeted.
Once the issues are resovled, we will plot [as a function of temperature] the system's magnetisation and hamiltonian, along with their averages and errors.
In [9]:
from copy import deepcopy
from math import pi, cos, sin, sqrt
import numpy as np
import operator
import random
nearest_neighbor_delta = [(-1, -1, 0), (-1, 0, -1), (-1, 0, 1), (-1, 1, 0),
(0, -1, -1), (0, -1, 1), (0, 1, -1), (0, 1, 1),
(1, -1, 0), (1, 0, -1), (1, 0, 1), (1, 1, 0)]
next_nearest_neighbor_delta = [(-2, 0, 0), (2, 0, 0),
(0, -2, 0), (0, 2, 0),
(0, 0, -2), (0, 0, 2)]
class Site:
def __init__(self, index, coords, theta=pi, phi=0):
self.index = index
self.coords = coords
self.theta = theta # spin
self.phi = phi # spin
self.nearest = [] # neighbors
self.next_nearest = [] # neighbors
def neighbors(self):
n = list(self.nearest)
n.extend(self.next_nearest)
return n
def with_attr(self, attr, value):
s = deepcopy(self)
setattr(s, attr, value)
return s
def __repr__(self):
return "Site({}, {}): {} | {}".format(
self.index, self.coords,
map(lambda n: n.index, self.nearest),
map(lambda n: n.index, self.next_nearest))
class Simulator:
def __init__(self, size):
self.size = size
self.num_sites = self.get_num_sites(size)
self.sites = []
self.coord_to_site = {}
self.magnetizations = []
self.records = {}
# Parameters
self.J_NN = 1
self.J_NNN = 1
self.H_EXTERNAL = 1
self.MAX_DELTA_THETA = 1
self.MAX_DELTA_PHI = 1
self.NUM_STEPS = 50
self.SIM_THRESHOLD = 0
self.SIM_SAMPLE_RES = 5
self.TEMP_LOWER = 0
self.TEMP_UPPER = 5
self.TEMP_STEP_SIZE = 1
self.num_temp_steps = ((self.TEMP_UPPER - self.TEMP_LOWER) /
self.TEMP_STEP_SIZE) + 1
self.populate_sites() # initialize lattice structure
def reset(self):
"""
Resets this simulator to a ready state.
"""
self.sites = []
self.coord_to_site = {}
self.populate_sites()
def run(self):
"""
Runs the simulation using the initialization parameters.
"""
self.reset()
for temp in np.linspace(self.TEMP_LOWER, self.TEMP_UPPER, self.num_temp_steps):
print 'evaluating temp: {}'.format(temp)
self.magnetizations = []
for step in range(self.NUM_STEPS): # num steps of simulation
beta = self.J_NN / temp
# Evaluate each site
for i in range(self.num_sites):
self.evaluate_and_update_site_spin(self.sites[i], beta)
# Sample from every nth step after the threshold.
if step > self.SIM_THRESHOLD and step % self.SIM_SAMPLE_RES == 0:
self.magnetizations.append(self.get_magnetization())
# Record state for each step of temperature
self.records[temp] = self.get_state()
def get_magnetization(self):
"""
Calculate the magnetization value for the system and
save it to self.magnetizations.
"""
#S = 7.0 / 2.0
x = reduce(lambda a, site: a + sin(site.theta) + cos(site.phi),
self.sites, 0)
y = reduce(lambda a, site: a + sin(site.theta) + sin(site.phi),
self.sites, 0)
z = reduce(lambda a, site: a + cos(site.theta),
self.sites, 0)
#x *= S
#y *= S
#z *= S
return sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2))
def get_state(self):
"""
Calculate magnetic moment and errors
"""
mags = self.magnetizations
num_mags = len(self.magnetizations)
# averages
m_avg = reduce(operator.add, mags, 0) / num_mags
# average errors
m_error = reduce(lambda a, m: a + pow(m - m_avg, 2), mags, 0)
m_error = sqrt(m_error / (num_mags * (num_mags - 1)))
return dict(mag_avg=m_avg, m_error=m_error)
def evaluate_and_update_site_spin(self, site, beta):
"""
Changes the given site's spin slightly and calculates the
change in hamiltonian. If the change falls within the threshold,
the site's new spin is saved.
"""
new_site = deepcopy(site)
# Create a new site with only its theta changed
delta_theta = self.MAX_DELTA_THETA * random.uniform(-1, 1)
site_theta = site.with_attr('theta', site.theta + delta_theta)
if random.random() < -beta * self.delta_hamiltonian(site, site_theta):
new_site.theta = site_theta.theta
# Create a new site with only its phi changed
delta_phi = self.MAX_DELTA_PHI * random.uniform(-1, 1)
site_phi = site.with_attr('phi', delta_phi)
if random.random() < -beta * self.delta_hamiltonian(site, site_phi):
new_site.phi = site_phi.phi
self.sites[site.index] = new_site
def delta_hamiltonian(self, old_site, new_site):
"""
Calculate the difference in Hamiltonian values
for the two systems.
"""
old_hamiltonian = 0.0
new_hamiltonian = 0.0
def hamiltonian(site, b):
def Hamlet(accum, neighbor):
return accum - b * \
(sin(site.theta) *
sin(neighbor.theta) *
cos(site.phi - neighbor.phi) +
cos(site.theta) *
cos(neighbor.theta))
return Hamlet
old_hamiltonian = reduce(hamiltonian(old_site, self.J_NN), # reduction function
old_site.nearest, # target neighbors
0) # starting value
old_hamiltonian = reduce(hamiltonian(old_site, self.J_NNN), # reduction function
old_site.next_nearest, # target neighbors
old_hamiltonian) # starting value
new_hamiltonian = reduce(hamiltonian(new_site, self.J_NN),
new_site.nearest,
0)
new_hamiltonian = reduce(hamiltonian(new_site, self.J_NNN),
new_site.next_nearest,
new_hamiltonian)
return new_hamiltonian - old_hamiltonian
def populate_sites(self):
"""
Creates a map of site indices to arrays of site indices
that represent the nearest and next-nearest neighbors for
each site.
"""
origin = Site(0, (0, 0, 0))
self.save_site(origin)
seen = [origin.index]
queue = [origin]
while queue:
s = self.add_neighbors(queue.pop(0))
s = self.save_site(s)
for n in s.neighbors():
if n.index not in seen:
seen.append(n.index)
queue.append(n)
def add_neighbors(self, site):
"""
Updates the site with lists of its nearest and next-nearest
neighbor sites.
"""
def tuple_add(a, b):
"""
Add two tuples together
Ex: (0, 0, 1) + (1, 0, 1) = (1, 0, 2)
See http://stackoverflow.com/a/497894
"""
return tuple(map(operator.add, a, b))
def add_neighbor_site(neighbor_delta, neighbor_sites_list):
"""
Create the neighboring site using the delta coords
and add it to the given list
"""
coords = tuple_add(neighbor_delta, site.coords)
neighbor_site = self.get_site(coords)
neighbor_sites_list.append(neighbor_site)
for neighbor_delta in nearest_neighbor_delta:
add_neighbor_site(neighbor_delta, site.nearest)
for neighbor_delta in next_nearest_neighbor_delta:
add_neighbor_site(neighbor_delta, site.next_nearest)
return site
def get_site(self, coords):
"""
Either fetches the site with the given coords
or creates one.
Makes sure coords are within the bounds of the fcc
and wraps if they are not.
get_site((0, 0, 0)) -> Site 0
get_site((1, 1, 1)) -> Site 5
get_site((-1, -1, -1)) -> Site 3
get_site((2, 2, 2)) -> Site 3
get_site((5, 5, 5)) -> Site 3
"""
coords = tuple(map(lambda c: c % self.size, coords))
if coords in self.coord_to_site:
return self.coord_to_site[coords]
site = Site(len(self.sites), coords)
return self.save_site(site)
def save_site(self, site):
# Saves the site to our list of sites
if site.index >= len(self.sites):
site.index = len(self.sites)
self.sites.append(site)
else:
self.sites[site.index] = site
self.coord_to_site[site.coords] = site
return site
def get_num_sites(self, size):
"""
Get the number of sites in a (n x n x n) fcc,
where n = size.
Explanation: When retrieving the nearest neighbor
for a given site, indexing out of the bounds of the
cube should wrap around and always provide a valid
site. Therefore, when looking at a 1x1x1 fcc, only
a single corner site along with 3 face sites are
necessary to model, since the rest of the sites
are effectively redundant, since all of the corner
sites are actually the same site, just wrapped.
"""
return pow(size, 3) * 4
s = Simulator(1) # Simulator(size); size = n, in (n x n x n)
s.run()
from pprint import pprint
pprint(s.records)
pprint(s.sites)
In [ ]: