# MonteCoffee Documentation¶

Welcome to the documentation of MonteCoffee!

For a general idea of what the code does and how it works, see Code overview.

For a minimal working example of how to perform kinetic Monte Carlo with MonteCoffee, see Quick start.

## Installation¶

To download MonteCoffee , simply clone the git:

git clone https://gitlab.com/ChemPhysChalmers/MonteCoffee.git


https://gitlab.com/ChemPhysChalmers/MonteCoffee/-/archive/master/MonteCoffee-master.zip


### Setup¶

MonteCoffee is presently installed, simply by copying (or linking) the base-file folder to the run-directory which contains the user-specific files.

To verify that the code runs successfully on your current workstation, you can run:

python3 test.py


in one of the tutorial folders.

### Dependencies¶

MonteCoffee depends on the following Python packages:
• Numpy used for array manipulations.
• Matplotlib for plotting results.
• ASE for handling atomic structures.
• h5py for handling file output and compression
• mpi4py for running various kMC simulations in parallel

## Code overview¶

### What is kinetic Monte Carlo?¶

Kinetic Monte Carlo (kMC) is a simulation technique that can be used to investigate the kinetics of chemical reactions. Kinetics can be seen as transitions between different chemical states, which obeys the master equation:

\begin{equation} \dfrac{\text{d}P_\alpha}{\text{d}t} = \sum_\beta W_{\alpha\beta}P_\beta - W_{\beta\alpha}P_\alpha \\ \end{equation}

where $$\alpha, \beta$$ are the states defined by the site-occupations (e.g. CO on site 1, CO on site 2, site 3 empty ,…) , $$W_\alpha\beta$$ is the transition rate from state $$\beta$$ to state $$\alpha$$, and $$P_\alpha$$ is the probability for being in state $$P_\alpha$$. The equation defines a system of coupled differential equations, with one equation for each $$\alpha$$.

KMC solves this system of equations by randomly generating transitions between states. The transitions are generated by reactive events, which for example can be $$\mathrm{O_2}$$ dissociative adsorption proceeding on sites number 1 and 3, where site 1 and 3 are neighboring sites. The time of occurrence of a reactive event (i) is in MonteCoffee generated according to the first-reaction method:

