Example of Kinetic Monte Carlo simulation of reaction in Langmuir binary gas
Download latest version here: test_kmc_binary_gas.py.
import numpy as np
from susmost.mc import metropolis
from susmost.mc.metropolis import kmc_step, kmc_rate
from susmost import make_single_point_zmatrix, make_square_plane_cell, \
normal_lattice_task, nearest_int, SiteStateType, register_property, sitestate, mc
from time import time
from math import exp
register_property('coverage_A','==', None, sitestate.join_cover) # counter for molecule A coverage
register_property('coverage_B','==', None, sitestate.join_cover) # counter for molecule B coverage
sitestate.SPMI_db['coverage'].merge_func = None # prevent merging of occupied and empty states due to equal lateral interactions
def binary_gas_KMC(rnd_seed, muA, muB, eAA, eBB, diff_Ea, react_Ea, Temperature):
"""
Do KMC simulation for a binary gas with symmetric lateral interactions between nearest neighbors.
Parameters:
:rnd_seed: random seed
:muA: chemical potential of type A molecules
:muB: chemical potential of type B molecules
:eAA: lateral interaction energy between type A molecules
:eBB: lateral interaction energy between type B molecules
:diff_Ea: diffusion barrier for A and B molecules
:react_Ea: activation energy for A + B reaction
:Temperature: temperature in k_B units
Returns:
Tuple of three elements: (A coverage, B coverage, number of reactions)
"""
INF_E = 1E6 # prohibitive energy level
#k_B = 8.3144598E-3 # kJ/(mol*K)
#k_B = 0.000086 # eV
k_B = 1.0
beta = 1.0 / (k_B * Temperature)
a = 1.0 # lattice parameter
r_cut_off= 1.0 + .0001
dop_int = {('molecule_A','molecule_B'):0.0,('molecule_A','molecule_A'):eAA,('molecule_B','molecule_B'):eBB}
lateral_interactions = lambda cc1,cc2: nearest_int(cc1,cc2,dop_int=dop_int)
#cell,atoms = make_struct.make_hexagonal_plane_cell(a)
cell,atoms = make_square_plane_cell(a)
empty_state = SiteStateType('Empty',make_single_point_zmatrix(), 0.0, ['H'], coverage=0., coverage_A=0., coverage_B=0.)
molA_state = SiteStateType('molecule_A', make_single_point_zmatrix(), -1., ['N'], coverage=1., coverage_A=1., coverage_B=0.)
molB_state = SiteStateType('molecule_B', make_single_point_zmatrix(), 0.0, ['O'], coverage=1., coverage_A=0., coverage_B=1.)
lt = normal_lattice_task([empty_state, molA_state, molB_state], cell, atoms, lateral_interactions, r_cut_off, INF_E)
assert lt.states_count == 3, lt
lt.set_ads_energy('molecule_A', muA)
lt.set_ads_energy('molecule_B', muB)
empty_idx = lt.find_state_index('name','Empty')
A_idx = lt.find_state_index('name','molecule_A')
B_idx = lt.find_state_index('name','molecule_B')
m = mc.make_metropolis(lt, 24, Temperature, k_B, seed=int(rnd_seed))
ads_des_rate = kmc_rate(1, 0, 1)
diff_rate = kmc_rate(1, diff_Ea, 1)
react_rate = kmc_rate(1, react_Ea, 0)
all_rates = []
all_rates.append(ads_des_rate.calc(beta, muA)) # dH=muA
all_rates.append(ads_des_rate.calc(beta, -muA)) # dH=muA
all_rates.append(ads_des_rate.calc(beta, muB)) # dH=muB
all_rates.append(ads_des_rate.calc(beta, -muB)) # dH=muB
all_rates.append(diff_rate.calc(beta, 0)) # dH=0
all_rates.append(react_rate.calc(beta, 0)) # dH=0
r0 = max(all_rates)
#print(all_rates, beta)
m.add_singe_site_step(empty_idx, kmc_step(ads_des_rate, A_idx, -1, 0)) # A adsorption
m.add_singe_site_step(empty_idx, kmc_step(ads_des_rate, B_idx, -1, 0)) # B adsorption
m.add_singe_site_step(A_idx, kmc_step(ads_des_rate, empty_idx, -1, 0)) # A desorption
m.add_singe_site_step(B_idx, kmc_step(ads_des_rate, empty_idx, -1, 0)) # B desorption
for ei in m.regr.nn_edges_range(): # nearest neighbors in all directions
m.add_pair_step(empty_idx, ei, A_idx, kmc_step(diff_rate, A_idx, empty_idx, 0)) # A diffusion
m.add_pair_step(A_idx, ei, empty_idx, kmc_step(diff_rate, empty_idx, A_idx, 0)) # A diffusion
m.add_pair_step(empty_idx, ei, B_idx, kmc_step(diff_rate, B_idx, empty_idx, 0)) # B diffusion
m.add_pair_step(B_idx, ei, empty_idx, kmc_step(diff_rate, empty_idx, B_idx, 0)) # B diffusion
m.add_pair_step(A_idx, ei, B_idx, kmc_step(react_rate, empty_idx, empty_idx, 1)) # A + B reaction
m.add_pair_step(B_idx, ei, A_idx, kmc_step(react_rate, empty_idx, empty_idx, 1)) # B + A reaction (same as above)
m.params_log_period = 100 # log parameters each 100 steps
m.setup_kmc(r0)
mc.run(m, 100, 100000) # do 100 iterations of 100000 steps each, saving surface snapshots in traj*.xyz file
#mc.stat_digest(m) # print detailed statistics about logged parameters
assert m.kmc_max_r0 <= r0, f"Selected r0 = {r0} turned to be exceeded by rate = {m.kmc_max_r0}"
A_cov_idx = m.param_idx('coverage_A')
B_cov_idx = m.param_idx('coverage_B')
means = np.mean(m.params_log, axis=0) # mean values of all parameters
return means[A_cov_idx], means[B_cov_idx], m.kmc_counters[1].accept
Ts = np.arange(.1, 3+1e-9, 0.1)
R = []
for T in Ts:
res = binary_gas_KMC(time(), -1., 0.0, 0.0, 1e-9, 0.0, 0.0, T) # 1e-9 - hack to prevent dropping of all lattice edges :-(, TODO: fix it
print (T, *res)
R.append(res)
print ("\nFinal result:\nT covA covB Rate")
for T, res in zip(Ts, R):
print (T, *res)
Result plot:
