import matplotlib.pyplot as plt
import matplotlib.animation as ani
import numpy as np
from numpy.random import default_rng
import copy
[docs]def print_progress(step, total):
"""
Prints a progress bar.
Args:
step (int): progress counter
total (int): counter at completion
"""
message = "simulation progress ["
total_bar_length = 60
percentage = int(step / total * 100)
bar_fill = int(step / total * total_bar_length)
for i in range(total_bar_length):
if i < bar_fill:
message += "|"
else:
message += " "
message += "] "+str(percentage)+" %"
if step < total:
print(message, end="\r")
else:
print(message)
[docs]def generate_system(lattice_size, cluster=False):
"""
Creates the system.
The system is represented by a 2D :math:`N \\times N` lattice
where each lattice site
may be either occupied (1) or unoccupied (0).
The system may be a surface, in which case one of its sides
has a row of empty sites, or a cluster, in which case all
sides have empty sites.
Args:
lattice_size (int): lattice size :math:`N`
cluster (bool): if True, the system is a cluster
Returns:
array: the system as an integer array
"""
grid = np.array( [ [1]*lattice_size ]*lattice_size )
if cluster: # create a cluster by clearing sites on all sides
grid[:,0] = 0
grid[:,-1] = 0
grid[0,:] = 0
grid[-1,:] = 0
else: # create a surface by clearing sites on one side
grid[-1,:] = 0
return grid
[docs]def show_system(grid):
"""
Draws the system.
Args:
grid (array): the system as an integer array
"""
plt.clf()
ax = plt.axes()
ax.set_aspect('equal')
plt.pcolormesh(grid, cmap='Greys')
plt.show()
[docs]def draw(frame, history):
"""
Draws the system for animation.
Args:
frame (int): index of the frame to draw
history (list): list of systems at different times
"""
plt.clf()
ax = plt.axes()
ax.set_aspect('equal')
plt.pcolormesh(history[frame], cmap='Greys')
[docs]def animate(history):
"""
Animate the simulation.
Args:
history (list): list of systems at different times
"""
nframes = len(history)
print("animating "+str(nframes)+" frames")
fig = plt.figure()
motion = ani.FuncAnimation(fig, draw, nframes, fargs=( history, ) )
plt.show()
[docs]def count_neighbors(i,j,grid):
"""
Count the number of occupied neighbors for the lattice site (i,j).
The system is represented as a square lattice so the first neighbors of
a given site are the 4 sites directly next to it
left, right, up and down.
The second neighbors of a site are the 4 sites diagonally next to it:
up-left, up-right, down-left, down-right.
This function counts how many of those sites are occupied (have the value 1).
The edges of the system wrap around.
It means that for a site that is, say, on the left edge of the system,
the left-side neighbor is a site on the right edge of the system.
Args:
i (int): row index
j (int): column index
grid (array): the system as an integer array
Returns:
int, int: number of first neighbors, number of second neighbors
"""
lattice_size = len(grid)
left = (j-1)%lattice_size # i-1 (mod n). This is n-1, if i = 0.
right = (j+1)%lattice_size # i+1 (mod n). This is 0, if i = n-1.
up = (i+1)%lattice_size
# The bottom layer is supposed to be bulk matter, so we count the atoms there
# to be their own neighbors. This makes them harder to remove so that surfaces
# are really etched from the top.
down = max( (i-1), 0 )
firsts = grid[i,left] + grid[i,right] + grid[up,j] + grid[down,j]
seconds = grid[up,left] + grid[up,right] + grid[down,left] + grid[down,right]
return firsts, seconds
[docs]def get_removal_rate( i, j, grid, unremovables = [] ):
"""
Calculates the removal rate for a site.
For normal sites, the removal rate is evaluated based on
how many 1st and 2nd neighbors the site has.
These are calculated with :meth:`count_neighbors`.
Certain sites may be listed as protected sites that may never be
removed. Such sites mimic the use of protective masks.
.. note ::
Edit to change the physics!
Args:
i (int): row index
j (int): column index
grid (array): the system as an integer array
unremovables (list): indices of sites that are never removed
Return:
float: the removal rate
"""
for protected in unremovables:
if i == protected[0] and j == protected[1]:
return 0
# If this is a normal site, we evaluate its rate
# of being removed according to how many neighbors it has.
n_1st, n_2nd = count_neighbors(i,j,grid)
# There are a few possible sets of removal rates below.
# kinks are much more soluble than planes,
# horizontal and diagonal planes are equally soluble
rates_A = np.array(
[ [10**20, 10**15, 10**10, 1],
[10**15, 10**10, 10**5, 0],
[10**10, 10**5, 100, 0],
[1, 0, 0, 0] ] )
# kinks are much more soluble than planes,
# horizontal planes are more soluble than diagonal
rates_B = np.array(
[ [10**20, 10**15, 10**10, 1],
[10**15, 10**10, 1000, 0],
[10**10, 10**5, 100, 0],
[1, 0, 0, 0] ] )
# kinks are much more soluble than planes,
# horizontal planes are less soluble than diagonal
rates_C = np.array(
[ [10**20, 10**15, 10**10, 1],
[10**15, 10**10, 10**5, 0],
[10**10, 1000, 100, 0],
[1, 0, 0, 0] ] )
# kinks and planes are equally soluble
rates_D = np.array(
[ [10**20, 10**15, 10**10, 1],
[10**15, 10**5, 10**5, 0],
[10**10, 10**5, 100, 0],
[1, 0, 0, 0] ] )
# Choose which set of rates to use:
rates = rates_A
return rates[n_1st-1, n_2nd-1]
[docs]def count_probabilities(grid):
"""
Calculates, for each site, the probability that this particular site is
the next to be removed.
The function loops over all occupied sites :math:`(i,j)` and calculates
the removal rate :math:`r_{i,j}` for each site.
The total removal rate is defined as
.. math ::
R = \\sum_{i,j} r_{i,j}
and the probability that site :math:`(i,j)` is the next to be removed, is
.. math ::
p(i,j) = \\frac{1}{R} r_{i,j}.
Finally, the function lists the coordinates of all occupied sites in a
single site vector :math:`S` and constructs the corresponding vectors
containing the site specific probabilities :math:`p` and accumulated
probabilities :math:`P`.
The elements of vector :math:`S` are pairs of integers, so you may have,
e.g., :math:`S_0 = (0,0), S_1 = (1,0), S_2 = (2,0)` etc.
The elements of vector :math:`p` are probabilities, so you may have,
e.g., :math:`p_0 = p(0,0), p_1 = p(1,0), p_2 = p(2,0)` etc.
The elements of vector :math:`P` are sums of probabilities,
.. math::
P_n = \\sum_{k = 0}^n p_k.
So :math:`P_0 = p_0, P_1 = p_0 + p_1, P_2 = p_0 + p_1 + p_2` etc.
Especially the probabilities must sum up to 1, so the last
element of vector :math:`P` is always 1.
Args:
grid (array): the system as an integer array
Returns:
array, array, array: :math:`p, P, S`
"""
accumulated_probabilities = []
probabilities = []
coordinates = []
total_rate = 0.0
# Loop over all sites except the bottom row.
# We assume the bottom row is connected to bulk
# and will not be removed.
for i in range(lattice_size-1):
for j in range(lattice_size):
# only consider occupied sites
if grid[i,j] > 0:
# calculate the rate and sum of rates
removal_rate = get_removal_rate( i, j, grid )
total_rate += removal_rate
# only consider sites that can be removed
if removal_rate > 0:
# save the rate, sum of rates and current coordinates
probabilities.append( removal_rate )
accumulated_probabilities.append( total_rate )
coordinates.append( [i,j] )
# Divide rates by the final total rate to get probabilities.
accumulated_probabilities = np.array( accumulated_probabilities ) / total_rate
probabilities = np.array( probabilities ) / total_rate
coordinates = np.array( coordinates )
return probabilities, accumulated_probabilities, coordinates
[docs]def pick_random_event(grid):
"""
Randomly choose the next event (removal of a site) to happen.
The random event is chosen as follows:
* The probability :math:`p` for choosing each site is calculated
with :meth:`count_probabilities`.
* Also the accumulated probabilities :math:`P` are calculated.
* A random number :math:`R \\in [0,1]` is drawn.
* The event :math:`n` for which :math:`P_{n-1} < R \le P_{n}` is chosen.
(For :math:`n = 0`, it is enough that :math:`R \le P_{0}`.)
* The function returns the coordinates of the site for this event, :math:`S_n`.
.. note ::
This function is incomplete!
Args:
grid (array): the system as an integer array
Returns:
array: coordinates :math:`(i,j)` of the chosen site
"""
randomizer = random.random()
prob, acc_prob, coords = count_probabilities(grid)
return coords[ int( randomizer*len(coords) ) ] # todo
[docs]def randomly_remove_atom(grid):
"""
Randomly choose one site for removal and make it unoccupied.
The site is chosen using :meth:`pick_random_event`.
The chosen site is marked as unoccupied by
setting the grid value to 0.
Args:
grid (array): the system as an integer array
"""
atom_coordinates = pick_random_event(grid)
grid[atom_coordinates[0], atom_coordinates[1]] = 0
[docs]def main(lattice_size, n_steps, n_plots, show_as_you_go = False, cluster = False):
"""
Main program.
Creates a system as a :math:`N \\times N` square lattice,
then removes atoms from it.
Args:
lattice_size (int): lattice size :math:`N`
n_steps (int): the total number of atoms to remove
n_plots (int): the number of images to create
show_as_you_go (bool): If True, will show images of the
system during simulation. If False, an animation is
created at the end. Only set to True for short tests.
cluster (bool): if True, the system is a cluster
"""
grid = generate_system(lattice_size, cluster)
history = [grid]
for i in range(n_steps):
randomly_remove_atom(grid)
print_progress(i+1, n_steps)
# take a look at the system every now and then
if i%(n_steps/n_plots) == 0:
# show progress as simulation proceeds
if show_as_you_go:
show_system(grid)
else:
history.append(copy.copy(grid))
if not show_as_you_go:
animate(history)
if __name__ == "__main__":
lattice_size = 40 # system size
n_steps = 100 # simulation length
n_plots = 50 # number of images to produce
show_as_you_go = False # draw the system during simulation or at the end?
random = default_rng()
main(lattice_size, n_steps, n_plots, show_as_you_go)