# This file was automatically generated by build-ipynb_FILES.sh from
# the CoveragesAndReactionRate.ipynb file. All changes to this file will be lost

# The scale disparity issue is common in Kinetic Monte Carlo simulations. It emerges
# as a result of the fact that some fundamental events or groups of them typically
# occur at vastly different time scales; in other words, their rate constants can
# span multiple orders of magnitude. In heterogeneous catalysis, there are typically
# two groups: 1) very fast events that correspond to the species' surface diffusions
# and 2) slow reactions that change their chemical identity. The latter group of events
# is usually the one of interest because it allows the evaluation of the material's
# catalytic activity. In contrast, the species' surface diffusion does not contribute
# significantly to the net evolution of the slow reactions. But, as becomes the most
# frequent step, it also becomes the limiting factor of the simulation progress,
# considerably increasing the computational cost. This tutorial shows how to speed up
# the calculation by several orders of magnitude without sacrificing precision by
# automatically detecting and scaling the rate constants of fast reactions.
#
# We will focus on the net reaction $\text{CO}+\frac{1}{2}\text{O}_2\longrightarrow \text{CO}_2$
# that takes place at a catalyst's surface and whose reaction mechanism is described by
# the Langmuir-Hinshelwood model. Because this model has four very fast processes
# ($CO$ and $O_2$ adsorption, and $O*$ and $CO*$ diffusion) and one slow process
# ($CO$ oxidation), it is an ideal prototype for demonstrating the benefits of the
# automated rescaling of rate constants technique. Our ultimate goal is to investigate
# how altering the relative percentage of the gas reactants $CO$ and $O_2$ (at a specific
# temperature and pressure) affect the rate of $CO_2$ production under steady state
# conditions. This example is inspired on Zacros tutorial
# [What's KMC All About and Why Bother](https://zacros.org/tutorials/12-about-kinetic-monte-carlo?showall=1).

# Let's start! The first step is to import all packages we need:

import multiprocessing
import numpy
import scm.plams
import scm.pyzacros as pz
import scm.pyzacros.models


# Then, we initialize the **pyZacros** environment.

scm.pyzacros.init()


# Notice this command created the directory where all **Zacros** input and output files
# will be stored if they are needed for future reference (``plams_workdir`` by default).
# Typically, the user doesn't need to use these files.

# This calculation necessitates a significant computational effort. On a typical laptop,
# it should take around 20 min to complete. So, in order to speed things up, we'll
# use the ``plams.JobRunner`` class to run as many parallel instances as possible. In this
# case, we choose to use the maximum number of simultaneous processes (``maxjobs``) equal
# to the number of processors in the machine.

maxjobs = multiprocessing.cpu_count()
scm.plams.config.default_jobrunner = scm.plams.JobRunner(parallel=True, maxjobs=maxjobs)
print("Running up to {} jobs in parallel simultaneously".format(maxjobs))


# First, we initialize our Langmuir-Hinshelwood model, which by luck is available as a
# predefined model in pyZacros,

lh = pz.models.LangmuirHinshelwood()


# Then, we must set up a ``ZacrosParametersScanJob`` calculation, which will allow
# us to scan the molar fraction of $CO$ as a parameter. However, this calculation
# requires the definition of a ``ZacrosSteadyStateJob``, that in turns requires a
# ``ZacrosJob``. So, We will go through them one at a time:

# **1. Setting up the ZacrosJob**
#
# For ``ZacrosJob``, all parameters are set using a ``Setting`` object. To begin,
# we define the physical parameters: ``temperature`` (in K), and ``pressure``
# (in bar). The calculation parameters are then set: ``species numbers`` (in s)
# determines how frequently information about the number of gas and surface species
# will be stored, ``max time`` (in s) specifies the maximum allowed simulated time,
# and "random seed" specifies the random seed to make the calculation precisely
# reproducible. Keep in mind that ``max time`` defines the calculation's stopping
# criterion, and it is the parameter that will be controlled later to achieve the
# steady-state configuration. Finally, we create the ``ZacrosJob``, which uses the
# parameters we just defined as well as the Langmuir-Hinshelwood model's lattice,
# mechanism, and cluster expansion. Notice we do not run this job, we use it as a
# reference for the steady-state calculation described below.

z_sett = pz.Settings()
z_sett.temperature = 500.0
z_sett.pressure = 1.000
z_sett.species_numbers = ("time", 1.0e-5)
z_sett.max_time = 100 * 1.0e-5
z_sett.random_seed = 1609

z_job = pz.ZacrosJob(
    settings=z_sett, lattice=lh.lattice, mechanism=lh.mechanism, cluster_expansion=lh.cluster_expansion
)


