Example of Kinetic Monte Carlo simulation of ZGB model

Download latest version here: test_zgb.py and ZGB model

Original ZGB paper: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.56.2553

from susmost import load_lattice_task
from susmost import mc
from susmost.mc.metropolis import kmc_step, kmc_rate
import numpy as np

lt = load_lattice_task('models/ZGB')
lt.set_property('O_cnt', {"O":1})       # counter of O coverage
lt.set_property('C_cnt', {"C":1})       # counter of CO coverage

assert lt.states_count == 3                     # 3 lattice site states
m = mc.make_metropolis(lt, 20, [1.0], 1.0) # lattice size is 20x20; temperature and k_B are irrelevant here so set to 1.0


empty_idx = 0
C_idx = lt.find_state_index('name','C')
O_idx = lt.find_state_index('name','O')

yCO = 0.54 # critical values y1 = 0.389, y2 = 0.525
rCO = 4*yCO/(1+3*yCO)  # reaction rate is specified separately for each oxygen dimer orientation, so it should be adjsuted
m.add_singe_site_step(empty_idx, kmc_step(kmc_rate(rCO, 0, 0), C_idx, -1, 0))   # CO adsorption, 0-th counter

for ei in m.regr.nn_edges_range():      # nearest neighbors in all directions
        m.add_pair_step(empty_idx, ei, empty_idx, kmc_step(kmc_rate(1.0 - rCO,0,0), O_idx, O_idx, 1))  # O2 adsorption, counter index == 1
        m.add_pair_step(O_idx, ei, C_idx, kmc_step(kmc_rate(1000.0,0,0), empty_idx, empty_idx, 2))              # O + CO reaction, 2nd counter
        m.add_pair_step(C_idx, ei, O_idx, kmc_step(kmc_rate(1000.0,0,0), empty_idx, empty_idx, 2))              # CO + O reaction, 2nd counter

m.setup_kmc(1000.0)     # r0 == 1000.0

m.params_log_period = 100
O_cnt_idx = m.param_idx('O_cnt')
C_cnt_idx = m.param_idx('C_cnt')

for i in range(1000):
        m.run(1000000)
        m.save_as_xyz("zgb.xyz", comment='E')
        # print mean O coverage, CO coverage, acceptance rate, 2nd counter (number of CO + O reactions
        print (i, np.mean(m.params_log[:,O_cnt_idx]), np.mean(m.params_log[:,C_cnt_idx]), np.mean(m.params_log[:,-4]), m.kmc_counters[2].accept)

Visualization of zgb.xyz with y_CO = 0.530: