A kinetic Monte Carlo model with Fast Multipole Method

This program was written by Ian Thompson (IRT) and uses a fast multipole electrostatic algorithm for the energy evaluations (as decribed in paper by Will et al.). FMM written by Will Saunders and derived from the PPMD framework

Installation

Make a Python virtualenv then install the following,

Install PPMD:

https://github.com/ppmd/ppmd

Install coulomb_kmc (an FMM-KMC implementation):

https://github.com/ppmd/coulomb_kmc

Install this package:

pip install --upgrade --no-cache-dir git+https://gitlab.com/ABW_bath_group/kmc_fmm@master

Running

After pip installing this package, the kmc-fmm executable can be ran with

kmc-fmm input.json labelname

where input.json and labelname are user inputs. An example input configuration is given in examples/aNPD.

Variable	Description
Ns	The number of hopping sites in the unit cell
Nb	The number of neighbours a site has
MAX_Nb	The maximum number of neighbours a site can have
Np	The number of particles in the simulation
L_U	The lengths of the unit cell
L_S	The lengths of the simulation cell
R_U	The number of repeat units of the simulation cell

General principle

A kMC simulation that uses Marcus rates to simulate the motion of charge carriers in a bulk (3D periodic) organic semiconductor system. The system must be charge neutral and the total simulation cell must be cubic. To achieve this the simulation supports regular tiling of a primitive (non-cubic) unit cell to create an overall cubic system. The output is a params file and a single trajectory file that is placed into a date stamped folder that is itself placed into a folder in saved_states/[voltage]V (by default).

There is a clear separation between recording the change in the system for the user (e.g. recording a trajectory) and modifying the system to reflect the event that has been picked (e.g. moving a particle or creating one). The fomrer is up to the user and can be done however as long as dats are only read and not modified. The latter must be done to the PPMD state and should use wrapper functions.

##A note on the backend

The guts are mostly done but to add events the user needs to modify a few kernels: this is due to the role of the allowed_move_mask. This is a Np x MAX_Nb dat that for each particle determines whether a move to all possible neighbours is allowed. The value of the mask determines the event type that occurs.

The copy_prop_kernel masks off neighbours that don't exist (where Nb < MAX_Nb) and sets the rest to 1. This happens at the start of every kMC step.

Then the exclude kernel sets the rest of the mask appropriately. It checks every proposed position against all other particles, the occupancy of the site determines what can happend. At the moment if neighbouring site n is empty we set the mask mask.i[n] =1, this is a hopping event. If it is occupied by a blocking particle then mask.i[n] = 0, this is a blocked event. If it is occupied by a particle that has an interaction with particle i then we set mask.i[n] = -1.

Mask = 0: Null event, nothing happens, no further calculations.
Mask > 0: An event that you want to calculate the difference in electrostatic energy for.
Mask < 0: An event that you do NOT want to calculate the difference in electrostatic energy for e.g. an event where two charges would be on top of one another and so the energy is infinite.

The user can decide what number is what type of event as long as they follow the above rule.

Finally the rate calculation kernel uses the mask to set which rate function should be used. The user needs to specify the rate for each type of event, due to the near limitless zoo of possibilities I have avoided providing presets.

Then in the main program after an event is picked the user can query the allowed_move_mask to find out what type of event has been picked, make the relevant changes to the system and record data.

Physical units

Lengths are in nanometres
Times are in picoseconds
Energies are in electron Volts

Inputs for the simulation

For each site in the unit cell the simulation needs to know:

The type of the hopping site. Ns array
The position of the hopping site in XYZ, in the unit cell in nanometres. Ns x 3 array
The number of neighbours the site has. Ns array
The indices of the neighbours in the unit cell (any values needed to pad the array to MAX_Nb can be any value). Ns x MAX_Nb array
The strength of the couplings between the site and its neighbours (this is expressed as a transfer integral but could be any quantity that is defined per pair). Again padding can be any value. Ns x MAX_Nb array
The site energies of each site. For a single energy landscape (e.g. HOMO or LUMO) this is a Ns array
The offsets from the site to each if its neighbours in X, Y and Z. These distances should NOT be periodically wrapped as particles can cross from one unit cell to another (i.e. use the short distance), the simulation takes care of the wrapping. Ns x MAX_Nb x 3 array. (Optional - can be calculated from other inputs using get_real_space_offsets)
The size of the unit cell in nanometres. 3 array
The number of repeats of the unit cell in each direction (can be 1, 1, 1). 3 array
The maximum distance between any two sites that can be occupied, far_cutoff.

These can be read in from a file anyway the user wants, as long as the numpy array exists of the right shape it can be converted into a PPMD dat.

##The format of the input files used by IRT are defined as:

###A JSON input file, the name is a command line argument

The JSON file contains run parameters such as the number of steps, the temperature, relative permittivity etc. The entries have nice, straighforward names and are parsed as a Python dictionary at the top of the main program. Entries can be added easily but beware of backwards compatibility.

geometry.dat

Each line is a site formatted in 5 columns as:

| geometry.dat |

By convention electron transporting material sites are type 1, hole transporting material sites are type 2, materials that can transport both are 4, positive dopants are type 101, negative dopants are type 102.

neighbours.dat

Each site has 2 lines of data: the first is a list of the neighbour indices, the second is the distance between the site and neighbour, formatted as:

| neighbours.dat |

energetics.dat

Each site has 2 lines of data: the first is a set of energies, the second is the couplings between the site and each neighbour, formatted as:

| energetics.dat |

The absorption coefficient is unused in the kMC version given. For other uses than organic semiconductors HOMO and LUMO can be changed.

device_size.dat

This can be in one of two formats, either:

Just the unit cell size in a column, e.g.:

Unit_cell_LX   
Unit_cell_LY   
Unit_cell_LZ

Or as the first column in file e.g.:

Unit_cell_LX Sim_cell_LX    
Unit_cell_LY Sim_cell_LY   
Unit_cell_LZ Sim_cell_LZ

pair_energies.dat (optional)

Each site has 2 lines of data: the first is a list of the difference in energies between the site and neighbours, the second is the reorganisation energy between the site and each neighbour, formatted as:

| pair_energies.dat |

This is optional as you can easily calculate the difference in site energies from the energetics.dat file and the simulation uses a constant value for the reorganisation energy anyway.

Output

The trajectory file contains a line per event and is formatted as:

| event_trajectory.traj |

For a hopping event there is no Particle2 but this is designed to allow recombination and fission. By convention electrons are type 1, holes are type 2. The duration is in picoseconds and the site indices are calculated according to: site_idx + (repeat_unit_idx * Ns), where site_idx is in the unit cell, so that each site in the simulation cell has a unique index.

The trajectory file can be read by the kMC_trajectory_reader program.

The results of the search are