# **2. Setting up the ZacrosSteadyStateJob**
#
# We also need to create a ``Setting`` object for ``ZacrosJob`` There, we ask for a
# steady-state configuration using a TOFs calculation with a 95% confidence level
# (``turnover frequency.confidence``), using four replicas to speed up the calculation
# (``turnover frequency.nreplicas``); for more information, see example
# **Ziff-Gulari-Barshad model: Steady State Conditions**. Then, we ask for the
# rate constants to be automatically  scaled (``scaling.enabled``) using an inspection
# time of 0.0006 s ("scaling.max time"). In a nutshell, the scaling algorithm uses
# this maximum time and the original rate constants to execute a probe simulation.
# From there, the occurrence rates for each reaction are determined and scaled
# appropriately. Only reactions that are proven to be quasi-equilibrated are modified.
# The actual simulation is then started from the beginning using the new reaction
# rates following this. The ``ZacrosSteadyStateJob.Parameters`` object allows to set the grid of maximum
# times to explore in order to reach the steady state (``ss_params``). Finally, we
# create ``ZacrosSteadyStateJob``, which references the ``ZacrosJob`` defined above
# (``z_job``) as well as the ``Settings`` object and parameters we just defined:

ss_sett = pz.Settings()
ss_sett.turnover_frequency.confidence = 0.95
ss_sett.turnover_frequency.nreplicas = 4
ss_sett.scaling.enabled = "T"
ss_sett.scaling.max_time = 60 * 1e-5

ss_params = pz.ZacrosSteadyStateJob.Parameters()
ss_params.add("max_time", "restart.max_time", 2 * z_sett.max_time * (numpy.arange(10) + 1) ** 2)

ss_job = pz.ZacrosSteadyStateJob(settings=ss_sett, reference=z_job, parameters=ss_params)


# **3. Setting up the ZacrosParametersScanJob**
#
# Although the ``ZacrosParametersScanJob`` does not require a ``Setting`` object,
# it does require a ``ZacrosSteadyStateJob.Parameters`` object to specify which
# parameters must be modified systematically. In this instance, all we need is a
# dependent parameter, the $O_2$ molar fraction ``x_O2``, and an independent
# parameter, the $CO$ molar fraction ``x_CO``, which ranges from 0.05 to 0.95.
# Keep in mind that the condition ``x_CO+x_O2=1`` must be met. These molar fractions
# will be used internally to replace ``molar fraction.CO`` and ``molar fraction.O2``
# in the Zacros input files. Then, using the ``ZacrosSteadyStateJob`` defined
# earlier (``ss job``) and the parameters we just defined (``ps params``), we
# create the ``ZacrosParametersScanJob``:

ps_params = pz.ZacrosParametersScanJob.Parameters()
ps_params.add("x_CO", "molar_fraction.CO", numpy.linspace(0.05, 0.95, 11))
ps_params.add("x_O2", "molar_fraction.O2", lambda params: 1.0 - params["x_CO"])

ps_job = pz.ZacrosParametersScanJob(reference=ss_job, parameters=ps_params)


# The parameters scan calculation setup is ready. Therefore, we can start it
# by invoking the function ``run()``, which will provide access to the results via the
# ``results`` variable after it has been completed. The sentence involving the method
# ``ok()``, verifies that the calculation was successfully executed, and waits
# for the completion of every executed thread. Go and grab a cup of coffee, this
# step will take around 15 mins!

results = ps_job.run()

if not ps_job.ok():
    print("Something went wrong!")


# If the execution got up to this point, everything worked as expected. Hooray!
#
# Finally, in the following lines, we just nicely print the results in a table. See
# the API documentation to learn more about how the ``results`` object is structured,
# and the available methods. In this case, we use the ``turnover_frequency()`` and
# ``average_coverage()`` methods to get the TOF for the gas species and average coverage
# for the surface species, respectively. Regarding the latter one, we use the last ten
# steps in the simulation to calculate the average coverages.

x_CO = []
ac_O = []
ac_CO = []
TOF_CO2 = []

results_dict = results.turnover_frequency()
results_dict = results.average_coverage(last=10, update=results_dict)

for i in range(len(results_dict)):
    x_CO.append(results_dict[i]["x_CO"])
    ac_O.append(results_dict[i]["average_coverage"]["O*"])
    ac_CO.append(results_dict[i]["average_coverage"]["CO*"])
    TOF_CO2.append(results_dict[i]["turnover_frequency"]["CO2"])

print("------------------------------------------------")
print("%4s" % "cond", "%8s" % "x_CO", "%10s" % "ac_O", "%10s" % "ac_CO", "%12s" % "TOF_CO2")
print("------------------------------------------------")
for i in range(len(x_CO)):
    print("%4d" % i, "%8.2f" % x_CO[i], "%10.6f" % ac_O[i], "%10.6f" % ac_CO[i], "%12.6f" % TOF_CO2[i])


# The results table above demonstrates that when $CO$ coverage rises, the net $CO$
# oxidation reaction tends to progress more quickly until it reaches a maximum of
# about 0.7. At this point, it begins to decline. Notice that the maximal $CO_2$
# generation occurs simultaneously as the $CO*$ coverage reaches parity with the
# $O*$ coverage. This makes intuitive sense since, at that point, we maximize the
# likelihood of discovering an $O*$ and a $CO*$ that are close enough to react in
# accordance with the proper reaction's stoichiometry. The number of $O*$ or $CO*$
# observers rises as we move away from that critical point, decreasing the likelihood
# that the $CO$ oxidation process could take place.

# Now, we can close the pyZacros environment:

scm.pyzacros.finish()