notebooks/care_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CARE demo\n",
    "\n",
    "This notebook covers the essential steps to generate and manipulate the chemical reaction networks (CRNs) generated with CARE. All the steps shown here are included in the main script `care.scripts.care_run`, which is the recommended way to generate the CRN.\n",
    "Install Ipykernel to run the notebook!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Chemical Reaction Network: Blueprint\n",
    "\n",
    "To generate the CRN blueprint, i.e. the listing of the potential species and elementary reactions in your CRN, use the function `gen_blueprint`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from care import gen_blueprint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inputs\n",
    "\n",
    "*ncc* = Network Carbon Cutoff, i.e. the maximum number of C atoms that a species in the CRN can possess  <br>\n",
    "*noc* = Network Oxygen Cutoff, i.e. the maximum number of C atoms that a species in the CRN can possess  <br>\n",
    "*cyclic* = Whether species with rings are should be included or not (True/False)  <br>\n",
    "*additional_rxns* = Whether to include intramolecular rearrengments (True/False)  <br>\n",
    "*electro* = Whether to generate a CRN in electrochemical conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ncc, noc = 2, 2\n",
    "cyclic = False\n",
    "rearr_rxns = True\n",
    "electro = False\n",
    "\n",
    "inters, rxns = gen_blueprint(ncc=ncc, \n",
    "                             noc=noc, \n",
    "                             cyclic=cyclic, \n",
    "                             additional_rxns=rearr_rxns, \n",
    "                             electro=electro)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of the function is a dict of `Intermediate` objects indexed by their InChIKey (plus a phase-identifier character, '*' for adsorbed, 'g' for gas-phase), and a list of `ElementaryReaction` instances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(type(inters)) \n",
    "print(type(inters[list(inters.keys())[0]]))\n",
    "print(type(rxns))\n",
    "print(type(rxns[-1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ElementaryReaction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = rxns[10]\n",
    "print(r)\n",
    "print(type(r))\n",
    "print(r.repr_hr)  # human-readable text representation\n",
    "print(r.r_type)  # bond-breaking type\n",
    "print(r.components)  # components of the reaction\n",
    "print(r.stoic)  # Soichiometry coefficients in the reaction\n",
    "print(r.e_rxn)  # reaction energy\n",
    "print(r.e_act)  # reaction activation energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observations:\n",
    "- species participating in the elementary reactions are displayed as '[n]InChIKeyP(formula)', with 'n' being the absolute stoichiometric coefficient, formula being the brute formula provided by ASE, and 'P' after the InChIKey being a phase-identifier character, '*' for adsorbed species and 'g' for gaseous species. <br>\n",
    "- ``e_rxn`` and ``e_act`` are ``None`` as these reaction properties still have to be evaluated\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Intermediate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = inters['WSFSSNUMVMOOMR-UHFFFAOYSA-N*']\n",
    "print(a)\n",
    "print(type(a))\n",
    "print(a.phase)\n",
    "print(a.smiles)\n",
    "print(a.ref_energy())\n",
    "print(a.is_closed_shell())\n",
    "print(a['C'])  # number of carbon atoms\n",
    "print(a['H'])  # number of hydrogen atoms\n",
    "print(a['O'])  # number of oxygen atoms\n",
    "print(a.ads_configs)  # adsorption configurations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observations:\n",
    "- ``Intermediate.ads_configs`` is a dictionary where the adsorption configurations will be stored. As the blueprint considers the network generation in gas-phase, is empty for now.\n",
    "- the ``Intermediate.code`` contains a phase-identifier character between the InChIKey and the brute formula, being '*' for adsorbed species and 'g' for gaseous species. Changing ``WSFSSNUMVMOOMR-UHFFFAOYSA-N*`` to ``WSFSSNUMVMOOMR-UHFFFAOYSA-Ng`` returns a copy of the same object, they will represent different species onces these are evaluated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Chemical Reaction Network: Energy Evaluation\n",
    "\n",
    "Now we define the surface under study and evaluate the energies of intermediates and reactions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load evaluator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from care.evaluators import load_surface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metal = 'Ru'\n",
    "facet = '0001'  # hkl notation\n",
    "\n",
    "surface = load_surface(metal, facet)\n",
    "\n",
    "print(surface)\n",
    "print(type(surface))\n",
    "print(surface.slab)  # ASE Atoms object\n",
    "print(surface.active_sites)  # Obtained with ACAT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Intermediate energy evaluator\n",
    "\n",
    "In our case, this is GAME-Net-UQ, a graph neural network trained on DFT data. CARE allows you to implement your own model via two defined interfaces, one for intermediates and the other for reaction properties (see next cells)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from care.evaluators.gamenet_uq import GameNetUQInter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gnn = GameNetUQInter(surface)\n",
    "\n",
    "print(gnn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine the domain of applicability of the evaluator in terms of potential adsorbate composition and surface characteristics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gnn.adsorbate_domain()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gnn.surface_domain()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate intermediates (sequential)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have 171 intermediates to evaluate. We will show here the sequential evaluation to give an idea. If you want to get the same output faster, skip this section and run the parallel evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inters_evaluated = {}\n",
    "for k, inter in inters.items():\n",
    "    print(inter)\n",
    "    inters_evaluated[k] = gnn.eval(inter)\n",
    "    print(inters_evaluated[k].ads_configs)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inters_evaluated['LHGAACIDYUKUTF-UHFFFAOYSA-N*'].ads_configs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate intermediates (parallel)\n",
    "\n",
    "Parallel execution across multiple CPU cores can speed up the evaluation process. In our case, the parallel execution took 1min 51s, while the sequential 20 minutes (10x speed up with 24 cores)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import multiprocessing as mp\n",
    "import resource\n",
    "\n",
    "_, hard = resource.getrlimit(resource.RLIMIT_NOFILE)\n",
    "resource.setrlimit(resource.RLIMIT_NOFILE, (hard, hard))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tasks = [(gnn, inter) for inter in inters.values()]\n",
    "print(f'{len(tasks)} intermediates to evaluate')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mp_func(gnn, inter):\n",
    "    print(inter.code)\n",
    "    gnn.eval(inter)\n",
    "    return inter\n",
    "\n",
    "with mp.Pool(mp.cpu_count()) as pool:\n",
    "    results = pool.starmap(mp_func, tasks)\n",
    "\n",
    "inters_evaluated = {i.code: i for i in results}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get random key and print the adsorption configurations\n",
    "key = list(inters_evaluated.keys())[18]\n",
    "inters_evaluated[key].ads_configs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate reaction properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from care.evaluators.gamenet_uq import GameNetUQRxn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gnn_rxn = GameNetUQRxn(inters_evaluated)\n",
    "print(gnn_rxn)\n",
    "\n",
    "print(f'{len(rxns)} reactions to evaluate')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rxns_evaluated = []\n",
    "for rxn in rxns:\n",
    "    print(rxn)\n",
    "    gnn_rxn.eval(rxn)\n",
    "\n",
    "rxns_evaluated = sorted(rxns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "idx = random.randint(0, len(rxns_evaluated))\n",
    "\n",
    "step_example = rxns_evaluated[idx]\n",
    "print(step_example, step_example.r_type)\n",
    "print(step_example.e_is)\n",
    "print(step_example.e_ts)\n",
    "print(step_example.e_fs)\n",
    "print(step_example.e_rxn)\n",
    "print(step_example.e_act)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore the CRN\n",
    "\n",
    "Everything is now evaluated. Let's assemble the reaction network.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from care import ReactionNetwork"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crn = ReactionNetwork(intermediates=inters_evaluated, \n",
    "                      reactions=rxns_evaluated,\n",
    "                      surface=surface,\n",
    "                      ncc=ncc,\n",
    "                        noc=noc,\n",
    "                      type='thermal')\n",
    "\n",
    "print(crn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get top-10 hubs in the network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crn.get_hubs(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iterate through elementary steps and show reaction energy and activation barrier. Each property is displayed as a tuple where the first item is the mean value and the second the standard deviation (uncertainty) of the prediction. Elementary reactions having same values for both reaction energy and barrier correspond to barrierless endothermic steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crn.get_reaction_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize elementary steps by providing the reaction index. Remember that the uncertainty of the properties is a relative measure of the model uncertainty across multiple predictions!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crn.visualize_reaction(23, show_uncertainty=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get number of closed-shell molecules in the network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crn.num_closed_shell_mols"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize intermediates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crn.visualize_intermediate('VRLIPUYDFBXWCH-UHFFFAOYSA-N*')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize graph representation (only when predictions are made with GAME-Net-UQ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from care.evaluators.gamenet_uq.graph_tools import graph_plotter\n",
    "\n",
    "adsorption_configs = crn.intermediates['VRLIPUYDFBXWCH-UHFFFAOYSA-N*'].ads_configs\n",
    "adsorption_configs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graph_plotter(adsorption_configs['15']['pyg'], node_index=False, dpi=150)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select elementary reactions belonging to a specific reaction type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crn.search_reaction(r_types=['adsorption'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select all the elementary steps involving a specific intermediate, e.g., CO2*. A current limitation is the fact that to select a specific molecule/fragment, you have to first find the corresponding code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crn.search_reaction(inters=['CURLTUGMZLYLDI-UHFFFAOYSA-N*'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select elementary reaction involving a specific intermediate in a specific reaction type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CO2* adsorption\n",
    "\n",
    "crn.search_reaction(r_types=['adsorption'], inters=['CURLTUGMZLYLDI-UHFFFAOYSA-N*'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Search for intermediates having specific stoichiometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crn.search_inter_by_elements({'C': 2, 'H': 1})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CH* C-C coupling\n",
    "\n",
    "crn.search_reaction(r_types=['C-C'], inters=['VRLIPUYDFBXWCH-UHFFFAOYSA-N*'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get NetworkX graph of the reaction network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crn_graph = crn.graph\n",
    "\n",
    "print(crn_graph)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot crn graph with NetworkX function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from networkx import draw\n",
    "\n",
    "draw(crn_graph, with_labels=False, node_size=15, edge_color='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Manipulating elementary reactions\n",
    "\n",
    "CRNs built with CARE contain all the possible elementary reactions (defined by the reaction templates) that could ideally occur within a defined chemical space. From this point, two options are available: (i) directly post-processing the generated CRN by means of microkinetic modeling, or (ii) selecting specific paths and extracting a subset of steps to work with based on your intuition/ideas/etc. Here we show how to manipulate elementary reactions separately. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "r1_idx = random.randint(0, len(crn) - 1)\n",
    "r2_idx = random.randint(0, len(crn) - 1)\n",
    "r3_idx = random.randint(0, len(crn) -1)\n",
    "\n",
    "r1 = crn[r1_idx]\n",
    "r2 = crn[r2_idx]\n",
    "r3 = crn[r3_idx]\n",
    "\n",
    "print(r1, r1.e_rxn)\n",
    "print(r2, r2.e_rxn)\n",
    "print(r3, r3.e_rxn)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiply an elementary reaction by positive scalar value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r4 = r1 * 3.14\n",
    "print(r4, r4.e_rxn) \n",
    "print(type(r4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiply an elementary reaction by negative scalar value a < 0, this equals to reversing the reaction and multiplying it by the absolute scalar value: a * R1 = |a| * (-R1) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r5 = r1 * (-3.14)\n",
    "print(r5, r5.e_rxn)\n",
    "print(type(r5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sum two elementary reactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r5 = r1 + r2\n",
    "print(r5, r5.e_rxn)\n",
    "print(type(r5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subtract one elementary reaction from the other (this equals summing the reverse step): R1 - R2 = R1 + (-R2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r6 = r1 - r2\n",
    "print(r6, r6.e_rxn)\n",
    "print(type(r6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate linear combination of elementary reactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r6 = 0.5 * r1 + 2 * r2 + 0.75 * r3\n",
    "print(r6, r6.e_rxn)\n",
    "print(type(r6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observations:\n",
    "- Applying linear combinations to elementary reactions produces a ``ReactionMechanism`` object, which is a subclass of ``ElementaryReaction``.\n",
    "- The global stoichiometry, as well as the global reaction path energy, are automatically evaluated.\n",
    "- No information about the sequence order of the steps is now included within the class, will be soon included.\n",
    "- One can select the desired steps and construct its own reaction mechanism starting from the pool of elementary reactions in the CRN. Selection of teh desired elementary steps is not trivial and will be improved in the future, for now make use of the ``ReactionNetwork.search_reaction()`` class method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extra: Breaking a molecule in all possible ways\n",
    "\n",
    "If you are interested in quantifying how many bond-breaking reactions could occur given a single starting molecule, this can be easily done with the `dissociate` function, which is the function used within CARE to generate the extended Chemical Space (eCS) in CARE.\n",
    "The required input is a list of smiles strings. Here below we show how to get all the potential bond-breaking reactions and fragments that can be obtained by breaking ethylene glycol (C2H6O2, smiles 'OCCO')."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from care import dissociate\n",
    "\n",
    "smiles_list = ['OCCO']  # you can include more than one smiles!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inters, rxns = dissociate(smiles_list)\n",
    "\n",
    "print('Breaking the following molecules: ', smiles_list)\n",
    "print('Generated intermediates: ', len(inters))\n",
    "print('Bond-breaking reactions: ', len(rxns))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Feel free to try other molecules (Note: not working with input aromatic molecules!)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "carep",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}