03_Introduction_to_ASE_Calculators.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ASE Calculators\n",
    "\n",
    "In this lesson, we will cover the basics of ASE calculators and how to use them. ASE has the ability to attach a \"calculator\" to an atoms object. At a fundamental level, these are just classes that will return energies and forces. That means it can be a DFT program, a quantum chemistry program, or even just a classical potential like \"Leonard-Jones\" There are [a ton](https://wiki.fysik.dtu.dk/ase/ase/calculators/calculators.html) of calculators implmented. Here we will just use a simple one, the [EMT](https://wiki.fysik.dtu.dk/ase/ase/calculators/emt.html) calculator.\n",
    "A calculator acts as a black box that takes in atomic numbers and positions to give out energy, forces and stresses. \n",
    "\n",
    "The `EMT` calculator is just a simple pair potential calculator for a few metals. In practice it is a toy calculator used for testing.\n",
    "\n",
    "You start by making an instance of your calculator, then using `set_calculator` to attach it to an atoms object. Once that is done you can call `get_potential_energy` and `get_forces` to calculate the energy and forces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.calculators.emt import EMT\n",
    "\n",
    "calc = EMT()\n",
    "water.set_calculator(calc)\n",
    "energy = water.get_potential_energy()\n",
    "forces = water.get_forces()\n",
    "\n",
    "print(energy)\n",
    "print(forces)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our group we use quatum espresso, VASP, and Abinit to do DFT in general.\n",
    "\n",
    "A primary use of calculators is to perform structural optimizations. This allows us to find the \"lowest energy configuration\" of a given structure. ASE has [tools](https://wiki.fysik.dtu.dk/ase/ase/optimize.html) to do this. Below we're using the `BFGSLineSearch` optimization method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.optimize import BFGSLineSearch\n",
    "\n",
    "relax = BFGSLineSearch(atoms = water)\n",
    "relax.run(fmax = 0.05) # relax the structure until the maximum force is 0.05 eV/A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view(water)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise:\n",
    "\n",
    "Use the EMT calculator to plot the potential energy as a function of distance between two H atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Solution to the exercise:\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "rs = np.linspace(0.4, 3, 50)\n",
    "engs = []\n",
    "## calculate energy as a function of r\n",
    "for r in rs:\n",
    "    atoms = molecule('H2')\n",
    "    atoms.set_distance(0,1,r)\n",
    "    atoms.set_calculator(EMT())\n",
    "    engs.append(atoms.get_potential_energy())\n",
    "    \n",
    "plt.plot(rs, engs)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating Reaction Energies\n",
    "\n",
    "Reaction energies are the core of computational catalysis and surface science because they provide fundamental information about the energy for a chemical reaction. This section is meant to cover the basics of how to calculate a reaction energy using DFT.\n",
    "\n",
    "The adsorption energy is defined as the energy difference between the combined system and the separate systems:\n",
    "\n",
    "$E_{rxn} = \\Sigma E_{products} - \\Sigma E_{reactants}$\n",
    "\n",
    "here we're going to calculate the reaction energy of forming water for O$_2$ and H$_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We'll start by making an O$_2$ molecule\n",
    "from ase.build import molecule\n",
    "\n",
    "O2 = molecule('O2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#let's visualize it next\n",
    "from ase.visualize import view\n",
    "\n",
    "view(O2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#To calculate the energy of this molecule, we're going to use the EMT calculator\n",
    "from ase.calculators.emt import EMT\n",
    "\n",
    "O2.set_calculator(EMT())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Next we need to optimize the structure, because chances are it's not in its lowest energy configuration\n",
    "from ase.optimize import QuasiNewton\n",
    "\n",
    "dyn = QuasiNewton(O2)\n",
    "dyn.run(fmax=0.05)\n",
    "\n",
    "\n",
    "E_O2 = O2.get_potential_energy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's build and optimize H$_2$ and H$_2$O. ASE has some nice tools for doing this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H2 = molecule('H2')\n",
    "view(H2, viewer='nglview')\n",
    "\n",
    "H2.set_calculator(EMT())\n",
    "dyn = QuasiNewton(H2)\n",
    "dyn.run(fmax=0.05)\n",
    "\n",
    "E_H2 = H2.get_potential_energy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H2O = molecule('H2O')\n",
    "view(H2O, viewer='nglview')\n",
    "\n",
    "H2O.set_calculator(EMT())\n",
    "dyn = QuasiNewton(H2O)\n",
    "dyn.run(fmax=0.05)\n",
    "\n",
    "E_H2O = H2O.get_potential_energy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E_rxn = E_H2O - E_H2 - 0.5*E_O2\n",
    "print(E_rxn) # answer is in eV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this is not very accurate! The real value is approximately -285.8261 kJ/mol. The reason is that we used the \"EMT\" calculator, which is not appropriate for the physics of molecules. However, it is very fast, so it allows us to see the principles behind how reaction energies can be computed. If we want accurate numbers, we need to use quantum mechanics, approximated through Density Functional Theory (DFT). This will be much more computationally expensive, and will require a supercomputer. You will learn more about this in future lessons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.units import kJ, mol\n",
    "true_value = -285.8261 * kJ / mol\n",
    "print(true_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error = (E_rxn-true_value)/true_value\n",
    "print(error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# More about ASE Calculators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Supported calculators in ASE\n",
    "\n",
    "ASE supports a variety of calculators. List of calculators either available through native ASE interfaces or python implementations in ASE package can be found [here](https://wiki.fysik.dtu.dk/ase/ase/calculators/calculators.html). The Calculator takes in keywords which can be parsed in by user based on the type of calculations they want to perform. Some example of keywords for a DFT calculator are: **xc** (exchange-correlation functional), **kpts** (k-points), **smearing**, **nbands**. Each calculator will have different native keywords specific to that calculator. These arguments can also be set or changed later by using the **set()** method. \n",
    "\n",
    "## Calculator interface in ASE\n",
    "\n",
    "An ASE calculator will store a copy of atoms object that is used for calculation. To print out the properties after a calculation has been performed, one can use *get_potential_energy()*, *get_forces()*, *get_stress()* methods. "
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}