\begin{equation} t^\text{occ}_i = t^\text{gen}_i-\dfrac{\text{ln}\,u}{k_i},\quad u \in [0,1[ \\ \end{equation}

where $$t^\text{occ}_i$$ is the time of occurrence, $$t^\text{gen}_i$$ is the time the event was generated (simulation time), $$k_i$$ is the rate constant of the reaction-step, and $$u$$ is a random uniform deviate.

### Structure of modules¶

MonteCoffee is written using object-oriented programming. The central modules are defined as shown in the figure below: The modules have different functions:
• Site: Defines a catalytic site with properties such as site-type, coordination numbers, and nearest neighbors.
• System: Defines the collection of catalytic sites.
• Events: Defines the possible events in the simulation, and their rate constants.
• NeighborKMC: Controls the actual simulation. Therefore, it is the central object of the simulation.
• Basin: Defines the acceleration in the simulation based on rates. (N.B., in the journal article M. Jørgensen and H. Grönbeck J. Chem. Phys. 149, 114101 (2018) this module belonged to NeighborKMC, but it was later separated as the program grew.)
• Logging: Handles logging of simulation by receiving input from NeighborKMC.

The site-connectivity is defined by each site’s neighborlist. The set of neighborlists for all sites defines a global connectivity pattern.

For more information about the classes, modules, and methods, please see User interface (API). The API is vital as MonteCoffee is designed as a programmable application. Thus, the user downloads the modules and prepares a simulation by making changes to the files named user_*.py.

### Implicit assumptions¶

MonteCoffee has a few implicit assumptions:

• The user masters the concept of object-oriented programming in Python.
• The chemical species are simply represented as integers for computational efficiency. The user decides the meaning of each integer.
• At most two sites are involved in binding adsorbates and reactions. (Coarse-grained sites can be assumed).
• Only sites that are in each others’ neighbor-list are connected.
• The event numbering is decided by the order of which the user loads the events (see the example in test.py).
• The model implemented by the user is thermodynamically consistent, and detailed balance is obeyed by the events.

### Code output¶

#### Simulation log file¶

MonteCoffee generates a .txt log file with the naming convention kMClog_#date_#time_.txt. This is handled by the module base.Logging.

The first output contains all the dictionary of parameters, passed as the input parameters when instantiating a NeighborKMC object. This section may look as:

Mikkel Jorgensen (2015-2019)
Noemi Bosio (since 2019)
Elisabeth M. Dietze (since 2020)
Chalmers University of Technology
Goteborg, Sweden
--------------------------------------------------------------------------------
Simulation parameters
pCO       :    2000.0
pO2       :    1000.0
T         :    800.0
Name      :    COOx Simulation
reverses  :    {0: 1, 2: 3, 4: 4, 5: 5}
tend      :    1e-07
Nsites    :    432
--------------------------------------------------------------------------------
No time acceleration used.


Next the log begins with printing the simulation step, the time the log-point was dumped, the simulation time, and the number of events fired, respectively:

Step          time[hr:min:s]             Sim time [s]          Events called
1000             14:47:40             4.84677160781e-10        ['37', '0', '3', '0', '960', '1', '0']
2000             14:47:41             7.10054219544e-10        ['50', '0', '4', '0', '1945', '2', '0']
3000             14:47:42             8.66046768567e-10        ['59', '0', '4', '0', '2936', '2', '0']


The ordering of Events Called is determined by the order of which the parameter events are passed upon instantiating a NeighborKMC object.

#### Results as .txt files¶

MonteCoffee generates its output as plain-text (.txt) file for basic output and as object in the HDF5 format for more advanced output, as determined by the method save_txt() in base.NeighborKMC.

If desired, this method can be altered to output other information during runtime. The output files are generated using numpy’s savetxt() method, and can therefore be loaded using np.loadtxt().

The code outputs the following .txt files once prior to simulation:

• siteids.txt: The index of each site, passed as the parameter ind when instantiating a NeighborKMC.user_sites.Site object. This is useful for storing an Ase.Atoms object.
• stypes.txt: The site types for each site, passed as the parameter stype when instantiating a NeighborKMC.user_sites.Site object.

During the simulation the following .txt files are updated with a frequency specified in Options:

• mcstep.txt: The Monte Carlo step corresponding to each line in the other .txt files. This is directly dependent on the updating frequency.

• time.txt: The simulation time in seconds for every logged Monte Carlo step.

• coverages.txt: The coverages at each time-step for each lattice-site. The coverages are structured as coverages[mcstep][site_number].

• evs_exec.txt: The total number of event-type executions. For example, to find the total number of executions of event number 0:

evs_exec = np.loadtxt("evs_exec.txt")
N1 = evs_exec

• detail_site_event_evol.hdf5: The detailed time evolution of the system, which saves every step. The file can be read in as follows:

with h5py.File('detail_site_event_evol.hdf5','a') as f2:
time_list = np.array(f2["time"])
event_list = np.array(f2["event"])
site_list = np.array(f2["site"])
othersite_list = np.array(f2["othersite"])

Working with this file, one has to be careful, because it can be rather big, despite the compression. In principle hdf5 can be read also with any other hdf5 supporting programming language.

For further information about analyzing output, see Analyze results and Calculating a turnover frequency.

## Tutorials¶

The tutorials assume that the different elements of MonteCoffee are split into the user_* files. All files described here are provided in the Tutorial folder which is downloaded from GIT. Additionally we provide information on how to run parallel simulations using MonteCoffee and give some ideas about special site rules for more complex systems.

### Quick start¶

This tutorial provides a minimal working example of how to run MonteCoffee. For simplicity this example puts all definitions in one file. This structure is different from the following tutorials.

For this very first, simple example, all free energy barriers are assumed constant. Two different events are implemented, and neighbors are calculated from inter-atomic distances in an ASE.Atoms object. For simplicity, this guide shows how to setup a simulation in a single python script. In principle it is similar to the second tutorial on diatomic adsorption and desorption (see Dissociative Adsorption).

Step 1. Implement the two event-types

In this step we import the template class NeighborKMC.base.events.EventBase, and derive two classes from this to enable two different types of reactions. All defined events must implement four methods:

The first event class is an adsorption, that is possible if the pair of neighbor sites is empty. Each eventclass can be given a name which is only used in the output. The adsorption class has a rate-constant that is linearly dependent on the pressure, and if executed, it covers the sites with species 1:

from base.events import EventBase

def __init__(self, params):

def possible(self, system, site, other_site):

if system.sites[site].covered == 0 and system.sites[other_site].covered == 0:
return True
else:
return False

def get_rate(self, system, site, other_site):
R = 1. * self.params["pA"]
return R

def do_event(self, system, site, other_site):
# Cover the two sites with species 1
system.sites[site].covered = 1
system.sites[other_site].covered = 1

def get_involve_other(self):
return True


Now we define the reverse desorption-event with a constant rate:

class Desorption(EventBase):

def __init__(self, params):
EventBase.__init__(self, params, name='Desorption')

def possible(self, system, site, other_site):

if system.sites[site].covered == 1 and system.sites[other_site].covered == 1:
return True
else:
return False

def get_rate(self, system, site, other_site):
R = 1.
return R

def do_event(self, system, site, other_site):
# empty the sites:
system.sites[site].covered = 0
system.sites[other_site].covered = 0

def get_involve_other(self):
return True


Now we will store references to the classes in a list:

events = [Adsorption, Desorption]


How to accelerate the kMC simulations is described in the tutorial on Accelerating kMC and not used in this simplest case.

Step 2. Define sites

In this step, the sites are defined from an ASE.Atoms object. We create one site for each atom in a 10x10 fcc(111) surface, all with the same site-type stype=0 and without any covering species covered=0:

from ase.build import fcc111
from base.sites import SiteBase

a0 = 4.00  # Lattice Parameter (not related to DFT!)
atoms = fcc111("Pt", size=(10, 10, 1), a=a0)
sites = []
# Define a site for each atom that is empty with no pre-defined neighbors:
for i in range(len(atoms)):
sites.append(SiteBase(stype=0, covered=0, ind=i))


Now we have a list of empty sites, which are used to instantiate a system.

Step 3. Instantiate system and neighborlists

Here, the system is created and the sites are connected by calculating a neighborlist. In this example we assign neighbors within one nearest-neighbor distance:

import numpy as np
from base.system import SystemBase

p = SystemBase(atoms=atoms, sites=sites)
Ncutoff = a0 / np.sqrt(2.) + 0.05  # Nearest neighbor cutoff

for i, s in enumerate(sites):
for j, sother in enumerate(sites):
dcur = atoms.get_distance(s.ind, sother.ind, mic=True)
if dcur < Ncutoff and j != i:
s.neighbors.append(j)


mic=True uses periodic boundary conditions to imply an infinite surface.

Step 4. Instantiate a NeighborKMC object and run

Now we are ready to instantiate a NeighborKMC.NeighborKMCBase object, which is connecting the ingredients created in the previous step. The main part of the kinetic Monte Carlo procedure is in the NeighborKMC.NeighborKMCBase, but some details and logging should be defined by the user. That is done here in the class: simple_NKMC.

from base.logging import Log
from base.kmc import NeighborKMCBase

class simple_NKMC(NeighborKMCBase):

# First we initialize the kMC simulation and load the parameters
def __init__(self, system, tend, parameters={}, events=[]):
self.events = [ev(parameters) for ev in events]
NeighborKMCBase.__init__(self, system=system,
tend=tend, parameters=parameters)

# We also define a run_kmc, which runs the actual kMC simulation. Also the logging is defined here.
def run_kmc(self):
logparams = {}
logparams.update(self.parameters)
logparams.update({"tend": self.tend,
"Nsites": self.system.Nsites,
"Number of events": len(self.events),
"Number of site-types (stypes)": len(list(set([m.stype for m in self. system.sites])))
})
log = Log(logparams)

stepN_CNT = 0  # Parameter to count LogSteps threshold
stepNMC = 0    # Parameter to count the number of executed kMC steps

while self.t < self.tend:
self.frm_step()         # Execute a kMC step

if stepN_CNT >= self.LogSteps:       # Only for Logging purposes
print("Time : ", self.t, "\t Covs :", self.system.get_coverages(self.Nspecies))
log.dump_point(stepNMC, self.t, self.evs_exec)
stepN_CNT = 0

stepN_CNT += 1
stepNMC += 1


So now after defining the simple_NKMC object, the parameters can be loaded:

parameters = {"pA": 10., "Name": "Quickstart simulation"}
sim = simple_NKMC(system=p,
tend=10.0, # end after 10.s.
parameters=parameters, # parameters for event rate-constants.
events=events) # the list of events


And finally, we can now run the simulation by invoking:

sim.run_kmc()


Then it is just to have a cup of coffee and wait.

Afterthoughts

While this example shows how simple it can be to run a simulation, in all following tutorials and examples the details of the simulations are stored in so-called user-files:

• user_kmc.py can be used to customize the kMC routine, especially run_kmc().
• user_events.py can be used to store the event-types.
• user_energy.py can be used to store functions for obtaining energies used to calculate event rate constants.
• user_entropy.py can be used to store entropy calculation functions.
• user_constants.py can be used to store global and physical constants.

### Adsorption of atomic species A on surface¶

In this tutorial we show how to set-up a MonteCoffee simulation, requiring only a view steps. The entire files needed for this tutorial are in test.py and the references to the other modules mentioned herein.

In this tutorial the simple adsorption/desorption of an atomic species on a plain surface is demonstrated and the results compared to the solution of a mean-field model. The reaction to simulate is:

$A + * \longleftrightarrow A^*$

and the time evolution of the coverage of species A according to the mean-field model (see e.g.: Fichthorn and Weinberg)

$\theta(t) = \frac{r_A}{r_A+r_D}(1-e^{-(r_A+r_D)*t})$

with $$\theta$$ the coverage of species A, $$t$$ the time and $$r_{A,D}$$ the rate of adsorption and desorption respectively.

Before the simulation, constants, reaction sites and system as well as the events have to defined which will be shown in the next steps.

First, various parameters must be set and stored in a parameter-dictionary:

tend = 10.  # End time of simulation (s)
latt_param = 4.00  # Lattice Parameter (not related directly to DFT (but could be) )
Ncutoff = a / np.sqrt(2.) + 0.05  # Nearest neighbor cutoff


#### Define sites and system¶

One site is defined for each surface atom using an Ase.Atoms object. We start with an empty list of sites and a fcc(100) surface in ASE. That a Pt surface is chosen as example has no special meaning and for this example any kind of atom could be chosen:

from ase.build import fcc(100)
from user_sites import Site

surface = fcc100("Pt", a=latt_param, size=(10,10,1))
sites = []


Now we can create a site, free of adsorbates, for each surface atom with a stype . In case of the (100) surface, all surface atoms are of the same type:

# Create a site for each surface-atom:
for i in range(len(atoms)):
sites.append(Site(stype=0,
covered=0, ind=i))


Here, the block ind=i stores the index of the atom in the ASE.Atoms object on the NeighborKMC.user_sites.Site object.

Finally, we need to define neighbor lists. It is simplest to define this according to the nearest neighbor distances:

# Set the neighbor list for each site using distances.

for i, s in enumerate(sites):
for j, sother in enumerate(sites):
# Length of distance vector:
dcur = self.atoms.get_distance(s.ind, sother.ind, mic = pbc)

# If the site is a neighbor:
if dcur < Ncutoff and j != i:
s.neighbors.append(j)


Now the NeighborKMC.user_system.System object can be defined from the collection of sites:

from user_system import System
p = System(atoms=atoms, # store ASE.Atoms as well
sites=sites)


#### Define events¶

Various event-types are defined, which are stored in user_events.py. For each possible type of event, a class is derived from NeighborKMC.base.events.EventBase. In this case, we need to define two different events, the adsorption of species A, and correspondingly the desorption.

First we import the necessary functions, classes, and constants:

from base.events import EventBase


Now we derive a class to contain the event:

class AAdsEvent(EventBase):
def __init__(self, params):
EventBase.__init__(self, params)


The constructor __init__(self,params) attaches relevant parameters to the object. We need a function possible(self,system, site, other_site) that returns True if the event is possible on the current site-pair. For single atom adsorption it does not matter if the other_site is covered or not. Thus we are only interested in the site itself.

def possible(self, system, site, other_site):
# If site is uncovered
if (system.sites[site].covered == 0):
return True
else:
return False


Thus, for the event to be possible, the site needs to be empty. Now we also need to define a function get_rate(self, system, i_site, other_site) that returns the rate constant. To keep this as simple as possible, the rate constant is chosen to be $$R=1$$.

def get_rate(self, system, i_site, other_site):
R = 1.
return R


Each event requires a method do_event(self,system, site, other_site) to perform modifications to the site-occupations when fired:

def do_event(self, system, site, other_site):
system.sites[site].covered = 1


In this case, upon adsorption the site is covered with the species A, represented by the number 1 within the code.

To take care of the correct time evolution in MonteCoffee we introduce an additional block which returns if either neighboring sites are involved or not. Here no neighboring sites are involved, thus we return False.

def get_involve_other(self):
return False


Finally, the events are stored in the main simulation file, in a list:

events = [AAdsEvent, ADesEvent]


The numbering of events is determined by the order in the list events defined here and the output is numbered accordingly.

#### Define and run simulation¶

Now the simulation object NeighborKMC.user_kmc.NeighborKMC can be defined and the simulation performed:

parameters = { "Name": "A ads/des Simulation",
"reverses ": reverse_events}

# Instantiate simulator object.
sim = NeighborKMC(system=p, tend=tend,
parameters=parameters,
events=events,
rev_events=reverse_events)
result = sim.run_kmc()
print("Simulation end time reached ! ! !")


#### Analyze results¶

The results are analyzed by reading in the code output. Here, we would like to calculate the A coverage as a function of time for the entire system:

import numpy as np
Nsites = float(len(covs))
cov_A = [sum([1 for val in covs[i] if val == 1]) / Nsites for i in range(len(covs))]


This can be plotted as done in the following example with matplotlib

import matplotlib.pyplot as plt
plt.plot(time, cov_A, '-k')
plt.xlabel("Time [s]")
plt.ylabel("Coverage")
plt.savefig('coverage_spec_A.pdf')


To compare the effect of the used simulation surface on the result and also compare to the result of the mean-field model in the following a plot is shown with surface sizes of (5x5), (10x10) and (100x10) corresponding to 25, 100 and 1000 surface sites respectively. If an increase in the number of sites is not possible, it is recommended that multiple identically prepared simulations are performed. (see example on Parallel simulations and calculating turnover frequencies).

This tutorials moves on from the very basic single atom adsorption to the dissociative adsorption of a diatomic molecule of type B2:

$B_2 + 2^* \longleftrightarrow B^* + B^*$

With respect to the case of the single atom adsorption the main changes are done in the definition of the events. They are presented together with some comparison to a mean-field model. The entire tutorial is shown in test.py and the references to the other modules mentioned therein.

#### Define events¶

As before, defined event-types are stored in user_events.py. For each possible type of event, a class is derived from NeighborKMC.base.events.EventBase. In this case, we again need to define two different events, the adsorption of species B2, and correspondingly the desorption.

First we import the necessary functions, classes, and constants:

from base.events import EventBase


Now we derive a class to contain the event:

class B2AdsEvent(EventBase):
def __init__(self, params):
EventBase.__init__(self, params)


The constructor __init__(self,params) attaches relevant parameters to the object. We need a function possible(self,system, site, other_site) that returns True if the event is possible on the current site-pair. For the dessoziate adsorption, both neighboring sites have to be empty. Thus now the other_site becomes important.

def possible(self, system, site, other_site):
# If site is uncovered
if (system.sites[site].covered == 0 and system.sites[other_site] == 0):
return True
else:
return False


Now we also need to define a function get_rate(self, system, i_site, other_site) that returns the rate constant. To keep this as simple as possible, the rate constant is chosen to be $$R=1$$.

def get_rate(self, system, i_site, other_site):
R = 1.
return R


Each event requires a method do_event(self,system, site, other_site) to perform modifications to the site-occupations when fired:

def do_event(self, system, site, other_site):
system.sites[site].covered = 1
system.sites[other_site].covered = 1


In this case, upon adsorption the site and also the other_site is covered with the species B, represented by the number 1 within the code.

To take care of the correct time evolution in MonteCoffee we introduce an additional block which returns if either neighboring sites are involved or not. Here the neighboring sites are involved, thus we return True.

def get_involve_other(self):
return True


Finally, the events are stored in the main simulation file, in a list:

events = [B2AdsEvent, B2DesEvent]


Thus to run a kinetic Monte Carlo simulation of dissoiative adsorption, only the user_event.py file has to be changed with respect to the single atom adsorption, and the imported events updated in test.py.

#### Analyze results¶

To compare with the mean-field model we solve the following coupled differential equations for the surface coverages $${\theta_i}$$:

$\begin{split}\frac{d\theta_B}{dt} & = k^{+}\theta_*^2 - k^-\theta_B^2 \\ \theta_* & = 1 - \theta_B\end{split}$

with $$k^{+,-}$$, being the rate of the forth and back reaction respectively. Comparing the mean-field results with kinetic Monte Carlo simulations is only in this very simple cases, which do not include any adsorbate-adsorbate interactions or diffusion limitations possible. Also one has to account in the mean-field model for the coordination number of the surface site over which the reaction takes place. Using the (100) surface, we have 4 possible pairs of neighbouring sites at which the adsorption can happen. In consequence, $$k^{+,-}$$ has to be multiplied by 4. In the following image, the time evolution for both models is shown for various system sizes in the case of the kinetic Monte Carlo simulation. As for the single atom adsorption, both models agree and an increase in surface size reduces the variations of the kinetic Monte Carlo simulation.

### Adsorption and reaction of species A and B2¶

Turning from simple adsorption events towards surface reactions, we demonstrate here the $$A+B_2$$ reaction on a fcc(100) surface and a truncated octahedral nanoparticle. The reaction is modeled using the following reactions:

$\begin{split}A(g) + ^* \longleftrightarrow A^* \\ A^* + ^* \longleftrightarrow * + A^* \\ B_2(g) + 2^* \longleftrightarrow 2B^* \\ B^* + ^* \longleftrightarrow ^*+B^* \\ A^*+B^* \longrightarrow AB(g) + 2^*\end{split}$

Thus additional to the already described events of single atom and dissociative adsorption, the diffusion of species $$A$$ and $$B$$, as well as the formation and desorption of $$AB$$ has to be implemented. In the following we present first the event classes which are additionally needed and are the same for the nanoparticle and the surface.

#### Define events¶

The defined events are stored in user_events.py. For each possible type of event, a class is derived from NeighborKMC.base.events.EventBase. In this case, we need to define seven different events, the adsorption, desorption and diffusion of species A and B2, and the reaction between both species. The adsorption and diffusion of A and B2 are defined as before. For the A diffusion we derive a class to contain the event:

class ADiffEvent(EventBase):
def __init__(self, params):
EventBase.__init__(self, params)


The constructor __init__(self,params) attaches relevant parameters to the object. We need a function possible(self,system, site, other_site) that returns True if the event is possible on the current site-pair. Here we assume that the species A is assigned to the coverage number:2. A diffusion event can take place if the site is covered with A and the neighbor site empty or vice versa.

def possible(self, system, site, other_site):
# If site is uncovered
if (system.sites[site].covered == 2 and system.sites[other_site].covered == 0) or (system.sites[site].covered == 0 and system.sites[other_site].covered == 2):
return True
else:
return False


Now we also need to define a function get_rate(self, system, i_site, other_site) that returns the rate constant. To keep this as simple as possible, the rate constant is chosen to be $$R=1$$.

def get_rate(self, system, i_site, other_site):
R = 1.
return R


Each event requires a method do_event(self,system, site, other_site) to perform modifications to the site-occupations when fired:

def do_event(self, system, site, other_site):
old_cov_site = system.sites[site].covered
old_cov_other_site = system.sites[other_site].covered
system.sites[site].covered = old_cov_other_site
system.sites[other_site].covered = old_cov_site


In this case, the coverage of the site and other_site are interchanged.

To take care of the correct time evolution of the MonteCoffee we introduce an additional block which returns if either neighboring sites are involved or not. Here no neighboring sites are involved, thus we return False.

def get_involve_other(self):
return True


The diffusion of species B is exactly as for species A implemented, with the difference that B is represented by 1 in the code.

For the reaction, A and B have to be present on either site or other_site. So the reaction becomes possible if:

def possible(self, system, site, other_site):
if (system.sites[site].covered == 1 and system.sites[other_site].covered == 2) or (system.sites[site].covered == 2 and system.sites[other_site].covered == 1):
return True
else:
return False


The rate constant can be chosen according to:

def get_rate(self, system, site, other_site):
R = .1
return R


and after the formation of AB, both sites: code:site and code:other_site are emptied:

def do_event(self, system, site, other_site):
system.sites[site].covered = 0
system.sites[other_site].covered = 0


In this example, it is assumed that the desorption of the formed product AB is instantaneous without the possibility to re-absorb and split into A and B.

#### Reaction over a (100) surface¶

For the reaction over a (100) surface we use ASE to define the sites and their neighbors:

from ase.build import fcc(100)
from user_sites import Site

surface = fcc100("Pt", a=latt_param, size=(20,20,1))
sites = []

# Create a site for each surface-atom:
for i in range(len(atoms)):
sites.append(Site(stype=0,
covered=0, ind=i))


A 20x20 lattice is used to ensure convergence of the coverage and each atom site is connected to its four neighbors. The in Fig. 1 shown coverage is based on the following rate constants: Adsorption, desorption and diffusion rate: 1 s-1 and the reaction rate: 0.1 s-1. It should be noted the mean-field rate of 1 s-1 for B2 adsorption/desorption corresponds to 0.5 s-1 in the kMC model. With the reaction being the rate limiting step, all coverages are 1/3 after an equilibration period as can be seen in Fig. 1. Time evolution of the coverage of species A and B on a 20x20 (100) surface.

#### Reaction over a nanoparticle¶

Similar to the (100) surface, we employ again ASE to define the sites, this time constructing a truncated octahedron consisting of 260 atoms of which 144 are exposed on the surface. In contrast to the surface, we are going to use two different type of sites, stype=0 representing the (111) facet sites and stype=1 representing the (100) facet sites, edges and corners. In the inset of Fig. 2 the different types of sites are visualized.

The atoms are defined as follows:

atoms = Octahedron("Pt", 8, cutoff=3, latticeconstant = a)
sites = []
write('trunc_octa.traj', atoms) # see how the nanoparticle looks like


but to assign the site types to the surface atoms, we use the coordination number to distinguish them, the (111) facet sites having a coordination number of 9. The coordination number of each atom is calculated and a list, surface_atom_ids created in which the surface atom ids’ (coordination number < 12) are stored.

 CNS = np.zeros(len(atoms))
for i, at in enumerate(atoms):
pcur = at.position
dp = np.sqrt([(p - pcur) ** 2. + (p - pcur) ** 2. +
(p - pcur) ** 2. for p in atoms.positions])
CNS[i] = len([val for val in dp if 0. < val < a / np.sqrt(2) + 0.01])
surface_atom_ids = [i for i in range(len(CNS)) if CNS[i] < 12]


In difference to the (100) surface, each site has now stype assigned which is either 0 or 1. For the (100) surface all sites have the same stype=0.

for i,indic in enumerate(surface_atom_ids):
if CNS[indic] == 9:
sstype = 0
else:
sstype = 1
sites.append(Site(stype=sstype, covered=0, ind=indic))


The neighbor list is calculated for the surface atom shell only (the atoms saved in the sites-list). All atoms for the nanoparticle do not have the same number of direct neighbors.

In the following we keep all rates fixed to the (100) surface ones, beside the rates for the A adsorption. For the A adsorption we are going to employ different, site-dependent adsorption rates. Therefore, beside ensuring that the site is empty, also the stype has to be considered to determine the adsorption rate. That is done in the following way in the AAdsEvent :

def get_rate(self, system, site, other_site):
if system.sites[site].stype == 0:
R = 1.
elif system.sites[site].stype == 1:
R = 10.
return  R


To see the effect of the rate of the A adsorption on the turn over frequency (TOF) of the simulation, we study four different combinations: First using either a rate constant of 1 s-1 or 10 s-1 on both sites and second by using the mixed cases, having 1 s-1 for stype=0 and 10 s-1 for stype=1 or vice versa. The results can be seen in Fig. 2. In the case of employing the same rate for the A adsorption as for the B adsorption the TOF is the highest, and with having a 10 times faster A adsorption than B adsorption, it being the lowest. In the case of high A adsorption, the sites are blocked leading to poisoning. For the mixed cases, the TOF is higher for the one with rate stype=0: 10 s-1. Not as many sites with stype=0 exist and therefore the A poisoning is less pronounce. Turn over frequency for different choices of A adsorption rates. The first number refers to the rate in s-1 of the site with stype=0 and the second number to the rate of stype=1. The inset shows the different type of sites on a truncated octahedron with orange being: stype=0 and blue stype=1.

### CO oxidation on Pt(111)¶

We present here the CO oxidation on Pt(111) which is published along with the oxidation on Pt nanoparticles in: M. Jørgensen and H. Grönbeck, ACS Catal., 7, 5054-5061 (2017) .

For the CO oxidation the following five chemical reactions have to be considered:

$\begin{split}O_2 + 2^* & \longleftrightarrow 2O^* \\ CO + * & \longleftrightarrow CO* \\ CO^* + O^* & \longrightarrow CO_2 \\ CO^* + * & \longleftrightarrow * + CO^* \\ O^* + * & \longleftrightarrow * + O^*\end{split}$

which are dissociative adsorption of oxygen, adsorption of CO, reaction of adsorbed O and CO to CO2 and CO and O diffusion. In the model it is assumed that the reaction to CO2 is not reversible and the reaction product is directly desorbing.

In this example, species 0 denotes empty sites, 1 is CO and 2 refers to O. All energies were obtained using density functional theory and are given in the paper mentioned above.

#### Define constants and parameters¶

Any global constants can be defined in user_constants.py. In this example this file stores physical constants but also the calculated vibrations of the various adsorbed species and of the CO2 formation transition state and their physical properties.

Next, the simulation parameters must be set and stored in a dictionary:

T = 800.  # Temperature (K)
pCO = 2E3  # CO pressure (Pa)
pO2 = 1E3  # O2 pressure (Pa)
tend = 1E-3  # End time of simulation (s)


Here the Temperature is set to 800 K, and the pressure ratio of CO to O2 to 2:1. The pressure is given in Pa. Also the end time is defined, but because we are interested mainly in steady-state properties, the simulation is run until steady-state is reached and usually stopped before the given end time reached.

#### Define sites and system¶

One site is defined for each surface atom using an Ase.Atoms object. We are using the ase.build class to construct our (111) surface. It consists of only one layer, because its only purpose is to simplify the site-connectivity set-up. Therefore, the used lattice-parameter is also not related to any results of density functional theory calculations or experiments. Thus starting with an empty list of sites we construct a $$10x10$$ surface consisting of 100 atoms:

from ase.build import fcc111
from user_sites import Site

sites = []
a = 4.00  # Lattice Parameter (not related to DFT!)
Ncutoff = a / np.sqrt(2.) + 0.05  # Nearest neighbor cutoff
atoms = fcc111("Pt", a = a, size = (10,10,1))
atoms.write('surface.traj')


We also have written out the prepared surface as .traj file to see if it looks as intended. As next step we have to add a site for each surface atom. Although we know that the oxygen preferably adsorbs on fcc sites and the CO on top sites, we assume here all sites to be equal (coarse-graining) but use the corresponding energies of the preferred sites.

# Create a site for each surface-atom:
for i in range(len(atoms)):
sites.append(Site(stype=0,
covered=0, ind=i))


Here, the block ind=i stores the index of the atom in the ASE.Atoms object on the NeighborKMC.user_sites.Site object, which can be later on used to write out .traj-files during the simulation.

Now the NeighborKMC.user_system.System object can be defined from the collection of sites:

from user_system import System
p = System(atoms=atoms, # store ASE.Atoms as well
sites=sites)


In addition to the overall system, we also need to define neighbor lists. It is simplest to define this according to the nearest neighbor distances. Thus we call in the main definition file of our kMC simulation:

p.set_neighbors(Ncutoff,pbc = True)


which sets the global neighbor list based on distances for our System. Because we are using a surface, periodic boundary conditions are desirable. The actual function is then defined in NeighborKMC.user_sites.Site and looks as follows:

def set_neighbors(self, Ncutoff, pbc=False):

if self.atoms is None:
raise Warning("Tried to set neighbor-distances in user_system.set_neighbors() with self.atom = None")

for i, s in enumerate(self.sites):  #For all sites
for j, sother in enumerate(self.sites): #Check all the other sites
dcur = self.atoms.get_distance(s.ind, sother.ind, mic=pbc) # use ase function get_distance
if dcur < Ncutoff and j != i:
s.neighbors.append(j)        #add site j to neighbor list of site i
if len(self.neighbors) == 0:             #check if neighbors exists -- otherwise the site will not interact with each other
self.neighbors = [s.neighbors for s in self.sites]
self.verify_nlist()


#### Define reaction energies and entropies¶

In this step, the reaction energies, or methods to calculate these, are defined in user_energy.py and the entropies in user_entropy.py. That makes it simple to do all the book-keeping accordingly.

We used the from density functional theory obtained energies for the adsorption of oxygen on the fcc site and of CO on the top site and diffusion constants for oxygen to diffuse from fcc to fcc site and for CO for from top to top site. In this example, the adsorption energies are defined to be positiv in contradiction to the more generally used negative notation in theoretical papers.

EadsCO = 1.36

EdiffCO = 0.08 # CO diffusion barrier
EdiffO = 0.58 # O diffusion barrier


From the adsorption energies the activation energy for the CO2 formation is calculated from the BEP relations according to:

def get_Ea(ECO, EO):
ETS = 0.824 * (-EO -ECO) + 0.168 + 0.47238  # How much larger is the energy of CO and O wrt Pt(111)
Ea = ETS + ECO + EO  # Translate the barriers relative to Pt(111)
return Ea


The reason why not one single activation energy is used, are the repulsive adsorbate-adsorbate interactions which depend locally on the coverage of the neighbors of the reaction site:

def get_repulsion(cov_self, cov_NN, stype):
repulsion = 0.
ECOCO = 0.19  #  How CO affects CO
EOO = 0.32  # How O affects O
ECOO = 0.3  # How CO affects O
EOCO = 0.3  # How O affects CO
HInttwo = [[0., 0., 0.], [0., ECOCO, EOCO],
[0., ECOO, EOO]]  # Two body interaction Hamiltonian 3x3 because 0 = empty.
for j in cov_NN:  # For each covered Neighbor, give a repulsion:
repulsion += HInttwo[cov_self][j]

return repulsion


The in user_entropy.py stored entropies are included in the calculation of the rates as will be shown in the next section. Generally we define here the entropy for gas-phase CO and oxygen, as well as a method to calculate harmonic adsorbate entropy. The definitions of this functions can be looked up in this file directly.

#### Define events¶

Here event-types are defined, which are stored in user_events.py. For each possible type of event, a class is derived from NeighborKMC.base.events.EventBase. For example take the $$\mathrm{CO_2}$$ formation. This event is defined in user_events.py as follows.

First we import the necessary functions, classes, and constants:

from base.events import EventBase
from user_entropy import get_Zvib, get_Z_CO, get_Z_O2
from user_constants import mCO, mO2, Asite, modes_COads, modes_Oads,\
modes_TS_COOx, modes_COgas, modes_O2gas, kB, eV2J, s0CO, s0O, h,
from user_energy import EadsCO, EadsO, get_Ea, get_repulsion, EdiffCO, EdiffO


Now we derive a class to contain the event:

class COOxEvent(EventBase):
def __init__(self, params):
Zts = get_Zvib(params["T"], modes_TS_COOx)
self.Zratio = Zts / Zads
EventBase.__init__(self, params)


The constructor __init__(self,params) attaches relevant parameters to the object, and self.Zratio is the ratio between the partition functions in the initial state and transition state, used to calculate the rate constant in transition state theory. We need a function possible(self,system, site, other_site) that returns True if the event is possible on the current site-pair:

def possible(self, system, site, other_site):
if (system.sites[site].covered == 1 and system.sites[other_site].covered == 2) or \
(system.sites[site].covered == 2 and system.sites[other_site].covered == 1):
return True
else:
return False


Thus, for the event to be possible, the site needs to be covered by 1 (CO) and the neighbor site by 2 (O) or the other way round. That is originated by the use of single neighbor site pairs. Thus a pair with the indexes (10,11) would be the same as (11,10) in the code to avoid double counting in the time list.

Next we need to define a function get_rate(self, system, i_site, other_site) that returns the rate constant:

def get_rate(self, system, i_site, other_site):
if system.sites[site].covered == 1:
Ncovs_CO = [system.sites[n].covered for n in system.neighbors[site] ]
Ncovs_O = [system.sites[n].covered for n in system.neighbors[other_site]]
else:
Ncovs_CO = [system.sites[n].covered for n in system.neighbors[other_site] ]
Ncovs_O = [system.sites[n].covered for n in system.neighbors[site]]
ECO -= get_repulsion(1, Ncovs_CO, 0)
EO -= get_repulsion(2, Ncovs_O, 0)
Ea = max(0., get_Ea(ECO, EO)) # No negative energy barriers
return self.alpha * self.Zratio * np.exp(-Ea/(kB * self.params['T'])) * kB * self.params['T'] / h


The adsorption energies are fixed, but weaken in the case of any repulsions. A call is made to get_Ea(ECO, EO) to obtain the reaction energy barrier. The rate constant is multiplied with self.alpha, because self.alpha is the slowing-down factor that is adjusted dynamically for each event during simulation.

Also the event requires a method do_event(self,system, site, other_site) to perform modifications to the site-occupations when fired:

def do_event(self, system, site, other_site):
system.sites[site].covered = 0
system.sites[other_site].covered = 0


In this case, the two sites containing CO and O are simply emptied. At the end we define the method get_involve_other(self) to define if the neighboring sites are actually involved in the event:

def get_involve_other(self):
return True


After giving this example for one event, the other events can be defined similarly. How to define single site and dissociative adsorption was shown before. Only the rates have to be adjusted according to transition state theory. Having all events for each type of reaction defined in this order:

1. NeighborKMC.user_events.COAdsEvent for CO adsorption.
1. NeighborKMC.user_events.CODesEvent for CO desorption.
1. NeighborKMC.user_events.OAdsEvent for O2 dissociative adsorption.
1. NeighborKMC.user_events.ODesEvent for O2 desorption.
1. NeighborKMC.user_events.CODiffEvent for CO diffusion.
1. NeighborKMC.user_events.ODiffEvent for O diffusion.
1. NeighborKMC.user_events.COOxEvent for CO+O -> CO2.

we can now store the event-class references in a list for the NeighborKMC simulations and define accordingly a list of reverse events. In this example the CO (O) adsorption and desorption are reverse to each other. In addition are the diffusion processes reverse to themselves. The CO oxidation itself is not reversible. Thus it isn’t defined in this list.

events = [COAdsEvent, CODesEvent, OAdsEvent,
ODesEvent, CODiffEvent,
ODiffEvent, COOxEvent]
reverse_events = {0: 1, 2: 3, 4: 4, 5: 5}


The numbering of events is determined by the order in the list events defined here.

#### Define and run simulation¶

Now the simulation object NeighborKMC.user_kmc.NeighborKMC can be defined and the simulation performed:

# Instantiate simulator object.
sim = NeighborKMC(system=p, tend=tend,
parameters=parameters,
events=events,
rev_events=reverse_events)
sim.run_kmc()
print("Simulation end time reached ! ! !")


#### Analyze results¶

The results are analyzed by reading in the code output. For example, if we need to calculate the CO and O coverage as a function of time for the entire system:

import numpy as np
Nsites = float(len(covs))
cov_CO = [sum([1 for val in covs[i] if val == 1]) / Nsites for i in range(len(covs))]
cov_O = [sum([1 for val in covs[i] if val == 2]) / Nsites for i in range(len(covs))]
cov_free = [sum([1 for val in covs[i] if val == 0]) / Nsites for i in range(len(covs))]


Typically, a turnover frequency is also relevant to calculate:

evs_exec = np.loadtxt("evs_exec.txt")
TOF = evs_exec[-1] / (Nsites*time[-1]) # How many CO+O->CO2 has fired per time and site.


Often it can be useful to discard points out of steady-state. The detailed evolution of the time with the sites can be found in detail_site_event_evol.hdf5. To draw statistically sound conclusions, it is recommended that multiple identically prepared simulations are performed and in this case at least 100 CO oxidation events are performed in the steady state.

### Parallel simulations¶

Kinetic MonteCarlo simulations are stochastic in nature, making it reasonable to perform multiple identically prepared simulations to assess the convergence. In the following, we assume that ASE and MPI4PY are installed.

To submit a simulation, defined in a file kmc_master_parallel.py, in an environment that implements SLURM, one can submit a bash script as

#!/usr/bin/env bash
#SBATCH -N 1
#SBATCH -n 10
#SBATCH -p PROJNAME
#SBATCH -A ACCNO
#SBATCH -o out.txt
#SBATCH -e err.txt
#SBATCH -t 96:15:00
#SBATCH -J project_dir/submit_dir
#SBATCH --mail-user=USER@UNI.se
#SBATCH --mail-type=END

cp *.py $TMPDIR cp kMC_options.cfg$TMPDIR
cp -r base $TMPDIR cd$TMPDIR

while sleep 1800; do
# This will be executed once per every 3600 seconds
rsync -a $TMPDIR/*$SLURM_SUBMIT_DIR
done &     # The &-sign after the done-keyword places
# the while-loop in a sub-shell in the background
LOOPPID=$! # Save the PID of the subshell running the loop mpirun -np 10 python kmc_master_parallel.py # Copy the files back after run: rsync -a$TMPDIR/* $SLURM_SUBMIT_DIR kill$LOOPPID


This script copies all python files in MonteCoffee to a compute node, and cds into the simulation directory (\$TMPDIR) on the node. Then a while loop copies all dirs called run_* back to the submission directory every half hour. The script named kmc_master_parallel.py is then executed with mpirun. The script kmc_master_parallel.py can at first be very similar to the quick-start example. For simplification the detailed definitions as in quick-start are not given here, so that one can see the differences. Mainly they are in the first part of the code.

import os
from ase.build import fcc111
from ase.parallel import MPI4PY
from base.sites import SiteBase
import numpy as np
from base.system import SystemBase
from base.logging import Log
from base.kmc import NeighborKMCBase

world = MPI4PY()
rank = world.rank # What number simulation copy am I?
size = world.size # How many total simulation copies?
rundir = "run_"+str(rank)  # Name of the dir I create

os.mkdir(rundir) # Create dir
os.system("cp kMC_options.cfg "+rundir) # copy kMC_options.cfg here
os.chdir(rundir)

a0 = 4.00  # Lattice Parameter (not related to DFT!)
Ncutoff = a0 / np.sqrt(2.) + 0.05  # Nearest neighbor cutoff

atoms = fcc111("Pt", size=(10, 10, 1), a=a0)
sites = []

for i, indic in enumerate(atoms):
sites.append(Site(stype=0, covered=0, ind=[i]))

events = [Adsorption, Desorption]

p = SystemBase(atoms=atoms, sites=sites)
p.set_neighbors(Ncutoff)

parameters = {"pA": 10., "Name": "Parallel Simulation"}

sim = simple_NKMC(system=p,
tend=10.,
parameters=parameters,
events=events)

sim.run_kmc()


For further explanations about using MPI4PY within ASE, please see the ASE documentation on parallel calculations.

In general, it can be useful to assign a large tend and let the bash-script runtime determine the end of simulation. Because the code itself writes out log-files regularly, one will not loose any informations by letting the script runtime determine the end of the simulation.

### Special site-rules¶

#### Custom site-attributes¶

Sites in MonteCoffee, per default only contain a variable that determine their type named stype. Stype is used to analyze the rates and coverages over different sites in a system. However, to calculate reaction energies, it can be good to attach a coordination number to the class as well. This can simply be done by adding a parameter coordination_number to the constructor in user_sites.py as

class Site(SiteBase):

def __init__(self, stype=0, covered=0, ind=None, lattice_pos=None, coordination_number=None):
SiteBase.__init__(self, stype=stype, covered=covered, ind=ind,
lattice_pos=lattice_pos)
self.cn = coordination_number


Then the get_rate() methods of user_events.py can access the coordination number as

class A(EventBase):
...
def get_rate(self, system, site, other_site):
return self.alpha*1000.*system.sites[site].cn


#### Using stype for everything¶

stype that belongs to a NeighborKMC.user_sites.Site object is useful to define special rules for performing events. For a binary alloy system with 10 different generalized coordination numbers we have 20 different types of sites, thus stype can take on the values from 0 to 19.

To use stype, let us assume that we have defined

from user_sites import Site
from user_system
s1 = Site(stype = 0, covered = 0, ind = 0)
s2 = Site(stype = 1, covered = 0, ind = 1)


One example of a special rule is to make events that are only possible if stype == 1:

class A(EventBase):
...

def possible(self, system, site, other_site):
if system.sites[site].stype == 1:
return True
else:
return False


Another special rule is to return a different rate-constant based on stype

class A(EventBase):
...

def get_rate(self, system, site, other_site):
R = 1000*system.sites[site].stype+0.1
return self.alpha*R


This can be useful when having multiple different sites on a nanoparticle. If we want to calculate the rate-constant based on transition state theory, a different reaction energy barrier can be defined for each site’s and neighbor-site’s type as

from user_constants import kB
import numpy as np
class A(EventBase):
...
def get_rate(self, system, site, other_site):
stype = system.sites[site].stype
stype_other = system.sites[other_site].stype
stype_avg = 0.5*(stype+stype_other)
Ea = 1.08 + (4-stype_avg)*0.1
return self.alpha*1E12*np.exp(-Ea/(kB*self.params["T"]))


Here, a transition state like rate constant is returned, with a pre-exponential factor of $$\dfrac{k_\mathrm{B}T}{h}\dfrac{Z^{ts}}{Z^{ini}} \approx 10^{12}\,\mathrm{s}^{-1}$$ and the energy barrier is based on the average site-type number of the two sites.

## Special features¶

We present in this section different special features needed to treat real systems. We show how the time acceleration works and how to converge the different parameters and present results obtained with the time acceleration. In addition, the importance of adsorbate-adsorbate interactions are pointed out and how turn-over frequencies are calculated.

### Accelerating kMC¶

In MonteCoffee three different time acceleration schemes are implemented. All of them are based on the super-basin hopping concept, e.g. explained here: Generalized temporal acceleration scheme (Dybeck et al). One can picture this form of time acceleration in the following way: While executing reactions a super-basin $$S$$ is explored. The already explored region is called $$S_A$$ and the region to be explored, $$S_B$$. The explored area is extended while executing reversible reactions. The whole super-basin $$S$$ can be left upon executing a non reversible reaction, belonging to $$S_N$$. The scaling is based on the average time needed to exit the current super-basin $$S$$. The following picture visualizes the described super-basin ‘hopping’. On the left side, the super-basin $$S$$ is shown and how it is exited, entering a new super-basin $$S'$$ on the right side, by an event belonging to $$S_N$$. The most basic is the scaling of equilibrated reactions $$R_\mathrm{eq}$$ with a constant factor $$N_\mathrm{f}$$:

$R_\mathrm{eq} = \frac{1}{N_\mathrm{f}} \cdot R.$

The other two methods are based on the rate constants or the rate respectively. The scaling of the rate is adapted from the Generalized temporal acceleration scheme by Dybeck et al and the rate constant based scaling follows the same principle. In general, the code determines which events are fast and thus likely quasi-equilibrated, and slows down these events periodically. This is done by comparing the rates (rate constants) of frequently executed quasi-equilibrated events with infrequently executed quasi-equilibrated events and non-equilibrated events. In that manner, the fastest events are continuously slowed down until a non-equilibrated event proceeds upon which the rate constants are unscaled again.

In detail, after a scaling period of $$N_s$$ steps, the algorithm determines if any new reactions have become quasi-equilibrated by comparing the number of times the event has proceeded forward and backward the last $$n_e$$ simulation steps. A reaction is deemed quasi-equilibrated if it has been executed at least $$n_e$$ times, and fulfills:

\begin{equation} \dfrac{|n_f - n_b|}{n_f + n_b} < \delta \end{equation}

where $$\delta \in [0,1]$$ is a tolerance for determining if the event is reversible and $$n_{f,b}$$ are the number of executed forward and backward events of one type of reaction. The quasi-equilibrated reaction scaling factor is updated every $$N_s$$ steps:

\begin{equation} \alpha_m = N_f\dfrac{2r_S}{r_{m,f}+r_{m,b}} \end{equation}

where $$N_f$$ is a factor to separate quasi-equilibrated and non-equilibrated events, $$r_S$$ is the sum of rates of non-equilibrated and non-sufficiently-executed events, and $$r_{m,f},r_{m,b}$$ are the forward and backward rates of the reaction in question. The factor of 2 accounts for the forward and backward reaction.

For the rate based scaling $$r_S$$ is defined as:

$r_S = \sum_{n\in S_A}\sum_{m\in S_B,S_N}r_{m,n}\cdot\Delta t_n$

with the sum over the observation time period ($$n$$) which the system spend in the superbasin $$S_A$$ and all processes ($$m$$) which are either non-reversible (in $$S_N$$) or not sufficiently executed (in $$S_B$$). The rate of the forward reaction $$r_{m,f}$$ is, similarly to the back reaction, ($$r_{m,b}$$) defined as:

$r_{m,f} = \sum_{n\in S_A}r_n^f\cdot\Delta t_n$

Here $$m$$ is the index of the equilibrated event of which the factor $$\alpha_m$$ is calculated.

For the scaling based on the rate constant, not the sum over the observation period is evaluated, but the mean of the rate constant over the time period spend in $$S_A$$. Changing the variable from $$r$$ to $$k$$ one gets:

$\begin{split}k_S &= \sum_{m\in S_B, S_N} \left\langle k_m \right\rangle_{n\in S_A}\\ k_{m,f} &= \left\langle k^f \right\rangle_{n\in S_A}\end{split}$

$$k_{m,b}$$ is similar to $$k_{m,f}$$, only for the corresponding back reaction.

In practice, to accelerate the MonteCoffee simulation, one needs to specify which events are each others reverse reactions. Assume that we have two event-classes named A and B that are reverse reactions to each other, and a irreversible event called Z. To accelerate the simulation, we instantiate the NeighborKMC.user_kmc.NeighborKMC object as follows

from user_kmc import NeighborKMC
from user_system import System
from user_events import A, B, Z
events = [A, B, Z]
reverse_events = {0:1}
sim = NeighborKMC(system=system,
tend=1E1,
parameters=parameters,
events=events,
rev_events=reverse_events)


Note one must ensure that the get_rate() method of all reversible events multiplies the return with self.alpha, for example as:

class A(EventBase):
...
def get_rate(self, system, site, other_site):
R = 1000. * self.params["pA"]
return self.alpha * R  # alpha needed for temporal acceleration.


After defining the different events and reversed events the question is which different parameters to use for the time acceleration. In principle we have a four-dimensional parameter space ($$\delta,n_e,N_s \mathrm{and} N_f$$). Nevertheless by some intuition and looking what this parameters actually are and do, we can reduce the necessity of converging all of them.

First we take a look on the parameter $$\delta$$. It defines when a process has reached equilibrium. Because of statistical fluctuations it is sensible to choose a value within $$\delta \in [0.1,0.3]$$, but in principle any of these values is fine. The second easily to determine parameter is $$n_e$$, which is the number of least executions of either the forward or the reverse reaction for an event to be registered as quasi-equilibrated. From a conceptional point of view it is reasonable to choose $$50 < n_e < 500$$. With 50 necessary executions before an event is accounted as quasi equilibrated one ensures the kinetic consistency but if $$n_e$$ is too large the simulation time is unnecessarily prolonged.

The other two parameters $$N_s$$ and $$N_f$$ can not be so easily guessed and thus it is recommended to converge them separately. Here $$N_s$$ is the number of steps after which we check if quasi-equilibrated events are sufficiently executed and if that is the case, update the time acceleration parameter self.alpha. In principle, one doesn’t expect too many changes, depending on $$N_s$$, thus the most important parameter to converge is $$N_f$$.

In the following we show the convergence of the time acceleration parameters for the three schemes and their efficiency compared to a normal kinetic MonteCarlo simulation at the example of CO oxidation over a Pt(111) surface. For the simulations we chose $$\delta = 0.2$$ and $$n_e = 100$$. The parameters for the CO oxidation are: $$T=800$$ K, $$p_{\mathrm{CO}}=2\cdot 10^3$$ Pa and $$p_{\mathrm{O}_2}=10^3$$ Pa. The Pt(111) surface is modeled using a $$14\times 14$$ surface, consisting of 196 surface sites and applying periodic boundary conditions. The reason to do time acceleration, is because of events on very different time scales, as can be seen in the following figure: The presented energetics in this example are based on M. Jørgensen and H. Grönbeck, ACS Catal., 7, 5054-5061 (2017) , with the diffusion energy of CO $$E_\mathrm{CO}^\mathrm{diff}=0.53$$ eV, instead of 0.08 eV to achieve comparability between the kinetic MonteCarlo and time accelerated kinetic MonteCarlo. For completeness we include in Fig. 4 also the results of the very low activation energy for diffusion of CO. In principle, the overall TOF is mainly determined by the slow events. Thus as long as CO diffusion is fast in comparison with with other events, the result is unaffected. Overall the rates for CO oxidation on Pt(111) are not so dissimilar. Thus the achieved efficiency by using time acceleration is not particularly high.

In the following, we show the convergence of the turn over frequency (TOF) per site per second with respect to $$N_f$$ and the corresponding efficiency which is defined as ratio of $$N^{80}_\mathrm{kMC}$$ to $$N^{80}_\mathrm{accel.}$$ for the three different time acceleration schemes (scale_constant, scale_rate_constant and scale_rate). The efficiency is defined so, that if $$N^{80}_\mathrm{kMC} > N^{80}_\mathrm{accel.}$$ then the accelerated simulation has a higher efficiency than the standard kinetic Monte Carlo and the efficiency is larger than 1.

Fig. 3 A shows the results obtained using the time acceleration option: scale_constant, which is the simplest scheme of the three implemented. It can be seen that up to $$N_f=100$$ the TOF obtained in the blue curve agrees well with the reference value. Only if the scaling is too harsh, the simulation becomes diffusion limited. For demonstration, the code was modified in a way to allow for overall scaling of reaction constants after their first equilibration (no rescaling). It can be seen that this method is extremely sensitive to the chosen $$N_f$$-value. Therefore, that possibility is not available generally. Fig. 3 B, gives the efficiency for the simulations resulting in the correct TOF. It can be seen that with increasing value of $$N_f$$, the number of kMC steps to form 80 CO2 molecules is drastically reduced. A) Convergence of the TOF with respect to $$N_f$$ for constant scaling with rescaling (blue) or without (orange). The black solid line gives the reference obtained from a normal kMC simulation and the dotted lines the corresponding error range. B) Speed-up of the simulation using constant scaling with rescaling compared after the execution of 80 CO2 formation reactions.

Fig. 4 A shows the convergence of the TOF with respect to $$N_f$$ for the scaling based on the rate constants. It can be seen that only for $$N_f \geq 50$$, the TOF is converged to the reference value. It should be noted, that the mean value, doesn’t hit exactly the black line, but the error of the kMC run is beside the long run time still quite large. Thus being close to the actual value is acceptable. In Fig. 4 B, the speed-up of the simulation compared to the normal kMC is shown. Clearly the effect of the used scaling is relatively small compared to Fig. 3 B. Thus for the here presented relatively simple CO2 formation, the scaling using a constant value is the most effective. Anyway that may not be the case for a complex reaction network including very different reactions. A) Convergence of the TOF with respect to $$N_f$$ for scaling based on the rate constants. The black solid line gives the reference obtained from a normal kMC simulation and the dotted lines the corresponding error range. B) Speed-up of the simulation using scaling based on the rate constants compared to the normal kMC simulation after the execution of 80 CO2 formation reactions.

The convergence with respect to the TOF using the scaling based on the rate is presented in Fig. 5 for A,B: $$N_f$$ and C,D: $$N_s$$. The TOF converges only for high $$N_f \geq 10^4$$. Nevertheless a speed-up of the simulation is still observed. Investigating for this example also the effect of $$N_s$$ on the TOF, it can bee seen that the overall TOF for $$N_s \geq 100$$ is not affected by any of the chosen values. A) Convergence of the TOF with respect to $$N_f$$ for scaling based on the reaction rate. The black solid line gives the reference obtained from a normal kMC simulation and the dotted lines the corresponding error range. B) Speed-up of the simulation using scaling based on the reaction rate compared to the normal kMC simulation after the execution of 80 CO2 formation reactions. C) Convergence of the TOF with respect to $$N_s$$ for scaling based on the reaction rate. The black solid line gives the reference obtained from a normal kMC simulation and the dotted lines the corresponding error range. D) Speed-up of the simulation using scaling based on the reaction rate compared to the normal kMC simulation for different $$N_s$$, after the execution of 80 CO2 formation reactions.

As here presented, the various acceleration schemes are very different in respect of their parameters, but having the same general effect: speeding up the simulation compared to a standard kMC simulation. Which acceleration scheme to use depends solely on the system in hand and personal preferences.

In the following we list some additional tips how to handle fast events:

Tip 1 to slow down identical reactions, say CO adsorption, on different types of sites separately, simply define two event-classes, for example COAdsCorner and COAdsEdge.

Tip 2 diffusion events are often fast. They are, in principle, their own reverse and can be added as

>>> reverse_events = {3:3}


Tip 3 although the acceleration scheme is implemented in MonteCoffee, it may be beneficial to add a constant offset to the diffusion barriers to slow them down further. This should, however, be done carefully.

Adsorbate-adsorbate interactions are implemented by checking site-occupations dynamically in the get_rate() methods of the events. As an example, here we define a method get_repulsion() in user_energy.py to return the repulsion between species 1 and 2, represented as a list of lists:

def get_repulsion(cov_self, cov_NN, stype):

stype_factor = 0.5 if stype in [0, 1] else 1.0
repulsion = 0.
ECOCO = 0.19  #  How CO affects CO (eV).
EOO = 0.32  # How O affects O (eV) - double since it is called from get barrier of O2.

ECOO = 0.3  # How CO affects O (eV).
EOCO = 0.3  # How O affects CO (eV).

HInttwo = [[0., 0., 0.], [0., ECOCO, EOCO],
[0., ECOO, EOO]]  # Two body interaction Hamiltonian 3x3 because 0 = empty.

for j in cov_NN:  # For each covered Neighbor, give a repulsion:
repulsion += HInttwo[cov_self][j]

repulsion *= stype_factor  # Half the repulsion if edge/corner.

return repulsion


Here, if stype is 0 or 1 (corner or edge), the repulsions are halved in strength. The repulsions can be implemented in the get_rate() methods of the events as

class COOxEvent(EventBase):
.
.
.
def get_rate(self, system, site, other_site):
stype = system.sites[site].stype
stype_other = system.sites[other_site].stype

# Find the repulsion
# from the site-occupations:
Ncovs = [system.sites[n].covered for n in
system.neighbors[site]]
Nothercovs = [system.sites[n].covered for n
in system.neighbors[other_site]]

# Subtract repulsion from binding energies:
ECO -= get_repulsion(1, Ncovs, stype)
EO -= get_repulsion(2, Nothercovs, stype_other)

# Get activation energy from BEP relation:
Ea = max(0., get_Ea(ECO, EO)) # not < 0.

return self.alpha * self.Zratio * np.exp(-Ea /
(kB * self.params['T'])) * kB * self.params['T'] / h


If next nearest neighbor interactions are to be implemented, this example should be extended to access the neighbors of the neighbors. If this is done, remember to change nninteractions in kMC_options.cfg to update newly enabled events properly after each event execution.

### Calculating a turnover frequency¶

The turnover frequency (TOF) can be calculated from the number of times a product is formed per site and second. The same procedure can be followed for the individual elementary step rates. The script analyze_tof.py provides a complete example of how the TOF can be calculated.

Assume that the event list is instantiated as

from user_events import A, B, X, Z
events = [A, B, X, Z]


Where X is the forward reaction of a step that generates the product molecule, and Z is the reverse reaction that consumes one product molecule. To calculate the system’s overall TOF, we load in the time, and events that were executed (see output)

import numpy as np
stypes = np.loadtxt("stypes.txt") # for number of sites here
Nsites = len(stypes)

TOF_global = (evs_exec-evs_exec)/(time[-1]-time)/float(Nsites)


We may want a TOF for each type of site and to discard the first half of the simulation, which may be out of steady-state:

Nevents = 4
sid_ev = np.loadtxt("sid_ev.txt").reshape(-1, Nsites, Nevents)

Nhalf = int(np.round(len(sid_ev)/2.,0))
sid_ev = sid_ev[Nhalf:]

dt =  time[-1]-time[int(np.round(len(time)/2.,0))]
tofs_st = np.zeros(len(list(set(stypes))))

for time_chunk in sid_ev:
for n, st in enumerate(list(set(stypes))):
sids_st = [i for i, s in enumerate(stypes) if s == st] #  site-indices of stype==st.
Nst = float(len(sids_st)) #  number of sites with the current stype == st.
tofs_st[n] += sum([time_chunk[i] - time_chunk[i]
for i in sids_st]) / (dt*Nst) # TOF of the stype


Here, the last half of the simulation is used, and for all unique types of sites the indices are noted. Then the number of sites with the current stype is noted, and finally the TOF is calculated.

Statistical averaging should be done to address the convergence of the TOF. This can be done by running multiple identical simulations in parallel.

### Restarting your simulation¶

To collect sufficient statistic for the kinetic MonteCarlo simulation it may be important to restart MonteCoffee after a certain time has elapsed, especially on computer clusters with short queue times.

Thus, the user part can be adjusted in the following way, first in the main-run file:

real_t_end = 10 #Real end time of simulation to restart in s
# Instantiate simulator object, now including the simulation end time.
sim = NeighborKMC(system=p, tend=tend,
real_t_end = real_t_end,
parameters=parameters,
events=events,
rev_events=reverse_events)


In the ‘user_kmc.py’ two new functions need to be defined, serialize and deserialize and the package pickle imported:

import pickle, os, time

def serialize(self,filename):
"""Ads the possibility to dump self object"""
with open(filename, 'wb') as f:
pickle.dump(self.__dict__,f)

def deserialize(self,filename):
"""Reads the self object from the file"""
with open(filename, 'rb') as f:


Additionally, the variable real_t_end has to be added to the __init__ of the simulation:

def __init__(self, system, tend, real_t_end = (96*60*60), parameters={}, events=[], rev_events={}):
self.events = [ev(parameters) for ev in events]
self.reverses = None # Set later
self.real_t_end = real_t_end


The time module is used to follow the real time of the simulation. To use the real time as second break condition of the simulation, it is included in the while-clause. At the end of the while-clause the self-object with the system state is dumped as pickle-file.

start_time = time.time()  # save start time of simulation
if os.path.exists('data.pck'):   # if restart file exists, load self-object
self.deserialize('data.pck')

log.dump_point(self.stepNMC, self.t, self.evs_exec)
while  time.time() < start_time + self.real_t_end and self.t < self.tend:
self.frm_step()

self.serialize('data.pck') # dump self-object


Please notice: The time used here is the bare simulation time. Thus it must be reduced by any pre-process time to initialize the system.

## Options¶

The configurations and options in MonteCoffee are specified in the file kMC_options.cfg, and are loaded by load_options(). The options contained in kMC_options.cfg are:

• nspecies (int): how many species are included in the simulation. This integer decides how many coverages are written to the log.
• nninteractions (int): describes the extent of the nearest neighbor interactions. Specifying 1 makes the local search for newly enabled events extent to second nearest neighbor.
• savesteps (int): Number of Monte Carlo steps between saving the .txt files.
• logsteps (int): Number of steps to skip when saving .txt files and printing the log. Setting it to 1000 makes the log output and code save every 1000th simulation step.
• tinfinity (float): The time to consider infinite. This time is set for the occurrence time of impossible events.
• savecovs (bool): If True, the site-occupations (coverages) are saved periodically to coverages.txt.
• verbose (bool): If True, prints verbose information.
• write_atoms (bool): If True, write out traj files of coverages (may have to be customized depending on the system).

Options solely related to the acceleration of kMC simulations (See accelerating):

• use_scaling_algorithm (str): If None or False, no time acceleration will be used. In the current version the following three schemes are available:

• scale_constant – scaling with a constant factor
• scale_rate – scaling based on the reaction rates
• scale_rate_constant – scaling based on the reaction rate constants
• delta (float): The reversibility tolerance for quasi-equilibrated events. This determines the criterion for whether reactions are quasi-equilibrated. Setting delta=0.25 implies that there must be no more than 25 percent difference in rates of a forward and backward reaction to deem it quasi-equilibrated.

• nf (int): A factor that separates fast and slow events. nf = 1000 means that quasi-equilibrated events are slowed down to maximally 1000 times the non-equilibrated events.

• ns (int): The frequency of adjusting barriers for quasi-equilibrated events.

• ne (int): The number of events executed for a reaction in forward and backward direction for the reaction to be considered equilibrated.

## User interface (API)¶

### NeighborKMC package¶

#### NeighborKMC.base package¶

The listed modules are the main part of the NeighborKMC package and should be only touched if one knows what one is doing. All adaptations to specific systems are done in the user files and the main controll file.

##### NeighborKMC.base.basin module¶

Contains methods for performing temporal acceleration of kMC simulations.

The methods here aid in performing the acceleration of the kMC simulations in MonteCoffee. This is mainly based on the work of Dybeck et al. (https://doi.org/10.1021/acs.jctc.6b00859)

NeighborKMC.base.basin.leave_superbasin(sim)[source]

Leaves the superbasin.

Resets all rate-scalings and statistics connected to the superbasin. The sufficiently executed event list is reset with rescaling.

Parameters: sim (NeighborKMC) – main simulator object to perform rescaling of events for.
NeighborKMC.base.basin.rescaling(sim)[source]

Rescales the times of occurrences for events.

Rescales the times according to each quasi-equilibrated and sufficiently executed events alpha.

Parameters: sim (NeighborKMC) – main simulator object to perform rescaling of events for.
NeighborKMC.base.basin.scale_constant(sim)[source]

Rate based on a constant frequency factor for deeceleration as used for example by Hoffmann and Bligaard

NeighborKMC.base.basin.scale_rate(sim)[source]

Rate based superbasin escape time.

Calculates superbasin escape time according to non-equilibrated event rates escaping the superbasin.

c.f. the generalized temporal acceleration scheme of Dybeck et al.

Parameters: sim (NeighborKMC) – main simulator object to perform rescaling of events for. noneqevents (list(int)) – The indices of events that are not in equilibrium, according to the loading order passed to sim. The estimated superbasin escape-rate. float
NeighborKMC.base.basin.scale_rate_constant(sim)[source]

Rates based on the mean of the current observation period

NeighborKMC.base.basin.superbasin(sim, evtype, dt)[source]

Scales rates or leaves the current superbasin.

Based on the current Monte Carlo step, the method determines if the superbasin is left or the fast quasi-equilibrated events should be slowed down.

Keeps track and performs barrier adjustments, of the generalized temporal acceleration scheme of Dybeck et al.

Parameters: sim (NeighborKMC) – main simulator object to perform rescaling of events for. evtype (int) – The index of the event-type of the currently attempted Monte Carlo step. dt (float) – The time-step of the currently attempted Monte Carlo step. do_scaling (bool) – Bool if the rate constants are to be rescaled or not. Warning – If the time-step is negative.
##### NeighborKMC.base.events module¶

Defines the EventBase class.

The EventBase class is defined here as a template-class to derive events that can be stored in user_events.py.

The class defines methods that will throw an error if implemented wrongly, or are not implemented, in the derived classes.

Module
user_events
class NeighborKMC.base.events.EventBase(params={}, name='no-name')[source]

Bases: object

Template class for events.

Stores a list of parameters related to reaction events.

The class is only used as a parent, and is in this sense purely abstract.

params

Parameter dict dumped at the beginning of the log.

Type: dict
alpha

The slowing down factor that is adsjusted when the reaction is accelerated. This factor is set to 1 upon instantiation and varies periodically during simulation. (See Module: NeighborKMC.base.basin)

Type: int
diffev

Is the event a diffusion event. This can be used to make special rules for diffusion events.

Type: bool
do_event(system, site, other_site)[source]

Template method to perform the event.

Method needs to be overridden in user_events.py. Should change system site coverages by changing system.sites[i_site].covered and system.sites[other_site].covered.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list.
get_involve_other()[source]

Template method to decide if neighboring atoms are effected by event.

Method needs to be overridden in user_events.py. Should return True if event effects neighboring atoms and False if not (e.g. single atom adsorption).

get_rate(system, site, other_site)[source]

Template method to determine the rate constant.

Method needs to be overridden in user_events.py. Should return the reaction rate on site number i_site, and i_other for multi-site reactions.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. Rate constant of event – The rate-constant given the current occupation patterns around the site-pair i_site and i_other. float
possible(system, site, other_site)[source]

Template method to determine if event is possible.

Method needs to be overridden in user_events.py. Should return True if an event is possible on site number i_site and possible a neighbor site i_other, given the current site occupations.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. possible – True if event is possible on site-pair i_site and i_other. False if event is impossible on site-pair i_site and i_other. bool
##### NeighborKMC.base.kmc module¶

Defines the NeighborKMCBase class.

The methods are used to perform kMC simulations with the first reaction method.

class NeighborKMC.base.kmc.NeighborKMCBase(system, tend, parameters={})[source]

Bases: object

Main class for performing MonteCoffee simulations.

Assigns a system to the simulation, stores parameters, and reads in software configuration from the separate file kMC_options.cfg. Then it sets the time equal to zero and prepares to perform frm kinetic Monte Carlo simulations.

Attributes:

system: System
The system instance to perform the simulation on.
tend: float
Simulation end-time, given in seconds.
parameters: dict
parameters used to calculate rate-constants and to dump to log files. Example: parameters = {‘pCO’:1E2,’T’:700,’Note’:’Test simulation’}
t: float
Current simulation time in seconds.

Attributes used to keep track of events (where and when they happen)

siteslist: list(int)
The list of sites for each specific event. This list has the length: len(self.events)*len(self.sites)*len(self.sites).
other_sitelist: list(int)
The list of neighbor sites for each specific event. This list has the length: len(self.events)*len(self.sites)*len(self.sites).
lastsel: int
The int of the last selected site.
lastother: int
The int of the last selected neighbor site.
rindex: list(list(list(int)))):
The index of the specific events in lists like self.frm_times. For example to find the indices of site no i and event no j and neighbor number k to site i, call rindex[i][j][k].
possible_evs: list(int):
List of events that are possible, used for superbasin algorithms. This list has the length: len(self.events)*len(self.sites)*len(self.sites).
evs: numpy.ndarray(int):
The event numbers for each specific event. This list has the length: len(self.events)*len(self.sites)*len(self.sites).
rs: numpy.ndarray(float)
Rate constants of specific events. This list has the length: len(self.events)*len(self.sites)*len(self.sites).
wheres: list(list(int)):
List of all the positions of the event-types in the lists with length len(self.events)*len(self.sites)*len(self.sites). To find all site-indices where event i happens, call wheres[i].
involve_other: bool:
-False if the event happens only on one specific site -True if the event modifies two or more sites

Statistics counting attributes used to log and write output

SaveSteps: int
The number of Monte Carlo steps between saving the .txt files.
LogSteps: int
The number of Monte Carlo steps between logging steps.
tinfinity: float
What time to put impossible events to.
Nspecies: int
How many different types of species are in the simulation. Used to print and log.
nninter: int
How deep is the nearest-neighbor interaction (depth of effect of event on neighbor properties)
verbose: bool
If True, the code prints verbose information.
save_coverages: bool
If True, coverages are saved to coverages.txt and the site, othersite and event evolution in detail_site_event_evol.hdf5. This can result in large files.
write_atoms: bool
If True, the surface atoms are written with the step number in the filename. It has to be adjusted for adsorption species individually.
times: list(float)
List of times for each logged monte carlo steps in self.MCstep
MCstep: list(int)
List of Monte Carlo step numbers logged.
Nsites: int
The number of sites in self.system
Nstypes: int
The number of distinct site-types.
covered: list(list(int))
A list of site-occupations, of each site for each logged step. To find the site-occupation at step no self.MCstep[i] and site j, call covered[i][j].
system_evolution: list(list())
List which contains a list of the fired event with at site, other site and time
used_ijk: list(tuples(site,event,othersite))
List of tuples representing the unique neighbor-event pairs avoiding double counting.

Superbasin attributes related to temporal acceleration

equilEV: list(int)
A list of the event-indices that are quasi-equilibrated.
Suffex: list(int)
A list of the event-indices that are quasi-equilibrated and sufficiently executed.
tgen: list(float)
A list of when each specific event was generated. This list has the length: len(self.events)*len(self.sites)*len(self.sites).
us: list(float)
A list of random deviates used when each specific event was generated. This list has the length: len(self.events)*len(self.sites)*len(self.sites).
rescale: bool
Defines if the rates have to be rescaled and frm timelist updated.
r_S: numpy.ndarray(float)
The cumulative rates in the current superbasin.
dt_S: list(float)
The time-steps taken in the current superbasin.
nem: numpy.ndarray(int)
The number of times an event was performed the last self.ne steps
use_scaling_algorithm: bool
Defines if the rate constants are scaled in any way or not
delta: float
Reversibility tolerance to determine if reactions have become quasi-equilibrated.
Nf: int
The average number of steps a quasi-equilibrated event should be observed in each superbasin.
Ns: int
The frequency of barrier scaling.
ne: int
The minimum number of times to see a quasi-equilibrated event in each superbasin.
isup: int
How many steps were taken in the current superbasin.

Module NeighborKMC.base.basin
for documentation about the superbasin.

user_kmc - files in the tutorial/examples folder for aditional specifications.

frm_init()[source]

Prepare to perform FRM simulation.

Initializes empty rate and event lists to bookkeep the FRM algorithm. The initial times of occurrence for each event is also calculated and stored.

frm_step()[source]

Takes a Monte Carlo Step.

Takes a monte carlo step by performing the chronologically next possible event, which has index self.frm_arg in the list self.frm_times.

Raises: Warning: – If an impossible event is attempted. Usually due to an infinite time-step.
frm_update()[source]

Updates the FRM related lists.

Method updates the event list locally around the site where the last event happened. This is done by determining if new events have become possible as a result of performing the last event.

Events that are no longer possible because of executring the previous event are flagged as impossibe and their time is set to infinity.

load_events()[source]

Loads events (abstract method).

This method must be overridden by the child class in user_kmc.NeighborKMC.

Raises: NotImplementedError: – If called.
load_options()[source]

Loads all options set in kMC_options.cfg.

Instantiates a configuration parser, and loads in all options from kMC_options.cfg.

run_kmc()[source]

Runs the kMC simulation (abstract method)

This method must be overridden by the child class in user_kmc.NeighborKMC.

Raises: NotImplementedError: – If called.
##### NeighborKMC.base.logging module¶

Defines the Log class to log results of a kMC run.

class NeighborKMC.base.logging.Log(parameters, accelparams={'on': False})[source]

Bases: object

Handles logging of kMC simulations.

Initializes a filename based on the CPU date and time. All passed parameters will be written to the log.

Parameters: parameters (dict) – parameters to dump at the beginning of a log. For example >>> dict = {'T':300, 'pCO':1E2} 

Examples

Simply instantiate a Log and write a line as

>>> log = Log(parameters = _params)
>>> log.write_line("This is a line!")


Dump a simulation point to the log:

>>> log.dump_point(step=100, sim_time=1E-9, ev_called=[10,90,0,0,0])


Module
NeighborKMC.base.kmc
dump_point(step, sim_time, ev_called)[source]

Writes a simulation point to the log.

Method writes the Monte Carlo step number step, the time in the MC simulation sim_time, and the number of event calls ev_called.

Parameters: step (int) – The Monte Carlo step number. sim_time (float) – The simulation time in seconds. ev_called (list(int)) – The number of times each event is called during simulation. For example [N_CO_ads,N_CO_des,…].
save_atoms(filename)[source]

Writes tagged ase.Atoms to file.

Writes self.atom_cfgs to file with path filename. The variable self.atom_cfgs can be tagged with coverages or augmented with molecules near the sites to visualize the reaction trajectory. This is currently not implemented.

Parameters: filename (str) – Path to file.
save_txt(save_step=1000.0)[source]

Saves txt files containing the simulation data.

Saves the number of events executed on the different types of sites, the time vs mcstep, the site-types, and optionally the coverages if self.covered is True.

Growing lists are cleaned from memory.

write_line(string)[source]

Writes a line to the log.

Appends a string to the end of the log on its own line.

Parameters: string (str) – string to write to log.
##### NeighborKMC.base.sites module¶

Defines the SiteBase Class.

This class is used as a parent for the Site class defined in user_sites.py.

Module
user_sites
class NeighborKMC.base.sites.SiteBase(stype=0, covered=0, ind=[], lattice_pos=None)[source]

Bases: object

Class that templates site objects.

Assigns a site type (stype) to the site, the species covers the site (covered), atomic indices that constitute the site (ind), and the sites that are nearest-neighbors (neighbors).

stype

The site type. The user must decide what that implies. Example: 0 ~ (111)-facet-ontop, 1 ~ edge-ontop …

Type: int
covered

The species that covers the site. The user must decide what the integer implies. Example: 0 ~ empty-site, 1 ~ Oxygen covered, 2 ~ CO covered.

Type: int
ind

The atomic-indices c.f. an ase.Atoms object that constitute the site. This is convenient to define for later visualization purposes.

Type: list(int)
lattice_pos

The lattice position of the site. Can be used for systems that obey periodic boundary conditions, and to determine neighbor-lists.

Type: list(int)
##### NeighborKMC.base.system module¶

Defines the SystemBase class.

The module defines a class used to template the System class defined in user_system.

Module
sites
Module
user_sites
class NeighborKMC.base.system.SystemBase(sites, atoms=None)[source]

Bases: object

Defines a system class to perform kMC.

Method assigns an ASE.Atoms object (atoms) to the object and assigns a list of sites (sites).

A neighbor list (neighbors) is initialized from the sites, which is checked for inconsistencies by the method verify_nlist().

atoms

Can be passed to connect an ASE atoms object to the system.

Type: ase.atoms
sites

The list of sites that constitute the system.

Type: Site
find_nn_recurse(sim, update_sites, recursion=0)[source]

Deep search of first nearest neighbors.

Calculates the first nearest neighbors for a list of site_indices (update_sites).

For example, when passing update_sites = [0,1,2], the method returns [0,1,2,N neighbor 0 of site 0, Neighbor 1 of site 0, …, Neighbor 0 of site 1, …].

The method is calling itself recursively until the lattice is updated, c.f. the locality of nearest neighbor interactions.

Parameters: update_sites (list(int)) – The site indices to return neighborlist of. recursion (int) – The recursive level of which function was called, because the method calls itself, for example in base.kmc.frm_update(). out – An update to the list update_sites where the neighbors to update_sites are added. list(int)

kmc.NeighborKMC.frm_update()

get_coverages(N_species)[source]

Gets the site-occupations at the present moment.

Returns: cov list(list(float)) (a list of site-occupations for each species) and all sites. Thus to find the coverage of species i on site number j one calls ret[i][j].
get_ncovs(i_site)[source]

Gets the coverage on nearest neighbor sites.

Retrieves and returns the occupations of the nearest neighbor sites to the site with index i_site in self.sites.

Parameters: i_site (int) – Index of the site in self.sites. covs – List of species occupying the nearest neighbor sites. list(int)
verify_nlist()[source]

Tests the neighborlist for inconsistency.

The method checks if neighborlists are assymetric: If A is a neighbor to B, then B must also be present in the neighborlist of A.

Raises: Warning: – If the neighborlist is assymetric.
##### Module contents¶

Base classes used to template and run behind the scenes.

#### NeighborKMC.tutorials¶

Overview of the source-code employed for the shown tutorials

##### Quick start¶
###### NeighborKMC.tutorials.quick.run_quick.py¶

Script that runs a very simple NN-kMC from one file only to have a quick start.

class NeighborKMC.tutorials.quick.run_quick.Adsorption(params)[source]
do_event(system, site, other_site)[source]

Template method to perform the event.

Method needs to be overridden in user_events.py. Should change system site coverages by changing system.sites[i_site].covered and system.sites[other_site].covered.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list.
get_involve_other()[source]

Template method to decide if neighboring atoms are effected by event.

Method needs to be overridden in user_events.py. Should return True if event effects neighboring atoms and False if not (e.g. single atom adsorption).

get_rate(system, site, other_site)[source]

Template method to determine the rate constant.

Method needs to be overridden in user_events.py. Should return the reaction rate on site number i_site, and i_other for multi-site reactions.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. Rate constant of event – The rate-constant given the current occupation patterns around the site-pair i_site and i_other. float
possible(system, site, other_site)[source]

Template method to determine if event is possible.

Method needs to be overridden in user_events.py. Should return True if an event is possible on site number i_site and possible a neighbor site i_other, given the current site occupations.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. possible – True if event is possible on site-pair i_site and i_other. False if event is impossible on site-pair i_site and i_other. bool
class NeighborKMC.tutorials.quick.run_quick.Desorption(params)[source]
do_event(system, site, other_site)[source]

Template method to perform the event.

Method needs to be overridden in user_events.py. Should change system site coverages by changing system.sites[i_site].covered and system.sites[other_site].covered.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list.
get_involve_other()[source]

Template method to decide if neighboring atoms are effected by event.

Method needs to be overridden in user_events.py. Should return True if event effects neighboring atoms and False if not (e.g. single atom adsorption).

get_rate(system, site, other_site)[source]

Template method to determine the rate constant.

Method needs to be overridden in user_events.py. Should return the reaction rate on site number i_site, and i_other for multi-site reactions.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. Rate constant of event – The rate-constant given the current occupation patterns around the site-pair i_site and i_other. float
possible(system, site, other_site)[source]

Template method to determine if event is possible.

Method needs to be overridden in user_events.py. Should return True if an event is possible on site number i_site and possible a neighbor site i_other, given the current site occupations.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. possible – True if event is possible on site-pair i_site and i_other. False if event is impossible on site-pair i_site and i_other. bool
class NeighborKMC.tutorials.quick.run_quick.simple_NKMC(system, tend, parameters={}, events=[])[source]
run_kmc()[source]

Runs the kMC simulation (abstract method)

This method must be overridden by the child class in user_kmc.NeighborKMC.

Raises: NotImplementedError: – If called.
##### Atomic adsorption and desorption¶

Script that runs the tutorial for A adsorption and desorption.

NeighborKMC.tutorials.A_ads.test.run_test()[source]

Runs the test of A adsorption and desorption over a surface.

First, constants are defined and old output files cleared. Next, the sites, events, system and simulation objects are loaded, and the simulation is performed.

Last, the results are read in from the generated.txt files, and plotted using matplotlib.

Defines the Site Class derived from base.site.SiteBase.

The site class is defined here as an interface to the base class in base.site.SiteBase, where the user can add custom tags. Custom tags can be used to evaluate the rate-constants in user_events.py

class NeighborKMC.tutorials.A_ads.user_sites.Site(stype=0, covered=0, ind=[], lattice_pos=None)[source]

A site object.

Method calls the base class constructor first. Then the user can attach custom variables to site objects, such as coordination numbers, positions, etc.

stype

The site type, user must decide what that implies. Here: 0 refers to on top adsorption

Type: int
covered

The species that covers the site, user must decide what the integer implies. Here: 0 = empty-site and 1 = A covered

Type: int
ind

The atomic-indices c.f. an ASE.Atoms object that constitute the site. This is can be used later for visualization.

Type: list(int)

Contains all user-defined event types. All user-defined events are defined here, which must be derived from the parent class EventBase.

Module NeighborKMC.base.events
for documentation about the methods possible(), get_rate(), do_event() and get_invovle_other().
class NeighborKMC.tutorials.A_ads.user_events.AAdsEvent(params)[source]

A adsorption event class.

The event is A(g) + * -> A*.

The event is possible if the site is empty.

Performing the event adds an A to the site.

class NeighborKMC.tutorials.A_ads.user_events.ADesEvent(params)[source]

A desorption event class.

The event is A* -> A(g) + * .

The event is possible if the site is A-covered.

Performing the event removes a A from the site.

class NeighborKMC.tutorials.A_ads.user_events.ADiffEvent(params)[source]

A diffusion event class.

The event is A* + * -> * + A* .

The event is possible if the site is occupied with A and the neighbouring site empty.

Performing the event changes the occupation between the site and the neighbouring site. This eventclass is not used to obtain the comparisson to the mean-field results.

Defines the NeighborKMC class used to run a MonteCoffee simulation.

The module defines the main simulation class (NeighborKMC), which is used to run the simulation. The main engine is found in base.kmc.

Module
base.kmc
class NeighborKMC.tutorials.A_ads.user_kmc.NeighborKMC(system, tend, parameters={}, events=[], rev_events={})[source]

Controls the kMC simulation.

Calls constructor of NeighborKMCBase objects, and loads in the user-specified event-list in load_events().

The variable self.evs_exec is initialized as a list to count the number of times each event-type is executed.

Parameters: system (System) – The system instance to perform the simulation on. tend (float) – Simulation end-time. parameters (dict) – Parameters used, which are dumped to the log file. Example: parameters = {‘pCO’:1E2,’T’:700,’Note’:’Test simulation’} events (list(classobj)) – A list pointing to the event classes that defines the events. The order of list is kept consistently throughout the simulation. For example, given the event classes: class AdsEvent(EventBase): def __init__(self): ... class DesEvent(EventBase): def __init__(self): ...  One should define events as a list of class names as >>> events = [AdsEvent, DesEvent]  rev_events (dict) – Specifying which events are reverse to each other, following the order self.events. This dict is used for temporal acceleration. For example, if we have an adsorption and desorption event that are each others reverse, a third non-reversible event, and a diffusion event that is its own reverse: >>> events = [AdsEvent, DesEvent, ThirdEvent, DiffusionEvent]  Then rev_events is defined as >>> rev_events = {0:1,3:3}. 

Example

Assume that we have defined a System object (system), a list of event classes (events), and the dict of reverse events (rev_events). Then a NeighborKMC object is instantiated and simulation is run as:

nkmc = NeighborKMC(system=system,
tend=1.,
parameters=params,
events=events,
rev_events=rev_events)
nkmc.run_kmc()


Module
base.kmc
Module
base.basin
load_reverses(rev_events)[source]

Prepares the reverse_event dict.

Method makes the dict self.reverses two sided, and performs a check that each event only has one reverse in the end.

Parameters: rev_events (dict) – Specifying which events are reverse to each other, as described in the constructor of NeighborKMC. Warning: – If an reversible event has more than one reverse.
run_kmc()[source]

Runs a kmc simulation.

Method starts the simulation by initializing the log, initializes lists to keep track of time and step numbers.

It saves information about the site-indices in siteids.txt, and the site-types in stypes.txt.

While the simulation runs (self.t < self.tend), Monte Carlo steps are performed by calling self.frm_step().

Every self.LogSteps, a line is added to the simulation log.

save_txt()[source]

Saves txt files containing the simulation data.

Calls the behind-the-scenes save_txt() method of the base class. The user should add any optional writes in this method, which is called every self.SaveSteps steps.

Example

Add the following line to the end of the method:

>>> np.savetxt("user_stype_ev.txt", self.stype_ev)


Defines the System Class, derived from base.system module.

The System is supposed to be a singleton that is passed to a singleton NeighborKMC object.

class NeighborKMC.tutorials.A_ads.user_system.System(atoms=None, sites=[])[source]

Class defines a collection of sites and connected atoms.

Calls the base class system.py constructor, sets the global neighborlist from the individual site’s neighborlist.

atoms

Can (optionally) be passed to connect an ASE.Atoms object to the system. This can be useful for visualization of the simulation, for example by setting the ase.Atoms tag according to coverages.

Type: ase.Atoms
sites

The sites that constitute the system.

Type: list(Site)
.. seealso::
cover_system(species, coverage)[source]

Covers the system with a certain species.

Randomly covers the system with a species, at a certain fractional coverage.

Parameters: species (int) – The species as defined by the user (e.g. empty=0,A=1). coverage (float) – The fractional coverage to load lattice with.
set_neighbors(Ncutoff, pbc=False)[source]

Sets neighborlists of self.sites by distances.

Loops through the sites and using self.atoms, the method adds neighbors to the sites that are within a neighbor-distance (Ncutoff).

Parameters: Ncutoff (float) – The cutoff distance for nearest-neighbors in angstroms pbc (bool) – If the neighborlist should be computed with periodic boundary conditions. To make a direction aperiodic, introduce a vacuum larger than Ncutoff in this direction in self.atoms. Warning: – If self.atoms is not set, because then distances cannot be used to determine neighbors.
##### Dissociative Adsorption of type B2¶

Script that runs the tutorial of B2 adsorption and desorption on a surface

NeighborKMC.tutorials.B2_ads.test.run_test()[source]

Runs the simulation of the dissociativ B2 adsorption and desorption over a surface.

First, constants are defined and old output files cleared. Next, the sites, events, system and simulation objects are loaded, and the simulation is performed.

Last, the results are read in from the generated.txt files and compared with a file containing the results from the mean-field model, and plotted using matplotlib.

Defines the Site Class derived from base.site.SiteBase.

The site class is defined here as an interface to the base class in base.site.SiteBase, where the user can add custom tags. Custom tags can be used to evaluate the rate-constants in user_events.py

class NeighborKMC.tutorials.B2_ads.user_sites.Site(stype=0, covered=0, ind=[], lattice_pos=None)[source]

A site object.

Method calls the base class constructor first. Then the user can attach custom variables to site objects, such as coordination numbers, positions, etc.

stype

The site type, user must decide what that implies. Here: 0 refers to an ontop adsorption

Type: int
covered

The species that covers the site, user must decide what the integer implies. Here: 0 = empty-site, 1 = B covered

Type: int
ind

The atomic-indices c.f. an ASE.Atoms object that constitute the site. This is can be used later for visualization.

Type: list(int)
.. seealso:: Module :py:mod:NeighborKMC.base.sites

Contains all user-defined event types.

All user-defined events are defined here, which must be derived from the parent class EventBase.

Module NeighborKMC.base.events
for documentation about the methods possible(), get_rate(), do_event() and get_invovle_other().
class NeighborKMC.tutorials.B2_ads.user_events.B2AdsEvent(params)[source]

B2 adsorption event class. The event is 2B(g) + 2* -> 2B*. The event is possible if two neighbouring sites are empty. The rate is set constant for comparison with mean-field model. Performing the event adds two B to two sites.

class NeighborKMC.tutorials.B2_ads.user_events.B2DesEvent(params)[source]

A desorption event class. The event is 2B* -> B2(g) + 2*. The event is possible if two neighbouring sites are B-covered. The rate is set constant to match the mean-field model. Performing the event removes two B from two sites.

Defines the NeighborKMC class used to run a MonteCoffee simulation.

The module defines the main simulation class (NeighborKMC), which is used to run the simulation. The main engine is found in base.kmc.

Module
base.kmc
class NeighborKMC.tutorials.B2_ads.user_kmc.NeighborKMC(system, tend, parameters={}, events=[], rev_events={})[source]

Controls the kMC simulation.

Calls constructor of NeighborKMCBase objects, and loads in the user-specified event-list in load_events().

The variable self.evs_exec is initialized as a list to count the number of times each event-type is executed.

Parameters: system (System) – The system instance to perform the simulation on. tend (float) – Simulation end-time. parameters (dict) – Parameters used, which are dumped to the log file. Example: parameters = {‘pCO’:1E2,’T’:700,’Note’:’Test simulation’} events (list(classobj)) – A list pointing to the event classes that defines the events. The order of list is kept consistently throughout the simulation. For example, given the event classes: class AdsEvent(EventBase): def __init__(self): ... class DesEvent(EventBase): def __init__(self): ...  One should define events as a list of class names as >>> events = [AdsEvent, DesEvent]  rev_events (dict) – Specifying which events are reverse to each other, following the order self.events. This dict is used for temporal acceleration. For example, if we have an adsorption and desorption event that are each others reverse, a third non-reversible event, and a diffusion event that is its own reverse: >>> events = [AdsEvent, DesEvent, ThirdEvent, DiffusionEvent]  Then rev_events is defined as >>> rev_events = {0:1,3:3}. 

Example

Assume that we have defined a System object (system), a list of event classes (events), and the dict of reverse events (rev_events). Then a NeighborKMC object is instantiated and simulation is run as:

nkmc = NeighborKMC(system=system,
tend=1.,
parameters=params,
events=events,
rev_events=rev_events)
nkmc.run_kmc()


Module
base.kmc
Module
base.basin
load_reverses(rev_events)[source]

Prepares the reverse_event dict.

Method makes the dict self.reverses two sided, and performs a check that each event only has one reverse in the end.

Parameters: rev_events (dict) – Specifying which events are reverse to each other, as described in the constructor of NeighborKMC. Warning: – If an reversible event has more than one reverse.
run_kmc()[source]

Runs a kmc simulation.

Method starts the simulation by initializing the log, initializes lists to keep track of time and step numbers.

It saves information about the site-indices in siteids.txt, and the site-types in stypes.txt.

While the simulation runs (self.t < self.tend), Monte Carlo steps are performed by calling self.frm_step().

Every self.LogSteps, a line is added to the simulation log.

save_txt()[source]

Saves txt files containing the simulation data.

Calls the behind-the-scenes save_txt() method of the base class. The user should add any optional writes in this method, which is called every self.SaveSteps steps.

Example

Add the following line to the end of the method:

>>> np.savetxt("user_stype_ev.txt", self.stype_ev)


Defines the System Class, derived from base.system module.

The System is supposed to be a singleton that is passed to a singleton NeighborKMC object.

class NeighborKMC.tutorials.B2_ads.user_system.System(atoms=None, sites=[])[source]

Class defines a collection of sites and connected atoms.

Calls the base class system.py constructor, sets the global neighborlist from the individual site’s neighborlist.

atoms

Can (optionally) be passed to connect an ASE.Atoms object to the system. This can be useful for visualization of the simulation, for example by setting the ase.Atoms tag according to coverages.

Type: ase.Atoms
sites

The sites that constitute the system.

Type: list(Site)

cover_system(species, coverage)[source]

Covers the system with a certain species.

Randomly covers the system with a species, at a certain fractional coverage.

Parameters: species (int) – The species as defined by the user (e.g. empty=0,B=1). coverage (float) – The fractional coverage to load lattice with.
set_neighbors(Ncutoff, pbc=False)[source]

Sets neighborlists of self.sites by distances.

Loops through the sites and using self.atoms, the method adds neighbors to the sites that are within a neighbor-distance (Ncutoff).

Parameters: Ncutoff (float) – The cutoff distance for nearest-neighbors in angstroms pbc (bool) – If the neighborlist should be computed with periodic boundary conditions. To make a direction aperiodic, introduce a vacuum larger than Ncutoff in this direction in self.atoms. Warning: – If self.atoms is not set, because then distances cannot be used to determine neighbors.
##### Reaction of A + B2 on a (100) surface¶
###### tutorials.A_B2_reaction.surface_100.test_surface.py¶

Script that runs a full example of CO oxidation.

NeighborKMC.tutorials.A_B2_reaction.surface_100.test_surface.run_test()[source]

Runs the test of A adsorption and desorption over a surface.

First, constants are defined and old output files cleared. Next, the sites, events, system and simulation objects are loaded, and the simulation is performed.

Last, the results are read in from the generated.txt files, and plotted using matplotlib.

###### tutorials.A_B2_reaction.surface_100.user_sites.py¶

Defines the Site Class derived from base.site.SiteBase.

The site class is defined here as an interface to the base class in base.site.SiteBase, where the user can add custom tags. Custom tags can be used to evaluate the rate-constants in user_events.py

Module
user_events
class NeighborKMC.tutorials.A_B2_reaction.surface_100.user_sites.Site(stype=0, covered=0, ind=[], lattice_pos=None)[source]

A site object.

Method calls the base class constructor first. Then the user can attach custom variables to site objects, such as coordination numbers, positions, etc.

stype

The site type, user must decide what that implies. Example: 0 ~ (111) facet ontop, 1 ~ Edge ontop …

Type: int
covered

The species that covers the site, user must decide what the integer implies. Example: 0 ~ empty-site, 1 = Oxygen covered, 2 ~ CO covered.

Type: int
ind

The atomic-indices c.f. an ASE.Atoms object that constitute the site. This is can be used later for visualization.

Type: list(int)

Module
base.sites
###### tutorials.A_B2_reaction.surface_100.user_events.py¶

Contains all user-defined event types.

All user-defined events are defined here, which must be derived from the parent class EventBase.

Module
base.events for documentation about the methods possible(), get_rate(), do_event() and get_invovle_other().
class NeighborKMC.tutorials.A_B2_reaction.surface_100.user_events.B2AdsEvent(params)[source]

B2 adsorption event class. The event is 2B(g) + 2* -> 2B*. The event is possible if two neighbouring sites are empty. The rate is set constant for comparison with mean-field model. Performing the event adds two B to two sites.

class NeighborKMC.tutorials.A_B2_reaction.surface_100.user_events.B2DesEvent(params)[source]

A desorption event class. The event is 2B* -> B2(g) + 2*. The event is possible if two neighbouring sites are B-covered. The rate comes from the forward rate and the equilibrium constant. Performing the event removes two B from two sites.

class NeighborKMC.tutorials.A_B2_reaction.surface_100.user_events.BDiffEvent(params)[source]

B diffusion event class. The event is B* + * -> * + B*.

class NeighborKMC.tutorials.A_B2_reaction.surface_100.user_events.AAdsEvent(params)[source]

A adsorption event class. The event is A(g) + * -> A*.

class NeighborKMC.tutorials.A_B2_reaction.surface_100.user_events.ADesEvent(params)[source]

A desorption event class. The event is A* -> A + *.

class NeighborKMC.tutorials.A_B2_reaction.surface_100.user_events.ADiffEvent(params)[source]

A diffusion event class. The event is A* + * -> * + A*.

class NeighborKMC.tutorials.A_B2_reaction.surface_100.user_events.ABreactEvent(params)[source]

A with B reaction event class. The event is A*+B* -> AB(g) + 2*.

###### tutorials.A_B2_reaction.surface_100.user_kmc.py¶

Defines the NeighborKMC class used to run a MonteCoffee simulation.

The module defines the main simulation class (NeighborKMC), which is used to run the simulation. The main engine is found in base.kmc.

Module
base.kmc
class NeighborKMC.tutorials.A_B2_reaction.surface_100.user_kmc.NeighborKMC(system, tend, parameters={}, events=[], rev_events={})[source]

Controls the kMC simulation.

Calls constructor of NeighborKMCBase objects, and loads in the user-specified event-list in load_events().

The variable self.evs_exec is initialized as a list to count the number of times each event-type is executed.

Parameters: system (System) – The system instance to perform the simulation on. tend (float) – Simulation end-time. parameters (dict) – Parameters used, which are dumped to the log file. Example: parameters = {‘pCO’:1E2,’T’:700,’Note’:’Test simulation’} events (list(classobj)) – A list pointing to the event classes that defines the events. The order of list is kept consistently throughout the simulation. For example, given the event classes: class AdsEvent(EventBase): def __init__(self): ... class DesEvent(EventBase): def __init__(self): ...  One should define events as a list of class names as >>> events = [AdsEvent, DesEvent]  rev_events (dict) – Specifying which events are reverse to each other, following the order self.events. This dict is used for temporal acceleration. For example, if we have an adsorption and desorption event that are each others reverse, a third non-reversible event, and a diffusion event that is its own reverse: >>> events = [AdsEvent, DesEvent, ThirdEvent, DiffusionEvent]  Then rev_events is defined as >>> rev_events = {0:1,3:3}. 

Example

Assume that we have defined a System object (system), a list of event classes (events), and the dict of reverse events (rev_events). Then a NeighborKMC object is instantiated and simulation is run as:

nkmc = NeighborKMC(system=system,
tend=1.,
parameters=params,
events=events,
rev_events=rev_events)
nkmc.run_kmc()


Module
base.kmc
Module
base.basin
load_reverses(rev_events)[source]

Prepares the reverse_event dict.

Method makes the dict self.reverses two sided, and performs a check that each event only has one reverse in the end.

Parameters: rev_events (dict) – Specifying which events are reverse to each other, as described in the constructor of NeighborKMC. Warning: – If an reversible event has more than one reverse.
run_kmc()[source]

Runs a kmc simulation.

Method starts the simulation by initializing the log, initializes lists to keep track of time and step numbers.

It saves information about the site-indices in siteids.txt, and the site-types in stypes.txt.

While the simulation runs (self.t < self.tend), Monte Carlo steps are performed by calling self.frm_step().

Every self.LogSteps, a line is added to the simulation log.

save_txt()[source]

Saves txt files containing the simulation data.

Calls the behind-the-scenes save_txt() method of the base class. The user should add any optional writes in this method, which is called every self.SaveSteps steps.

Example

Add the following line to the end of the method:

>>> np.savetxt("user_stype_ev.txt", self.stype_ev)

###### tutorials.A_B2_reaction.surface_100.user_system.py¶

Defines the System Class, derived from base.system module.

The System is supposed to be a singleton that is passed to a singleton NeighborKMC object.

Module
base.system
Module
user_sites
class NeighborKMC.tutorials.A_B2_reaction.surface_100.user_system.System(atoms=None, sites=[])[source]

Class defines a collection of sites and connected atoms.

Calls the base class system.py constructor, sets the global neighborlist from the individual site’s neighborlist.

atoms

Can (optionally) be passed to connect an ASE.Atoms object to the system. This can be useful for visualization of the simulation, for example by setting the ase.Atoms tag according to coverages.

Type: ase.Atoms
sites

The sites that constitute the system.

Type: list(Site)

Module
base.system
cover_system(species, coverage)[source]

Covers the system with a certain species.

Randomly covers the system with a species, at a certain fractional coverage.

Parameters: species (int) – The species as defined by the user (e.g. empty=0,CO=1). coverage (float) – The fractional coverage to load lattice with.
set_neighbors(Ncutoff, pbc=False)[source]

Sets neighborlists of self.sites by distances.

Loops through the sites and using self.atoms, the method adds neighbors to the sites that are within a neighbor-distance (Ncutoff).

Parameters: Ncutoff (float) – The cutoff distance for nearest-neighbors in angstroms pbc (bool) – If the neighborlist should be computed with periodic boundary conditions. To make a direction aperiodic, introduce a vacuum larger than Ncutoff in this direction in self.atoms. Warning: – If self.atoms is not set, because then distances cannot be used to determine neighbors.
##### Reaction of A + B2 on a NP¶
###### tutorials.A_B2_reaction.truncated_octahedron.test.py¶

Script that runs A+B2 reaction over a truncated Octahedron

NeighborKMC.tutorials.A_B2_reaction.truncated_octahedron.test.run_test()[source]

Runs the test of A adsorption and desorption over a surface.

First, constants are defined and old output files cleared. Next, the sites, events, system and simulation objects are loaded, and the simulation is performed.

Last, the results are read in from the generated.txt files, and plotted using matplotlib.

###### tutorials.A_B2_reaction.truncated_octahedron.user_sites.py¶

Defines the Site Class derived from base.site.SiteBase.

The site class is defined here as an interface to the base class in base.site.SiteBase, where the user can add custom tags. Custom tags can be used to evaluate the rate-constants in user_events.py

Module
user_events
class NeighborKMC.tutorials.A_B2_reaction.truncated_octahedron.user_sites.Site(stype=0, covered=0, ind=[], lattice_pos=None)[source]

A site object.

Method calls the base class constructor first. Then the user can attach custom variables to site objects, such as coordination numbers, positions, etc.

stype

The site type, user must decide what that implies. Example: 0 ~ (111) facet ontop, 1 ~ Edge ontop …

Type: int
covered

The species that covers the site, user must decide what the integer implies. Example: 0 ~ empty-site, 1 = Oxygen covered, 2 ~ CO covered.

Type: int
ind

The atomic-indices c.f. an ASE.Atoms object that constitute the site. This is can be used later for visualization.

Type: list(int)

Module
base.sites
###### tutorials.A_B2_reaction.truncated_octahedron.user_events.py¶

Contains all user-defined event types.

All user-defined events are defined here, which must be derived from the parent class EventBase.

Module
base.events for documentation about the methods possible(), get_rate(), do_event() and get_invovle_other().
class NeighborKMC.tutorials.A_B2_reaction.truncated_octahedron.user_events.B2AdsEvent(params)[source]

B2 adsorption event class. The event is 2B(g) + 2* -> 2B*. The event is possible if two neighbouring sites are empty. The rate is set constant for comparison with mean-field model. Performing the event adds two B to two sites.

class NeighborKMC.tutorials.A_B2_reaction.truncated_octahedron.user_events.B2DesEvent(params)[source]

A desorption event class. The event is 2B* -> B2(g) + 2*. The event is possible if two neighbouring sites are B-covered. The rate comes from the forward rate and the equilibrium constant. Performing the event removes two B from two sites.

class NeighborKMC.tutorials.A_B2_reaction.truncated_octahedron.user_events.BDiffEvent(params)[source]

B diffusion event class. The event is B* + * -> * + B*.

class NeighborKMC.tutorials.A_B2_reaction.truncated_octahedron.user_events.AAdsEvent(params)[source]

A adsorption event class. The event is A(g) + * -> A*.

class NeighborKMC.tutorials.A_B2_reaction.truncated_octahedron.user_events.ADesEvent(params)[source]

A desorption event class. The event is A* -> A + *.

class NeighborKMC.tutorials.A_B2_reaction.truncated_octahedron.user_events.ADiffEvent(params)[source]

A diffusion event class. The event is A* + * -> * + A*.

class NeighborKMC.tutorials.A_B2_reaction.truncated_octahedron.user_events.ABreactEvent(params)[source]

A with B reaction event class. The event is A*+B* -> AB(g) + 2*.

###### tutorials.A_B2_reaction.truncated_octahedron.user_kmc.py¶

Defines the NeighborKMC class used to run a MonteCoffee simulation.

The module defines the main simulation class (NeighborKMC), which is used to run the simulation. The main engine is found in base.kmc.

Module
base.kmc
class NeighborKMC.tutorials.A_B2_reaction.truncated_octahedron.user_kmc.NeighborKMC(system, tend, parameters={}, events=[], rev_events={})[source]

Controls the kMC simulation.

Calls constructor of NeighborKMCBase objects, and loads in the user-specified event-list in load_events().

The variable self.evs_exec is initialized as a list to count the number of times each event-type is executed.

Parameters: system (System) – The system instance to perform the simulation on. tend (float) – Simulation end-time. parameters (dict) – Parameters used, which are dumped to the log file. Example: parameters = {‘pCO’:1E2,’T’:700,’Note’:’Test simulation’} events (list(classobj)) – A list pointing to the event classes that defines the events. The order of list is kept consistently throughout the simulation. For example, given the event classes: class AdsEvent(EventBase): def __init__(self): ... class DesEvent(EventBase): def __init__(self): ...  One should define events as a list of class names as >>> events = [AdsEvent, DesEvent]  rev_events (dict) – Specifying which events are reverse to each other, following the order self.events. This dict is used for temporal acceleration. For example, if we have an adsorption and desorption event that are each others reverse, a third non-reversible event, and a diffusion event that is its own reverse: >>> events = [AdsEvent, DesEvent, ThirdEvent, DiffusionEvent]  Then rev_events is defined as >>> rev_events = {0:1,3:3}. 

Example

Assume that we have defined a System object (system), a list of event classes (events), and the dict of reverse events (rev_events). Then a NeighborKMC object is instantiated and simulation is run as:

nkmc = NeighborKMC(system=system,
tend=1.,
parameters=params,
events=events,
rev_events=rev_events)
nkmc.run_kmc()


Module
base.kmc
Module
base.basin
load_reverses(rev_events)[source]

Prepares the reverse_event dict.

Method makes the dict self.reverses two sided, and performs a check that each event only has one reverse in the end.

Parameters: rev_events (dict) – Specifying which events are reverse to each other, as described in the constructor of NeighborKMC. Warning: – If an reversible event has more than one reverse.
run_kmc()[source]

Runs a kmc simulation.

Method starts the simulation by initializing the log, initializes lists to keep track of time and step numbers.

It saves information about the site-indices in siteids.txt, and the site-types in stypes.txt.

While the simulation runs (self.t < self.tend), Monte Carlo steps are performed by calling self.frm_step().

Every self.LogSteps, a line is added to the simulation log.

save_txt()[source]

Saves txt files containing the simulation data.

Calls the behind-the-scenes save_txt() method of the base class. The user should add any optional writes in this method, which is called every self.SaveSteps steps.

Example

Add the following line to the end of the method:

>>> np.savetxt("user_stype_ev.txt", self.stype_ev)

###### tutorials.A_B2_reaction.truncated_octahedron.user_system.py¶

Defines the System Class, derived from base.system module.

The System is supposed to be a singleton that is passed to a singleton NeighborKMC object.

Module
base.system
class NeighborKMC.tutorials.A_B2_reaction.truncated_octahedron.user_system.System(atoms=None, sites=[])[source]

Class defines a collection of sites and connected atoms.

Calls the base class system.py constructor, sets the global neighborlist from the individual site’s neighborlist.

atoms

Can (optionally) be passed to connect an ASE.Atoms object to the system. This can be useful for visualization of the simulation, for example by setting the ase.Atoms tag according to coverages.

Type: ase.Atoms
sites

The sites that constitute the system.

Type: list(Site)

Module
base.system
cover_system(species, coverage)[source]

Covers the system with a certain species.

Randomly covers the system with a species, at a certain fractional coverage.

Parameters: species (int) – The species as defined by the user (e.g. empty=0,CO=1). coverage (float) – The fractional coverage to load lattice with.
set_neighbors(Ncutoff, pbc=False)[source]

Sets neighborlists of self.sites by distances.

Loops through the sites and using self.atoms, the method adds neighbors to the sites that are within a neighbor-distance (Ncutoff).

Parameters: Ncutoff (float) – The cutoff distance for nearest-neighbors in angstroms pbc (bool) – If the neighborlist should be computed with periodic boundary conditions. To make a direction aperiodic, introduce a vacuum larger than Ncutoff in this direction in self.atoms. Warning: – If self.atoms is not set, because then distances cannot be used to determine neighbors.

#### NeighborKMC.examples.surface package¶

##### examples.surface.test module¶

Script that runs a full example of CO oxidation.

NeighborKMC.examples.Pt_111_COOx.test.run_test()[source]

Runs the test of A adsorption and desorption over a surface.

First, constants are defined and old output files cleared. Next, the sites, events, system and simulation objects are loaded, and the simulation is performed.

Last, the results are read in from the generated.txt files, and plotted using matplotlib.

##### examples.surface.user_constants module¶

Contains physical and user-defined global constants.

These constants should be imported into other modules, e.g. user_energy.py and user_events.py to calculate rate-constant, energy and entropy.

Note

This module is optional to use and can be customized to contain all physical and global constants. Constants module has a main() method for building documentation.

NeighborKMC.examples.Pt_111_COOx.user_constants.main()[source]

Main method of user_constants.py.

Main is defined to make API documentation build.

##### examples.surface.user_energy module¶

Defines constants and methods related to the reaction energy landscapes.

This module is not mandatory to use, but is suitable to used to keep track of reaction energies that are used in user_events.py

Module
user_events
Module
user_entropy

Note

This module is optional to use and can be customized to contain all reaction energies.

NeighborKMC.examples.Pt_111_COOx.user_energy.get_Ea(ECO, EO)[source]

Computes the reaction energy barrier for CO oxidation.

Uses a BEP relation for CO oxidation, relative to Pt(111), to calculate the reaction energy barrier.

Parameters: ECO (float) – Adsorption energy of CO in eV. EO (float) – Adsorption energy of O in eV. Ea – Reaction energy barrier of CO*+O*->CO2(g) in eV. float
NeighborKMC.examples.Pt_111_COOx.user_energy.get_repulsion(cov_self, cov_NN, stype)[source]

Calculates the energy perturbation to the adsorption energy of CO or O based on: the coverage of the site (cov_self) that determines if it is calculated for CO or O, the nearest neighbor coverages (cov_NN), and site-type (stype).

Parameters: cov_self (int) – Species on site to calculate the repulsion for. cov_NN (list(int):) – Nearest neighbor occupations. stype (int) – The site-type index. repulsion – Adsorbate-adsorbate repulsion in eV. float
##### examples.surface.user_entropy module¶

Defines methods to calculate entropy of reaction rate constants.

The methods and constants defined herein are meant to be imported into events, which may be defined in user_events.py

Module
user_events

Note

This module is optional to use and can be customized to contain all reaction entropies.

NeighborKMC.examples.Pt_111_COOx.user_entropy.get_Z_CO(T, pCO)[source]

Calculates the vibration of one ideal gas CO molecule

NeighborKMC.examples.Pt_111_COOx.user_entropy.get_Z_O2(T, pO2)[source]

Calculates the vibration of ideal gas O2 moldecule.

NeighborKMC.examples.Pt_111_COOx.user_entropy.get_Ztrans(T, V, m)[source]

Calculates ideal gas translational partition function.

Method calculates the partition function for 3 free translational degrees of freedom. The input (V) is the mean-free volume, which determines the gas-pressure.

Parameters: T (float) – The temperature in K. V (float) – The ideal-gas (V = mkBT/p) volume in m^3. m (float) – Mass of the molecule. Ztrans – The free 3D-translational partition function. float
NeighborKMC.examples.Pt_111_COOx.user_entropy.get_Zvib(T, modes)[source]

Calculates quantum harmonic vibrational partition function.

Method calculates the partition function for a quantum harmonic oscillator with vibration energies modes, as per the harmonic approximation.

Parameters: T (float) – The temperature in K. modes (list(float)) – Vibrational modes in eV. Zvib – The vibrational, quantum mechanical, partition function. float
##### examples.surface.user_events module¶

Contains all user-defined event types.

All user-defined events are defined here, which must be derived from the parent class EventBase.

Module
base.events for documentation about the methods possible(), get_rate(), and do_event().
class NeighborKMC.examples.Pt_111_COOx.user_events.COAdsEvent(params)[source]

Bases: NeighborKMC.examples.Pt_111_COOx.base.events.EventBase

CO adsorption event class. The event is CO(g) + * -> CO*. The event is possible if the site is empty. The rate comes from collision theory. Performing the event adds a CO to the site.

do_event(system, site, other_site)[source]

Template method to perform the event.

Method needs to be overridden in user_events.py. Should change system site coverages by changing system.sites[i_site].covered and system.sites[other_site].covered.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list.
get_involve_other()[source]

Template method to decide if neighboring atoms are effected by event.

Method needs to be overridden in user_events.py. Should return True if event effects neighboring atoms and False if not (e.g. single atom adsorption).

get_rate(system, site, other_site)[source]

Template method to determine the rate constant.

Method needs to be overridden in user_events.py. Should return the reaction rate on site number i_site, and i_other for multi-site reactions.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. Rate constant of event – The rate-constant given the current occupation patterns around the site-pair i_site and i_other. float
possible(system, site, other_site)[source]

Template method to determine if event is possible.

Method needs to be overridden in user_events.py. Should return True if an event is possible on site number i_site and possible a neighbor site i_other, given the current site occupations.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. possible – True if event is possible on site-pair i_site and i_other. False if event is impossible on site-pair i_site and i_other. bool
class NeighborKMC.examples.Pt_111_COOx.user_events.CODesEvent(params)[source]

Bases: NeighborKMC.examples.Pt_111_COOx.base.events.EventBase

CO desorption event class. The event is CO* -> CO(g) + *. The event is possible if the site is CO-covered. The rate comes from the forward rate and the equilibrium constant. Performing the event removes a CO from the site.

do_event(system, site, other_site)[source]

Template method to perform the event.

Method needs to be overridden in user_events.py. Should change system site coverages by changing system.sites[i_site].covered and system.sites[other_site].covered.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list.
get_involve_other()[source]

Template method to decide if neighboring atoms are effected by event.

Method needs to be overridden in user_events.py. Should return True if event effects neighboring atoms and False if not (e.g. single atom adsorption).

get_rate(system, site, other_site)[source]

Template method to determine the rate constant.

Method needs to be overridden in user_events.py. Should return the reaction rate on site number i_site, and i_other for multi-site reactions.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. Rate constant of event – The rate-constant given the current occupation patterns around the site-pair i_site and i_other. float
possible(system, site, other_site)[source]

Template method to determine if event is possible.

Method needs to be overridden in user_events.py. Should return True if an event is possible on site number i_site and possible a neighbor site i_other, given the current site occupations.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. possible – True if event is possible on site-pair i_site and i_other. False if event is impossible on site-pair i_site and i_other. bool
class NeighborKMC.examples.Pt_111_COOx.user_events.CODiffEvent(params)[source]

Bases: NeighborKMC.examples.Pt_111_COOx.base.events.EventBase

CO diffusion event class. The event is CO* + * -> * + CO*. The event is possible if the site is CO-covered, and the neighbor site is empty. The rate comes from transition state theory. Performing the event removes a CO from the site, and adds it to the neighbor site.

do_event(system, site, other_site)[source]

Template method to perform the event.

Method needs to be overridden in user_events.py. Should change system site coverages by changing system.sites[i_site].covered and system.sites[other_site].covered.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list.
get_involve_other()[source]

Template method to decide if neighboring atoms are effected by event.

Method needs to be overridden in user_events.py. Should return True if event effects neighboring atoms and False if not (e.g. single atom adsorption).

get_rate(system, site, other_site)[source]

Template method to determine the rate constant.

Method needs to be overridden in user_events.py. Should return the reaction rate on site number i_site, and i_other for multi-site reactions.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. Rate constant of event – The rate-constant given the current occupation patterns around the site-pair i_site and i_other. float
possible(system, site, other_site)[source]

Template method to determine if event is possible.

Method needs to be overridden in user_events.py. Should return True if an event is possible on site number i_site and possible a neighbor site i_other, given the current site occupations.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. possible – True if event is possible on site-pair i_site and i_other. False if event is impossible on site-pair i_site and i_other. bool
class NeighborKMC.examples.Pt_111_COOx.user_events.COOxEvent(params)[source]

Bases: NeighborKMC.examples.Pt_111_COOx.base.events.EventBase

CO oxidation event class. The event is CO* + O* -> CO2(g)+2*. The event is possible if the site is CO-covered and the neighbor is O-covered. The rate comes from transition state theory. Performing the event removes a CO+O from the sites.

do_event(system, site, other_site)[source]

Template method to perform the event.

Method needs to be overridden in user_events.py. Should change system site coverages by changing system.sites[i_site].covered and system.sites[other_site].covered.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list.
get_involve_other()[source]

Template method to decide if neighboring atoms are effected by event.

Method needs to be overridden in user_events.py. Should return True if event effects neighboring atoms and False if not (e.g. single atom adsorption).

get_rate(system, site, other_site)[source]

Template method to determine the rate constant.

Method needs to be overridden in user_events.py. Should return the reaction rate on site number i_site, and i_other for multi-site reactions.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. Rate constant of event – The rate-constant given the current occupation patterns around the site-pair i_site and i_other. float
possible(system, site, other_site)[source]

Template method to determine if event is possible.

Method needs to be overridden in user_events.py. Should return True if an event is possible on site number i_site and possible a neighbor site i_other, given the current site occupations.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. possible – True if event is possible on site-pair i_site and i_other. False if event is impossible on site-pair i_site and i_other. bool
class NeighborKMC.examples.Pt_111_COOx.user_events.OAdsEvent(params)[source]

Bases: NeighborKMC.examples.Pt_111_COOx.base.events.EventBase

Oxygen adsorption event class. The event is O2(g) + 2* -> 2O*. The event is possible if two neighbor sites are empty. The rate comes from collision theory and time 0.5 because of two atoms produced. Performing the event adds O to the two empty neighbor sites.

do_event(system, site, other_site)[source]

Template method to perform the event.

Method needs to be overridden in user_events.py. Should change system site coverages by changing system.sites[i_site].covered and system.sites[other_site].covered.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list.
get_involve_other()[source]

Template method to decide if neighboring atoms are effected by event.

Method needs to be overridden in user_events.py. Should return True if event effects neighboring atoms and False if not (e.g. single atom adsorption).

get_rate(system, site, other_site)[source]

Template method to determine the rate constant.

Method needs to be overridden in user_events.py. Should return the reaction rate on site number i_site, and i_other for multi-site reactions.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. Rate constant of event – The rate-constant given the current occupation patterns around the site-pair i_site and i_other. float
possible(system, site, other_site)[source]

Template method to determine if event is possible.

Method needs to be overridden in user_events.py. Should return True if an event is possible on site number i_site and possible a neighbor site i_other, given the current site occupations.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. possible – True if event is possible on site-pair i_site and i_other. False if event is impossible on site-pair i_site and i_other. bool
class NeighborKMC.examples.Pt_111_COOx.user_events.ODesEvent(params)[source]

Bases: NeighborKMC.examples.Pt_111_COOx.base.events.EventBase

Oxygen adsorption event class. The event is 2O* -> O2(g) + 2*. The event is possible if two neighbor sites are covered with species 2 (O). The rate comes from the forward rate and the equilibrium constant. Performing the event empties the two sites by setting covered to 0.

do_event(system, site, other_site)[source]

Template method to perform the event.

Method needs to be overridden in user_events.py. Should change system site coverages by changing system.sites[i_site].covered and system.sites[other_site].covered.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list.
get_involve_other()[source]

Template method to decide if neighboring atoms are effected by event.

Method needs to be overridden in user_events.py. Should return True if event effects neighboring atoms and False if not (e.g. single atom adsorption).

get_rate(system, site, other_site)[source]

Template method to determine the rate constant.

Method needs to be overridden in user_events.py. Should return the reaction rate on site number i_site, and i_other for multi-site reactions.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. Rate constant of event – The rate-constant given the current occupation patterns around the site-pair i_site and i_other. float
possible(system, site, other_site)[source]

Template method to determine if event is possible.

Method needs to be overridden in user_events.py. Should return True if an event is possible on site number i_site and possible a neighbor site i_other, given the current site occupations.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. possible – True if event is possible on site-pair i_site and i_other. False if event is impossible on site-pair i_site and i_other. bool
class NeighborKMC.examples.Pt_111_COOx.user_events.ODiffEvent(params)[source]

Bases: NeighborKMC.examples.Pt_111_COOx.base.events.EventBase

O diffusion event class. The event is O* + * -> * + O*. The event is possible if the site is O-covered, and the neighbor site is empty. The rate comes from transition state theory. Performing the event removes a O from the site, and adds it to the other site.

do_event(system, site, other_site)[source]

Template method to perform the event.

Method needs to be overridden in user_events.py. Should change system site coverages by changing system.sites[i_site].covered and system.sites[other_site].covered.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list.
get_involve_other()[source]

Template method to decide if neighboring atoms are effected by event.

Method needs to be overridden in user_events.py. Should return True if event effects neighboring atoms and False if not (e.g. single atom adsorption).

get_rate(system, site, other_site)[source]

Template method to determine the rate constant.

Method needs to be overridden in user_events.py. Should return the reaction rate on site number i_site, and i_other for multi-site reactions.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. Rate constant of event – The rate-constant given the current occupation patterns around the site-pair i_site and i_other. float
possible(system, site, other_site)[source]

Template method to determine if event is possible.

Method needs to be overridden in user_events.py. Should return True if an event is possible on site number i_site and possible a neighbor site i_other, given the current site occupations.

Parameters: system (System) – The system, which the simulation is performed on. i_site (int) – Index of site in the system.sites list. i_other (int) – Index of other/neighbor site in the system.sites list. possible – True if event is possible on site-pair i_site and i_other. False if event is impossible on site-pair i_site and i_other. bool
##### examples.surface.user_kmc module¶

Defines the NeighborKMC class used to run a MonteCoffee simulation.

The module defines the main simulation class (NeighborKMC), which is used to run the simulation. The main engine is found in base.kmc.

Module
base.kmc
class NeighborKMC.examples.Pt_111_COOx.user_kmc.NeighborKMC(system, tend, parameters={}, events=[], rev_events={})[source]

Bases: NeighborKMC.examples.Pt_111_COOx.base.kmc.NeighborKMCBase

Controls the kMC simulation.

Calls constructor of NeighborKMCBase objects, and loads in the user-specified event-list in load_events().

The variable self.evs_exec is initialized as a list to count the number of times each event-type is executed.

Parameters: system (System) – The system instance to perform the simulation on. tend (float) – Simulation end-time. parameters (dict) – Parameters used, which are dumped to the log file. Example: parameters = {‘pCO’:1E2,’T’:700,’Note’:’Test simulation’} events (list(classobj)) – A list pointing to the event classes that defines the events. The order of list is kept consistently throughout the simulation. For example, given the event classes: class AdsEvent(EventBase): def __init__(self): ... class DesEvent(EventBase): def __init__(self): ...  One should define events as a list of class names as >>> events = [AdsEvent, DesEvent]  rev_events (dict) – Specifying which events are reverse to each other, following the order self.events. This dict is used for temporal acceleration. For example, if we have an adsorption and desorption event that are each others reverse, a third non-reversible event, and a diffusion event that is its own reverse: >>> events = [AdsEvent, DesEvent, ThirdEvent, DiffusionEvent]  Then rev_events is defined as >>> rev_events = {0:1,3:3}. 

Example

Assume that we have defined a System object (system), a list of event classes (events), and the dict of reverse events (rev_events). Then a NeighborKMC object is instantiated and simulation is run as:

nkmc = NeighborKMC(system=system,
tend=1.,
parameters=params,
events=events,
rev_events=rev_events)
nkmc.run_kmc()


Module
base.kmc
Module
base.basin
load_reverses(rev_events)[source]

Prepares the reverse_event dict.

Method makes the dict self.reverses two sided, and performs a check that each event only has one reverse in the end.

Parameters: rev_events (dict) – Specifying which events are reverse to each other, as described in the constructor of NeighborKMC. Warning: – If an reversible event has more than one reverse.
run_kmc()[source]

Runs a kmc simulation.

Method starts the simulation by initializing the log, initializes lists to keep track of time and step numbers.

It saves information about the site-indices in siteids.txt, and the site-types in stypes.txt.

While the simulation runs (self.t < self.tend), Monte Carlo steps are performed by calling self.frm_step().

Every self.LogSteps, a line is added to the simulation log.

save_txt()[source]

Saves txt files containing the simulation data.

Calls the behind-the-scenes save_txt() method of the base class. The user should add any optional writes in this method, which is called every self.SaveSteps steps.

Example

Add the following line to the end of the method:

>>> np.savetxt("user_stype_ev.txt", self.stype_ev)

##### examples.surface.user_sites module¶

Defines the Site Class derived from base.site.SiteBase.

The site class is defined here as an interface to the base class in base.site.SiteBase, where the user can add custom tags. Custom tags can be used to evaluate the rate-constants in user_events.py

Module
user_events
class NeighborKMC.examples.Pt_111_COOx.user_sites.Site(stype=0, covered=0, ind=[], lattice_pos=None)[source]

Bases: NeighborKMC.examples.Pt_111_COOx.base.sites.SiteBase

A site object.

Method calls the base class constructor first. Then the user can attach custom variables to site objects, such as coordination numbers, positions, etc.

stype

The site type, user must decide what that implies. Example: 0 ~ (111) facet ontop, 1 ~ Edge ontop …

Type: int
covered

The species that covers the site, user must decide what the integer implies. Example: 0 ~ empty-site, 1 = Oxygen covered, 2 ~ CO covered.

Type: int
ind

The atomic-indices c.f. an ASE.Atoms object that constitute the site. This is can be used later for visualization.

Type: list(int)

Module
base.sites
##### examples.surface.user_system module¶

Defines the System Class, derived from base.system module.

The System is supposed to be a singleton that is passed to a singleton NeighborKMC object.

Module
base.system
Module
user_sites
class NeighborKMC.examples.Pt_111_COOx.user_system.System(atoms=None, sites=[])[source]

Bases: NeighborKMC.examples.Pt_111_COOx.base.system.SystemBase

Class defines a collection of sites and connected atoms.

Calls the base class system.py constructor, sets the global neighborlist from the individual site’s neighborlist.

atoms

Can (optionally) be passed to connect an ASE.Atoms object to the system. This can be useful for visualization of the simulation, for example by setting the ase.Atoms tag according to coverages.

Type: ase.Atoms
sites

The sites that constitute the system.

Type: list(Site)

Module
base.system
cover_system(species, coverage)[source]

Covers the system with a certain species.

Randomly covers the system with a species, at a certain fractional coverage.

Parameters: species (int) – The species as defined by the user (e.g. empty=0,CO=1). coverage (float) – The fractional coverage to load lattice with.
set_neighbors(Ncutoff, pbc=False)[source]

Sets neighborlists of self.sites by distances.

Loops through the sites and using self.atoms, the method adds neighbors to the sites that are within a neighbor-distance (Ncutoff).

Parameters: Ncutoff (float) – The cutoff distance for nearest-neighbors in angstroms pbc (bool) – If the neighborlist should be computed with periodic boundary conditions. To make a direction aperiodic, introduce a vacuum larger than Ncutoff in this direction in self.atoms. Warning: – If self.atoms is not set, because then distances cannot be used to determine neighbors.

The user interface is essential to MonteCoffee as the code works as a programmable application.

## How to cite¶

MonteCoffee is developed at Chalmers University of Technology at the division of Chemical Physics at the group of Henrik Grönbeck. The first version was written by Mikkel Jørgensen during 2015-2019 as part of his PhD project. The code has since then developed by contributions from Noemi Bosio (since 2019) and Elisabeth M. Dietze (since 2020).

If you find MonteCoffee useful, the following journal article should be cited:

M. Jørgensen and H. Grönbeck, J. Chem. Phys. 149, 114101 (2018)

The documentation in its first version was written by Mikkel Jørgensen 2019 and since then updated by Elisabeth M. Dietze.