{'fraction_non_alphanumeric': 0.02071145499150783, 'fraction_numerical': 0.004930175504812229, 'mean_word_length': 4.720280663878019, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
KMCLib: A general framework for lattice kinetic Monte Carlo (KMC) simulations

KMCLib is a general framework for lattice kinetic Monte Carlo (KMC) simulations. The program can handle simulations of the diffusion and reaction of millions of particles in one, two, or three dimensions, and is designed to be easily extended and customized by the user to allow for the development of complex custom KMC models for specific systems without having to modify the core functionality of the program. Analysis modules and on-the-fly elementary step diffusion rate calculations can be implemented as plugins following a well-defined API. The plugin modules are loosely coupled to the core KMCLib program via the Python scripting language. KMCLib is written as a Python module with a backend C++ library. After initial compilation of the backend library KMCLib is used as a Python module; input to the program is given as a Python script executed using a standard Python interpreter. We give a detailed description of the features and implementation of the code and demonstrate its scaling behavior and parallel performance with a simple one-dimensional A-B-C lattice KMC model and a more complex three-dimensional lattice KMC model of oxygen-vacancy diffusion in a fluorite structured metal oxide. KMCLib can keep track of individual particle movements and includes tools for mean square displacement analysis, and is therefore particularly well suited for studying diffusion processes at surfaces and in solids.


Introduction
Atomic scale simulations of slowly evolving dynamics pose a great challenge in many areas of computational materials science and catalysis. The dynamics of processes such as chemical reactions on surfaces and slow diffusion events in solids often involve time-scales that are impossible to access with conventional molecular dynamics (MD) simulations. When atomic resolution is demanded to properly describe the investigated physical system, but the time-scale of the events of interest renders MD simulations unsuitable, some other method working efficiently with atomic resolution on longer time scales must be employed. There exists today a range of accelerated dynamics methods handling the dynamics in one way or another to allow for the rare events to occur much more frequently in the simulation [1]. Kinetic Monte Carlo (KMC) [2], [3] is one such method.
The idea behind KMC is to utilize the separation of time scales between on the one hand the fast molecular vibrations, and on the other hand the rare events of interest evolving the slow dynamics. Only the latter is explicitly simulated, while the former is treated as thermally equilibrated at the time one of the rare events takes place. This is theoretically justified when the separation of the vibrations and rare events time scales is large [4], [5]. For many systems of high scientific and technological interest this condition is known to be well satisfied, paving the ground for the popularity of different flavors of the KMC method [5], [6], [7], [8], [9], with common applications ranging from biological systems [8] and heterogeneous catalysis [10] to the electrochemistry in fuel-cells [11] and atomic transport mechanisms in alloys [12].
We will in this work focus on the lattice KMC method, a particularly efficient form of KMC for simulating ordered -as opposed to amorphous -systems, where the geometry of the simulated system can be conveniently described on a regular grid in space [5]. KMCLib implements the lattice KMC algorithm in the form of the algorithm known as the variable step size method (VSSM) or the n-fold way [2], [3], [5]. Although several other methods are possible [13] the VSSM methods is theoretically well founded [14] and particularly well suited for an efficient implementation for simulations of large systems.
Despite its maturity and the popularity of the approach there has been surprisingly few publicly available computer programs implementing the lattice KMC method [13], [15], [16], [17]. Although this situation is improving [18], [19], [20], most available codes are still particularly geared towards surface studies [13], [15], [20], or closed source [15] preventing straightforward customization. Due to this shortage of publicly available open sourced software together with the relative simplicity of the algorithm and the need for significant method customizations in many cases, is has not been uncommon that a whole new KMC code gets written for each new set of physical problems to solve. Many research groups each have their own custom code aimed at simulating a particular type of systems, and it is not uncommon that these codes never make it to the public. We hope to reverse this trend by providing a fast and efficient general framework for lattice KMC simulations that is easy to customize and to integrate with other software. Our hope is that KMCLib will be a useful platform for implementing custom KMC models, meeting the highest expectations for flexibility, robustness and efficiency.
KMCLib consists of a highly optimized C++ core program with a well-defined front-end API in Python, and is after compilation of the C++ code used as a Python module. Code extensions can be written in Python with no need to recompile the underlying C++ code. We have introduced the possibility for using custom rate calculations on-the-fly during the simulation, via a Python interface, that facilitates usage of any external code for calculating elementary process rates. This makes it possible to include e.g. long-range interactions or site and geometry dependent process-rates in the simulations, without modifying the main program. A plugin API for analysis onthe-fly during the simulation has also been specified to allow for user-defined analysis tools to be integrated with the code. The code is also equipped with an implementation of a recent mean square displacement (MDS) algorithm for non-equidistant time step KMC simulations [21] to facilitate diffusion studies with the code. This paper describes the main simulation algorithms implemented in KMCLib with and without the custom rate calculation modifications. The code is then described in more detail with focus on technical solutions. The scaling behavior of the implementation is demonstrated using a simple one-dimensional A-B-C model as well as with a more complex three-dimensional fluorite structure metal oxide vacancydiffusion model. Parallel performance is finally discussed and demonstrated. The KMC formalism itself will not be described in detail since it is well-known and has been described elsewhere in the literature (see e.g. [2], [3], [5], [6]). Figure 1 shows a schematic illustration of the version of the VSSM or n-fold way lattice KMC algorithm implemented in KMCLib.

Outline of the main KMC algorithm
1. Start with a set of particles at a set of points in space (lattice sites), and a list of all possible elementary processes e.g. diffusion or reaction events that can take the system of particles from the present state to the next state. Each elementary process p i is associated with a rate constant r i . 2. Match all processes with the geometry at each lattice site to construct a table of the availability of each process and obtain the total rate of each process R i by multiplying the process' rate r i with the number of sites where the process is available N i For the list of processes a vector P is constructed with the incremental total rates, such that for the j:th process the incremental value is 3. A process to perform is chosen at random, proportionally to the relative total rate of the process. This is done by drawing a random number ρ between 0 and the last value in P, P k and a binary search is performed over the elements in P to determine which process to take. Starting from the lower end, the last process p j for which P j < ρ is chosen. 4. A site to perform the process at is chosen by drawing from the list of sites where the process is available, with each site being equally likely. 5. The chosen process is performed at the chosen site by updating the geometry according to the specifications of the particular process. Each process holds information about the local geometry around a site before and after the process is performed. The before geometry is used for matching the local geometry of a lattice site with the process, while the after geometry is used to update the local geometry around the lattice site where the process is performed. 6. Time is propagated by drawing a random number between 0 and 1, and the time step ∆t to add to the simulation time is determined by with P k being the total available rate in the system. 7. Repeat from step 2 until the maximum allowed number of steps is reached.
The simulation is run for a fixed number of steps. Geometry and time step data can either be stored in a trajectory file for post-processing and analysis, or analysis can be performed onthe-fly during the simulation. Update: Collect data for on-the-fly analysis and trajectory storage.

Algorithm modifications due to custom rate calculations
An important feature of KMCLib is its ability to handle processes where the rates are not given at the start of the simulation, but are instead calculated and updated on-the-fly during the simulation. The two different ways to determine the process rates (at simulation start-up or during the simulation) can furthermore be combined, such that some of the processes gets their rates updated during the simulation while others do not. The way the rate update framework is implemented allows for great flexibility allowing the user to define the desired behavior via Python code. It is e.g. possible to couple the update of the process rates to global or local geometrical parameters, such that a temperature or pressure gradient can affect the local rates of diffusion processes differently at different sites, or that the process rates are being modified at sites in the vicinity of certain atom types or geometrical features that may arise during the simulation. The handling of custom rates requires some modifications to the above described algorithm: 1. At start-up, not only are the initial set of particles and their geometry and all possible processes given, but also a custom rate calculator. The custom rate calculator is a Python class that can receive information about a process and the local geometry around a lattice site and return an updated rate constant for the process at that particular site. The custom rate calculator has to follow a well defined API to assure compatibility with the program and is described in more detail in section 3.3.2 below, as well as in the program user manual. 2. All processes are matched with the geometry at each site to construct a table of the availability and total rate of each process. At each site where a process is possible the rate is also updated using the custom rate calculator. The total rate of each process R i is now obtained by summing the updated local rate for each site where the process is available. Note that due to the use of custom rates the same elementary process will now have a different rate at different sites. With N i being the available sites for process i, and r j being the locally updated rate at site j the total rate for process i is given by The vector P with the incremental total rate for all processes is constructed according to (2). A similar vector S of incremental rates is also constructed for each process with the available sites and their respective rates where r i is the rate at site i of the process. 3. A process to perform is chosen at random, proportionally to the relative total rate of the process. This is done in the same way as described previously by drawing a random number between 0 and P k and performing a binary search over the P vector. 4. A site to perform the chosen process at, is chosen by drawing from the list of sites where the chosen process is available. Since the rate now may vary over the sites at which the process is available it becomes necessary to draw the site with a similar procedure as when the process is chosen. A random number ρ is generated between 0 and the last value in S, S j , and a binary search is performed over the elements in S to determine which site to take. Starting from the lower end, the last site i for which S i < ρ is chosen.
The reminder of the algorithm (perform the process, propagate simulation time and collect data) is identical with the procedure described in section 2.1 above.

Outline of the mean square displacement algorithm
To collect mean square displacement (MSD) data on the fly during the simulation a history buffer is used where a fixed number of coordinates and corresponding simulation times are stored for each particle of interest. Squared displacements of the coordinates are calculated and stored in a histogram.
1. Before the simulation starts the history buffer is populated with the initial coordinate and time for each particle. 2. When a particle is moved its new coordinates along with the current time of the simulation is stored at the first position in the history buffer for that particle. Space is made for the new data in the buffer by moving each coordinate and time already in the buffer to the next position in memory. If the buffer is full the last value is thrown away. 3. The squared displacements between the first (latest added) coordinate in the buffer and all subsequent coordinates stored are calculated and binned in a squareddisplacements histogram. The number of values added to each bin is collected in a MSD bin count histogram. 4. After the simulation is finished the MSD data is obtained by dividing the value at each bin in the squareddisplacements histogram with the corresponding bin value in the MSD bin count histogram and the standard deviation of the final MSD data is estimated.
A more detailed account of the algorithm as well as the somewhat involved procedure for deriving the correct error estimate of the final MSD data has been described in detail elsewhere [21].

Description of the program
KMCLib is a general platform from which custom lattice KMC models can be implemented. The Python programming language is used to define a well tested and well documented user-interface that is flexible enough to easily define custom KMC models. This is combined with a well-structured performance optimized core program written in C++. The potential need for customization of the underlying C++ code has largely been eliminated through the use of a custom rate calculator interface that allows for algorithmic modifications of the rate determination without the need to re-compile the underlying C++ code. Analysis can be performed and data collected during the simulation, with custom analysis modules loosely coupled to KMCLib via an analysis API in Python.   Figure 2 shows a schematic picture of how the code is organized. The code consists of two major parts. One part is the deployment code for actually performing KMC simulations (figure 2 left), while the other part consists of code for testing and includes also the documentation of the Python user-interface (figure 2 right).

Code overview
The deployment code consists of a set of Python classes that defines the user-interface to the KMCLib module (including the custom rate calculator API and the on-the-fly analysis API), as well as a C++ backend library wrapped to Python using the SWIG framework. The python frontend is responsible for handling user-input and all other communication with the outside world, as well as controlling the flow of the main Monte Carlo loop and communicating with the C++ backend via the SWIG interface. All computationally heavy tasks, from step #2 through to step #5 in figure 1, including matching, picking and updating, are implemented in the C++ backend.
All software needs testing to ensure its correct functioning. If the tests are performed systematically and at all levels of the code the risk of having a malfunctioning program can be greatly reduced [22]. The automated KMCLib tests are split up in three parts: functionality tests, written using the Python unittest module testing the whole program at an aggregate level; Python unittests (also written using the Python unittest module) testing the Python code at a detailed level; and the C++ unittests implementing detailed tests of the C++ code. The C++ unittests are setup using the CPPUNIT framework. The unittests are intended to test each building block of the code at as low level as possible, while the functionality tests set up realistic but small simulations using the KMCLib Python interface much in the same way a user would do. The functionality tests represent good usage examples and can as such be used by new users to get their first simulations up and running.
The Python code is marked up with doc-strings using the Sphinx documentation format. A manual for the Python interface can be auto generated from the source code and the doc-strings using Sphinx and the provided makefile. The documentation included with KMCLib also describes the installation procedure and the code usage. All header files in C++ are documented using the doxygen format to aid in any deeper level customization and development work.

Input structure
Input to KMCLib is written as a Python script and executed using a standard Python interpreter. The input script should start with loading the KMCLib module using a from KMCLib import * statement. This will load all API classes that are needed to setup and run a lattice KMC calculation with KM-CLib. Input is then built up by specifying the geometry and all processes via the KMCConfiguration and KMCInteractions classes (see the reference manual for details).
On-the-fly custom rate calculations can at this stage be introduced simply by providing an instance of a user-given class inheriting from the KMCRateCalculatorPlugin interface when setting up the KMCInteractions object. If provided, the matching procedure later in the algorithm will use this calculator to determine an updated rate for each process at each site where the process is available. The KMCRateCalculatorPlugin Python class in turn inherits from the RateCalculator C++ class. This is achieved by declaring the C++ base class a "director" in SWIG. When a process is matched the matching routine in C++ calls the user-specified Python CustomRateCalculator via its inherited C++ interface. The custom rates interface thus allows the user to provide an altered rate for the process at the specific site in question. The capacity to use custom rates greatly increases the flexibility of the program and opens up the possibilities for performing KMC simulations with site specific rates or to take into account longer-range interactions when determining the rate of each elementary process. Note that the user is only required to write a Python extension to the program and there is no need to re-compile the underlying C++ code. Details of the custom rate calculator interface can be found in the reference manual.
After specifying the geometry and all interactions the input script sets up an instance of the main simulation object the KMCModel, and the model is run by calling the run method on the KMCModel object with an instance of the KMCControlParameters class as input holding all information on the number of steps in the simulation and when analysis and trajectory data should be gathered. If on-the-fly analysis is required a list of user-specified analysis classes, each inheriting from the KMCAnalysisPlugin class, is given as an argument to the run method. Details of how to set up the necessary objects with the correct input can be found in the program reference manual.

Features and implementation details 3.3.1. Match and update algorithm
Since the matching of local geometries with possible processes plays a central role in the performance of any lattice KMC implementation we will below describe our matching algorithm in some detail.
The KMCLib backend uses a type of neighbor lists that we will refer to as match-lists. They are used for representing local geometries around lattice sites in the configuration. When the underlying C++ configuration is constructed from the usergiven Python input a match-list is generated for each lattice site. Each match-list holds the lattice-site type information (used for labeling different elements or particle types in the simulation), the corresponding lattice site indices and the relative local coordinates for all neighboring lattices sites. Each lattice site matchlists is sorted according to its relative coordinates. The size of the match-lists is determined from the user-defined processes. The process with the largest included neighborhood will determine the cutoff used for setting up all lattice site match-lists.
KMCLib has a dual representation of the geometry. The configuration is represented on a lattice and there is information on each lattice-point about the lattice-point type. But apart from the lattice geometry we also have a list of atom-ID tags with their corresponding current types and Cartesian coordinates. The two representations are connected by the lattice that keeps track of which atom that sits on each lattice-point. This might seem trivial at first, but provides a powerful tool for keeping track of the individual particles on the lattice, which in turn is a prerequisite for any efficient diffusion calculation.
Each user given process holds information about the local configuration around a lattice site before and after the process is performed, as well as information about how each of the involved particles should be moved from the initial to the final state of the process. The local geometry information is stored in match-lists similar to the ones used at the lattice-sites. We will refer to these match-lists as the before and after process matchlists. The before process match-list will be used for matching against local geometries on the lattice, while the after process match-list will be used to update the lattice when the process is performed. Each process match-list is sorted in the same way as the lattice site match-lists. If the user has allowed match-list wild-cards the coordinates of the process match-list are compared to the corresponding set of local coordinates at the lattice and any missing entry in the process match-list is inserted with a wild-card identifier as the type. The movements of the individual particles are represented as a set of move-vectors, which in turn will be applied to the Cartesian coordinates of the atom-IDs involved in the process.
The matching of a process with a lattice site is performed by comparing the site match-list with the before process matchlist. Since the lists are sorted in the same way only the type information needs to be compared. If wild-cards are not allowed in the process match-lists the coordinates must also be compared. Performing a process is done by comparing the site match-list and the process after match-list. Any mismatch between the two will result in the corresponding types information being updated on the lattice. The atom-IDs are updated on the lattice sites according to the specifications given from the process move-vectors and the move-vectors are furthermore added to the corresponding atom-ID coordinates. In this way the lattice configuration and the atom-ID coordinates are guaranteed to be synchronized, with the important difference that the atom-ID coordinates are never checked for periodic boundary conditions.
Before the simulation enters the Monte Carlo loop the matching must be done between all sites and all processes. This results in an O(n L × n P × m) scaling behavior, where n L is the number of lattice sites, n P is the number of processes and m is the average number of match-list entries that is looped over for each match. After a process has been performed, however, there is no need to re-match the processes at all sites. Only sites that are in the vicinity of the site where the process was performed are affected by the move and a number of affected sites n A n L , which is constant with the size of the simulation will require a re-match. This results in an O(n A × n P × m) scaling behavior of the matching in the Monte Carlo loop, and the matching is thus constant O(1) with respect to the size of the simulation.

On-the-fly rate calculations
To use the custom rate calculation abilities with KM-CLib it is enough to define a class that inherits from CustomRateCalculator, and overload its rate method. The custom rate calculator is given to the interactions object using its setRateCalculator method before the interactions object is used to construct the KMC model.
When a rate update is required for a process at a specific site the custom rate calculator's rate method will be called with parameters specifying the global coordinates of the site and its local geometry and elements (types) before and after the process is applied. When the rate method is called with a site and a process it is guaranteed that the process is available at the particular site (i.e. that the site and the process match). The rate constant associated with the process is also given and it is then fully up to the user to specify any modifications to this rate from the given geometry and process input. Details of the custom rate interface can be found in the reference manual.
It should be noted that the flexibility of using custom rate updates inevitably comes with a performance penalty. It is therefore not generally recommended to implement models using a custom rate calculator that can easily be implemented with processes with fixed rates. As an illustrating example we can look the Ising model of [2] that is provided with the functionality tests, where it is implemented both using a custom rate calculator in Python and using processes with fixed rates. The equivalent simulation with fixed rates runs roughly four times faster.

Trajectory format
A simulation trajectory can be produced in two formats using KMCLib, depending on the type of simulation. For pure lattice KMC simulations that are only concerned with lattice sites and their types, without keeping track of the movement of specific atoms, a lattice trajectory format is used where the type at each lattice site is stored. Since the lattice positions are fixed during the simulation it is enough to print them only once to the trajectory file.
For simulations where move-vectors have been provided with the processes, so that the program can keep track of the individual atoms during the simulation, it is also possible to use a simple xyz format that prints the type and coordinate for each atom along the simulation. Note that although the lattice sites do not move during the simulation the individual atoms do move according to the provided move vectors.
The user can set the interval between steps to save to the trajectory. To avoid unnecessary file IO operations the trajectory data is stored in an internal buffer. The buffer size and time since last file IO is saved and the trajectory buffer is printed to file only when the buffer has reached its size limit or when a long enough time has passed since last print to file.

On-the-fly analysis
Instead of storing trajectory data to file it can be handy to run analysis directly while running the simulation. KMCLib provides the possibility to perform on-the-fly analysis via an analysis plugin interface. The user defines an object that inherits from the KMCAnalysisPlugin class and overloads the setup, registerStep and finalize methods for desired behavior.
Right before the Monte Carlo loop is entered each userprovided analysis object will in turn get access to the initialized system. This is done by calling each analysis object's setup method with the configuration and simulation start time as arguments, allowing for the analysis objects to set up internal memory and collect initial data. After a step in the Monte Carlo loop has been performed program flow control is returned to Python where all required analysis and trajectory handling is performed (the "collect data" step in figure 1). Analysis is performed by letting each user-provided analysis object in turn get access to all details of the present state of the simulation by calling each analysis object's registerStep method with the configuration and simulation time as arguments. It is then up to each analysis object to extract relevant information at this stage in the simulation. At the end of the simulation the finalize method is called on each analysis object to allow for post-processing of collected data before the program stops.
The mean square displacement analysis described in section 3.3.5 below is implemented using the analysis plugin interface. A detailed account of the analysis plugin interface can be found in the reference manual.

On-the-fly mean square displacements
To facilitate diffusion studies with KMCLib the algorithm described in [21] for calculating mean square displacements with correct error estimates on-the-fly during a simulation is implemented. To run a simulation with mean square displacement analysis it is required that move vectors are provided with the processes to allow the program to keep track of the individual atoms. An atom type to track must be provided by the user as well as parameters specifying the number of history steps to use and the histogram spacing. It is also important that the primitive unit cell of the configuration is given in correct units, since the transformation to Cartesian coordinates is done using this information before summing contributions for each time-step.
The mean square displacement analysis is implemented as an analysis plugin as described in section 3.3.4 above and several examples of its use are included in the functionality tests.

Parallelization strategy
The MPI parallelization scheme in KMCLib is based on the assumption that the matching of a process with a local geometry, including rate evaluation when required, is by far the most time-consuming part of each step in the Monte Carlo loop. This condition is particularly well satisfied for all cases where custom rate calculations involve anything more than a simple table look-up.
When the list of sites to match against the list of processes is generated after a process is performed, the process-site pairs to match are distributed over all ranks. The matching is performed in parallel and the matching result is communicated to all ranks. If a custom rate calculator is used the list of process-site pairs that need their rates to be updated is also split up over all ranks and calculated in parallel, and the updated rates are communicated to all ranks. This parallelization strategy is similar to the approach in aKMC [6] where the expensive saddle-point searches are performed in parallel. Contrary to parallelization schemes that rely on a spatial decomposition of the simulation box with several processes performed simultaneously on different ranks, our approach does not affect the overall serial algorithm. Running in serial or parallel is thus guaranteed to generate exactly the same simulation results.
To achieve a good work-load using this parallelization scheme it is necessary that the number of process-site pairs to parallelize over is much longer than, or exactly matches, the number of MPI ranks. As will be apparent in section 4.3 below, the parallelization will furthermore, as expected, work more efficiently for simulations with heavier custom rate calculations.

Random numbers
Good quality random numbers are absolutely crucial for results from any Monte Carlo method to be valid [23]. Over the years, the quality of available pseudo random numbers has improved significantly and today the Mersenne-Twister [24] algorithm is considered the gold standard for non-cryptographic applications. A publicly available library Mersenne-Twister implementation [25] is therefore included in KMCLib.
To assure that the included random number generator works correctly we tested it against the standard C++11 implementation provided with the gcc (Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3 compiler. Both pseudo random number generators generated identical results for several million consecutive numbers. Interchanging the included random number generator with the standard C++11 implementation in KMCLib furthermore produced identical simulation results. An internal wrapper was finally implemented around the pseudo random number generator in KMCLib to facilitate use of any other pseudo random number generator implementation, or the standard C++11 implementation, with minimal modifications to the code.

Timings and scaling
We use two different systems for demonstrating the scaling behavior of KMCLib. The first system is a simple onedimensional A-B-C model, where particles of type A, B or C sit on a one-dimensional lattice. Only two processes are used for this system, the flipping of A to B and B to A, leaving the C particles inert. This allows us to vary the total size of the system while keeping the number of active A+B sites constant, simply by varying the number of C particles. Both included processes have the same rate and the number of A particles was always the same as the number of B particles at the start of each simulation.
The second system is a more complex model of oxygenvacancy diffusion in a fluorite metal oxide structure with a large fraction of dopants on the metal sites and vacancies in the oxygen sub-lattice. The primitive cell is based on the cubic oxygen sub-lattice. There are two basis points in the cell, (0,0,0) and (0.5,0.5,0.5). Since we use the primitive cell of the oxygen sub-lattice to represent a fluorite metal oxide structure half the (0.5,0.5,0.5) sites will be occupied with metal ions and half will be marked as empty. Oxygen-vacancy diffusion was described with six processes describing nearest-neighbor hops along the Cartesian axes. The system was run with two different settings. In the first setting oxygen diffusion was modeled with equal rates for all six processes. In the second setting the system was run with a custom rate calculator that modified the process rates based on the local distribution of dopants within three shells of primitive cells (258 neighbor sites) from the vacancy site of interest.

Serial performance
All serial performance tests were carried out on a laptop computer running Ubuntu-12.04 LTS, with an Intel i7-3517U 1.90 GHz CPU with 4 MB cache, and with 8 GB RAM. Only one of the CPU's cores was used.
In figure 3 we show the scaling with system size for the sim-  Figure 4 shows the scaling with system size for the fluorite metal oxide oxygen-vacancy diffusion system in the first setting. Both the initialization time and memory usage (4a and c) scales linearly with system size as expected. Contrary to the constant scaling for the A-B-C model, the time spent for each step in the loop (figure 4b) also grows linearly with system size. This linear behavior arises since the system size is scaled homogeneously so that the number of active sites (sites where the processes can be applied) grows linearly with the total size of the system. When a site is matched against a process the list of available sites for that process must be searched through for the site of interest to determine if the site should be added to (or removed from) this list. When the number of available sites grows with system size the list to search through grows accordingly, which gives rise to the observed linear scaling behavior.
The situation is quite different for the fluorite metal oxide oxygen-vacancy diffusion system in the second setting ( figure  5). While the initialization time and memory usage (figure 5a and c) still scales linearly, the time per step in the loop (figure 5b) is close to constant, with only a slight increase with system size. This behavior is expected and can be explained by noting the difference in scale on the y-axis in figures 4b and 5b. The linear component from figure 4b is still present in figure 5b, but the time spent in the loop is completely dominated by the rate evaluation resulting in close to constant O(1) scaling.
We expect this type of roughly constant O(1) scaling for most realistic applications even when custom rates are not used, since the matching step is typically far more time consuming than the search through the availability vectors of the processes.

System size and simulation time
A peculiarity when it comes to the scaling behavior of any KMC simulation implementing the VSSM algorithm, which at first can be easy to overlook and therefore well worth noting, is how the simulation time (as opposed to the wall-clock time measured in the timings above) varies with the system size. We recall that the simulation time is propagated according to the expression (3) where the sum of the rates in the system appears in the denominator. The length of each time-step in the simulation is thus inversely proportional to the total rate in the system. This means that doubling the size of the system (and the number of available sites) and running the simulation for the same number of elementary steps will only propagate the simulation time half as far.

Parallel performance
All parallel performance tests were run on a CentOS Linuxbased cluster with HP Proliant SL2x170z nodes, each equipped with two quad-core Intel Xeon E5520 CPUs and 36 GB RAM, with infiniband interconnect. One core was used for each MPI process.
As discussed in section 3.4 above the parallel performance is expected to depend strongly on the amount of work needed during the custom rate calculations. To run the same test system for different times spend in the rate calculation function a sleep statement was introduced. Three different cases were considered, waiting an additional 1, 10 and 100 milliseconds respectively each time the rate calculation function was called. Figure 6 shows the scaling calculated as the total time running on one MPI process over the total time running in parallel. As expected, the more time spent in the custom rate calculation the better the scaling. When only 1 millisecond is added to the rate calculation the scaling flattens out above 16 processes, but with 100 additional milliseconds in the rate calculation the scaling is excellent all the way up to the 128 processes.
Clearly, the heavier the custom rate calculations the more is to gain from running in parallel, but the number of sites to rematch at each step will also influence the parallel performance. The test system used in figure 6 is three-dimensional with two basis sites for each primitive unit cell and a large cutoff which resulted in more than 500 process-site pairs to re-evaluate at each step. For a good load balance it is necessary that the number of process-site pairs to re-match is either much larger than, or an exact multiple of, the number of MPI processes. Running efficiently in parallel will therefore be limited to a few MPI processes for systems with few neighbors (with low dimensionality) and short cutoff.

Caching of calculated rates
In those cases where the on-the-fly rates are kept constant over time, i.e., when no time-evolution is allowed in the custom rate expressions, it is possible to use caching techniques to save calculated rates for later retrieval when needed. The efficiency of such a caching scheme is highly dependent on the details of the simulated system and on the caching algorithm used, but it can in some cases provide a significant speedup. For one of the simpler systems investigated, the two-dimensional Ising spin model (see section 3.3.2), the custom rate implementation can be sped up almost to the level of the fixed rates implementation using a prototype caching scheme where all the calculated custom rates are saved.
For systems with larger number of possible lateral interactions saving all calculated rates is not feasible due to memory limitations. A caching scheme with several internal hash tables with fixed maximum size, where the oldest table filled is emptied when the latest table reaches its maximum size limit would resolve the memory issue. No cashing mechanism is implemented in the present version of the code but will be further investigated for future releases, however a prototype implementation is available in a development branch of the KMCLib git repository. The structure of the custom rate framework meanwhile makes it easy for the user to implement her own cashing mechanism.

Concluding remarks
We have presented a general framework for lattice KMC simulations designed to easily be customized, to facilitate implementation of custom KMC models without having to re-do all the work of implementing and testing the simulation engine for each new set of problems to investigate. We have made our code freely available since we believe this will be of great benefit for the larger research community, by saving time and letting researchers focus their effort on the design of new models to solve challenging physical problems, rather than spending time re-implementing a well-known method.
This paper has focused mainly on technical details and on the modifications we have made to the standard lattice KMC algorithm. We have demonstrated the scaling behavior in serial and shown how the parallel performance depends on the simulated system. The code itself is separately documented from a users perspective, and includes several usage examples in the functionality tests. We have therefore not given a full usage example here, but trust that the reader will find the code manual and examples enough to get started. To run the tests and set up new calculations require some knowledge of Python programming but following the examples and the documentation this barrier should be easy to overcome.
The modular plugin design of the program with a welldefined API is ideal for sharing programming work with the community. Custom rate calculators, analysis modules and scripts for setting up geometries and processes and to visualize results can be written and published separately, and we hope to see many such contributions from the research community, to include in future releases.