scripts/AEC_procedure.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f2fbb9dd",
   "metadata": {},
   "source": [
    "This notebook has 2 main steps\n",
    "- First, it generates an `input.yml` file for ARC to run single point calculations on the experimental geometries for the 16 reference species in our database.\n",
    "- After the ARC calculations have completed, the notebook uses Arkane to perform AEC fitting.\n",
    "\n",
    "Briefly, [ARC](https://github.com/ReactionMechanismGenerator/ARC) (Automated Rate Calculator) is a software developed by the RMG team for automating electronic structure calculations, which are then used to calculate thermochemical and kinetic parameters. Installation instructions can be found on ARC's documentation [link](https://reactionmechanismgenerator.github.io/ARC/installation.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bc597a2",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import yaml\n",
    "\n",
    "from arc.species.converter import xyz_to_str\n",
    "\n",
    "from arkane.encorr.reference import ReferenceDatabase\n",
    "from arkane.encorr.ae import SPECIES_LABELS, AEJob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82bc18fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "############# USER INPUT #############\n",
    "LOT = {'method': 'CCSD(T)-F12',\n",
    "       'basis':'cc-pVTZ-F12',\n",
    "       'software':'molpro',\n",
    "       }\n",
    "arc_project_name = 'aec_bac'\n",
    "#######################################\n",
    "model_chemistry = '/'.join([LOT.get(key, '') for key in ['method', 'basis']])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93711d41",
   "metadata": {},
   "source": [
    "# Step 1: Write ARC input file to conveniently run single point calculations in preparation for AEC fitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e8b1e58",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# load database\n",
    "database = ReferenceDatabase()\n",
    "database.load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c6bf5c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create dictionary to map species labels to folder names from ARC's output\n",
    "SPECIES_LABELS = {\n",
    "    'Br2': 'Dibromine',\n",
    "    'BrH': 'Hydrogen bromide',\n",
    "    'CH3': 'Methyl',\n",
    "    'CH3Cl': 'Chloromethane',\n",
    "    'CH4': 'Methane',\n",
    "    'Cl2': 'Dichlorine',\n",
    "    'ClH': 'Hydrogen chloride',\n",
    "    'F2': 'Difluorine',\n",
    "    'FH': 'Hydrogen fluoride',\n",
    "    'H2': 'Dihydrogen',\n",
    "    'H2O': 'Water',\n",
    "    'H2S': 'Hydrogen sulfide',\n",
    "    'H3N': 'Ammonia',\n",
    "    'N2': 'Dinitrogen',\n",
    "    'O2': 'Dioxygen',\n",
    "    'S2': 'Disulfur',\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1f6836c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load experimental geometries from reference database\n",
    "exp_geometries = {}\n",
    "for label in SPECIES_LABELS.values():\n",
    "    ref_spec = database.get_species_from_label(label)[0]\n",
    "    xyz = ref_spec.reference_data['CCCBDB'].xyz_dict\n",
    "    exp_geometries[label] = xyz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36a7dc6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define input file structure\n",
    "input_file = f\"\"\"project: {arc_project_name}\n",
    "sp_level: {model_chemistry}\n",
    "\n",
    "compute_thermo: false\n",
    "\n",
    "job_types:\n",
    "  conformers: false\n",
    "  opt: false\n",
    "  fine_grid: false\n",
    "  freq: false\n",
    "  sp: true\n",
    "  rotors: false\n",
    "\n",
    "species:\"\"\"\n",
    "\n",
    "species = \"\"\"\n",
    "  - label: {formula}\n",
    "    smiles: '{smiles}'\n",
    "    charge: {charge}\n",
    "    multiplicity: {multiplicity}\n",
    "    xyz: |\n",
    "      {xyz}\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3fdcec5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create input yaml file for ARC\n",
    "for spc in exp_geometries.keys():\n",
    "    label = spc    \n",
    "    xyz = xyz_to_str(exp_geometries[label])\n",
    "\n",
    "    xyz_indented = ''\n",
    "    for i, line in enumerate(xyz.split('\\n')):\n",
    "        if i == 0:\n",
    "            xyz_indented += f'{line}\\n'\n",
    "        else:\n",
    "            xyz_indented += f'      {line}\\n'\n",
    "    \n",
    "    ref_spec = database.get_species_from_label(label)[0]\n",
    "    charge = ref_spec.charge\n",
    "    multiplicity = ref_spec.multiplicity\n",
    "    formula = ref_spec.formula\n",
    "    smiles = ref_spec.smiles\n",
    "    \n",
    "    input_file += species.format(formula=formula,\n",
    "                                 smiles=smiles,\n",
    "                                 charge=charge,\n",
    "                                 multiplicity=multiplicity,\n",
    "                                 xyz=xyz_indented)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "637dd103",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"input.yml\", \"w\") as f:\n",
    "    f.write(input_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccba638e",
   "metadata": {},
   "source": [
    "Now use ARC to run this input file and obtain single point energies for the experimental geometries"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ece58cce",
   "metadata": {},
   "source": [
    "# Step 2: After running sp calculations, get AEC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "097921c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from arkane.ess.factory import ess_factory\n",
    "from arkane.modelchem import LevelOfTheory\n",
    "import rmgpy.constants as constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67255b69",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define path to the \"Species\" folder from ARC's output\n",
    "sp_dir = 'calcs/Species/'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45eb4483",
   "metadata": {},
   "source": [
    "`arkane/encorr.ae.py` says:\n",
    "\n",
    "```\n",
    "Notes:\n",
    "    The species energies should be provided as a dictionary\n",
    "    containing the species labels as keys and their single-\n",
    "    point electronic energies in Hartree as values. The\n",
    "    energies should be calculated using the experimental\n",
    "    geometry provided for the species in the reference\n",
    "    database, and **the zero-point energy should not be included\n",
    "    in the electronic energy.**\n",
    "            \n",
    "```\n",
    "\n",
    "The `load_energy()` method on Arkane adapters adapter returns the energy in J/mol and does NOT include the zero-point energy. So just convert the energy from J/mol to Hartrees before fitting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a7abb11",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the energy for each QM sp job\n",
    "species_energies = {}\n",
    "for root, dirs, files in os.walk(sp_dir):\n",
    "    dirs.sort()\n",
    "    print(len(dirs))\n",
    "    print(dirs)\n",
    "    for directory in dirs:\n",
    "        for root1, dirs1, files1 in os.walk(os.path.join(root, directory)):\n",
    "            for d1 in dirs1:\n",
    "                for root2, dirs2, files2 in os.walk(os.path.join(root1, d1)):\n",
    "                    files2.sort()\n",
    "                    # extract the output.out file\n",
    "                    output = os.path.join(root2, 'output.out')\n",
    "                    log = ess_factory(output)\n",
    "                    energy = log.load_energy()\n",
    "                    energy = energy / (constants.E_h * constants.Na)  # convert from J/mol to Hartrees\n",
    "                    species_energies.update({SPECIES_LABELS[directory]: energy})\n",
    "    break           \n",
    "species_energies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4c13cf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "lot = LevelOfTheory(**LOT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6367075",
   "metadata": {},
   "outputs": [],
   "source": [
    "ae = AEJob(species_energies=species_energies,\n",
    "           level_of_theory=lot,\n",
    "           )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e48d0389",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_file = 'AEC_' + '_'.join([LOT.get(key, '') for key in ['method', 'basis']])\n",
    "ae.execute(output_file=f'{output_file}.out')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1c34b28",
   "metadata": {},
   "source": [
    "Copy paste the values from the output file to `RMG-database/input/quantum_corrections/data.py`"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}