engn2770_computational_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ENGN 2770: Atomistic Reaction Engineering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computational Tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the pdf version. The interactive version (which includes necessary environment files) can be accessed from using the link provided: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://mybinder.org/v2/gh/bjkreitz/ase-tutorial-engn2770/main"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Note: This can take a minute or two to load and download required environment files)\n",
    "\n",
    "Once you open the link, you can download the \"engn2770_computational_tutorial.ipynb\" and save it on your local computer. \n",
    "\n",
    "Note, you will need jupyter notebook along with other necessary packages if you wish to use it on your local machine(ASE, python, matplotlib, numpy). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A) Intro to linux bash commands (https://files.fosswire.com/2007/08/fwunixref.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Image(filename='extras/linux_commands.PNG')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image(filename='extras/ccv_commands.PNG')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Setting up a Linux Virtual Machine**\n",
    "\n",
    "For the project, it will be easier to work on a Linux system. You don't have to install Linux on a separate partition (although you can if you want), but use a Virtual Machine instead. A guide to set up a Linux Virtual Machine can be found here: http://reactionmechanismgenerator.github.io/RMG-Py/users/rmg/installation/virtualMachineSetup.html \n",
    "\n",
    "I have good experience with VirtualBox and the setup for this is described in 3.1.1.5. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## B) Python "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Basic arithmetic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Basic arithmetic\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# square of a number\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# division & % \n",
    "\n",
    "# remainder of a division\n",
    "\n",
    "# If this is 0, then 8 is divisible by 2!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# functions\n",
    "\n",
    "# define add function\n",
    "\n",
    "# define power function\n",
    "\n",
    "# call add function\n",
    "\n",
    "# call power function\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Python has a lot of simple builtin functions\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### numpy, list & arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# using numpy (Numerical Python)\n",
    "import numpy as np\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# list \n",
    "\n",
    "for thing in things:\n",
    "    print(thing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Basic list commands\n",
    "\n",
    "_list = [1,2,3,4]\n",
    "\n",
    "\n",
    "print(_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# return evenly spaced numbers over a specified interval\n",
    "\n",
    "\n",
    "\n",
    "# convert a 1D array into 2D \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### For loop & conditional statements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for loops\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For loops with conditional statements\n",
    "\n",
    "dice = [1,2,3,4,5,6]\n",
    "\n",
    "type_list = []\n",
    "\n",
    "print(type_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num = np.linspace(0,10,11)\n",
    "print(num)\n",
    "\n",
    "doubled_list = []\n",
    "\n",
    "for element in num:\n",
    "    if element > 5:\n",
    "        doubled_list.append(element*2)\n",
    "\n",
    "print(doubled_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For loop with conditional statement in one line\n",
    "\n",
    "num = np.linspace(0,10,11)\n",
    "print(num)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "long_words = ['platinum', 'hi', 'carbon', 'welcome', 'hydrogen', 'neptune', 'covid', 'python']\n",
    "\n",
    "# checks and stores into new list based on length of a word\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## C) Atomic Simutlation Environment "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some basic packages we will need throughout\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "\n",
    "import ase\n",
    "from ase.parallel import paropen\n",
    "from ase.io import read, write\n",
    "from ase.visualize import view\n",
    "\n",
    "# Plotting and image visualizing in Jupyter Notebook\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import Image\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module 1: Atom / Atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the module (use tab tab for auto-completion)\n",
    "from ase import Atom, Atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use ? to know more\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "#Print properties of this atom\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create an N2 molecule\n",
    "d = 1.104 # N2 bondlength\n",
    "\n",
    "# The following three are equivalent\n",
    "\n",
    "#N2 = Atoms(numbers=[7, 7], positions=[(0, 0, 0), (0, 0, d)])\n",
    "\n",
    "# to save and view a snapsho\n",
    "\n",
    "\n",
    "# to have a 3D view in notebook\n",
    "#view(N2, viewer='x3d')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = 2.9\n",
    "L = 10.0\n",
    "\n",
    "# construct a Au wire\n",
    "wire = Atoms('Au')\n",
    "\n",
    "wire.get_chemical_symbols()\n",
    "view(wire)\n",
    "\n",
    "# to save and view a snapshot\n",
    "write('extras/wire.png', wire)\n",
    "Image(filename='extras/wire.png') \n",
    "\n",
    "\n",
    "#view(wire*(6,1,1), viewer='x3d')\n",
    "\n",
    "#view(wire*(6,1,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module 2: Build (molecule, bulk, surfaces)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.build import molecule\n",
    "\n",
    "\n",
    "# to save and view a snapshot\n",
    "write('extras/ethanol.png', m, rotation='-20x')\n",
    "Image(filename='extras/ethanol.png',width=150)\n",
    "\n",
    "#view(m, viewer='x3d')\n",
    "\n",
    "#view(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimize lattice constant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.build import bulk\n",
    "from ase.calculators.emt import EMT\n",
    "from ase.io import Trajectory\n",
    "\n",
    "a0 = 3.6\n",
    "\n",
    "# to save and view a snapshot\n",
    "write('extras/cu.png', cu)\n",
    "Image(filename='extras/cu.png',width=150)\n",
    "\n",
    "#view(cu, viewer='x3d')\n",
    "\n",
    "#view(cu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create an empty list to store energy\n",
    "energy_list = []\n",
    "\n",
    "# start with a good guess\n",
    "# and create a list with fluctuations in a.\n",
    "\n",
    "a0 = 3.6\n",
    "a_list = a0 + np.linspace(-0.2, 0.2, 101)\n",
    "\n",
    "# For loop to use each and every value of a\n",
    "# and calculate & store the potential energy\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting using matplotlib (https://matplotlib.org/)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(7,7))\n",
    "plt.scatter(a_list, energy_list, marker='o', s=10)\n",
    "plt.xlabel('$a$, $\\mathrm{\\AA}$', fontsize=15)\n",
    "plt.xticks(fontsize=15)\n",
    "plt.ylabel('potential enegry, eV', fontsize=15)\n",
    "plt.yticks(fontsize=15)\n",
    "fig.savefig('extras/lattice_constant.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# index corresponding to the lowest energy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module 3:  Surface, Constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.build import fcc111\n",
    "from ase.constraints import FixAtoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fcc111? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a slab using build\n",
    "\n",
    "# Use of constraints, here, Fix the bottom most layer\n",
    "\n",
    "# Note that the constraint is not visible through \n",
    "# the following method\n",
    "\n",
    "write('extras/cu_slab.png', slab, rotation='10z,-80x')\n",
    "Image(filename='extras/cu_slab.png') \n",
    "\n",
    "#view(slab, viewer='x3d')\n",
    "#view(slab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module 4: Add an adsorbate "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### method 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.build import add_adsorbate\n",
    "\n",
    "# create a slab\n",
    "slab = fcc111('Cu', size=(3,3,3), a=3.588, vacuum=10)\n",
    "slab.set_pbc((1,1,0))\n",
    "\n",
    "# create an atom and use add_adsorbate to add it to the slab\n",
    "\n",
    "write('extras/slab_with_adsorbate.png', slab * (2, 2, 1))\n",
    "Image(filename='extras/slab_with_adsorbate.png') \n",
    "\n",
    "#view(slab, viewer='x3d')\n",
    "#view(slab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### method 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 3.558\n",
    "\n",
    "# create a slab\n",
    "slab = fcc111('Cu', size=(3,3,3), a=3.588, vacuum=10)\n",
    "slab.set_pbc((1,1,0))\n",
    "\n",
    "# use the slab atom position to add it to the slab\n",
    "\n",
    "write('extras/slab_with_adsorbate.png', slab_with_add)\n",
    "Image(filename='extras/slab_with_adsorbate.png') \n",
    "\n",
    "#view(slab_with_add, viewer='x3d')\n",
    "#view(slab_with_add)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module 5: Structure optimization using EMT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.constraints import FixAtoms\n",
    "from ase.build import fcc111, add_adsorbate\n",
    "from ase.calculators.emt import EMT\n",
    "from ase.optimize import BFGS\n",
    "\n",
    "# Construct a slab with optimized lattice constant\n",
    "slab = fcc111('Cu', size=(3,3,3), a=3.588, vacuum=10)\n",
    "# Periodic boundary conditions\n",
    "slab.set_pbc((1,1,0))\n",
    "\n",
    "# Construct & add an adosorbate\n",
    "adsorbate = Atom('O')\n",
    "add_adsorbate(slab, adsorbate, 1.8, 'ontop')\n",
    "\n",
    "# Fix atoms\n",
    "indices=[atom.index for atom in slab if atom.tag == 3]\n",
    "constraint = FixAtoms(indices)\n",
    "slab.set_constraint(constraint)\n",
    "\n",
    "# Calculator\n",
    "\n",
    "\n",
    "# Optimization\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Open the relaxed slab\n",
    "#relaxed_slab = read('qn.traj')\n",
    "#view(relaxed_slab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Read the log file\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Put the same adsorbate in a different position\n",
    "\n",
    "from ase.constraints import FixAtoms\n",
    "from ase.build import fcc111, add_adsorbate\n",
    "from ase.calculators.emt import EMT\n",
    "from ase.optimize import BFGS\n",
    "\n",
    "# Construct a slab with optimized lattice constant\n",
    "slab = fcc111('Cu', size=(3,3,3), a=3.588, vacuum=10)\n",
    "# periodic boundary condition \n",
    "slab.set_pbc((1,1,0))\n",
    "\n",
    "# Construct & add an adosorbate\n",
    "adsorbate = Atom('O')\n",
    "add_adsorbate(slab, adsorbate, 1.8, 'hcp')\n",
    "\n",
    "# Fix atoms constraint to fix the last atomic layer\n",
    "indices=[atom.index for atom in slab if atom.tag == 3]\n",
    "constraint = FixAtoms(indices)\n",
    "slab.set_constraint(constraint)\n",
    "\n",
    "# Calculator\n",
    "calc = EMT()\n",
    "slab.set_calculator(calc)\n",
    "\n",
    "if os.path.exists('qn.log'):\n",
    "    os.remove('qn.log')\n",
    "\n",
    "# Optimization\n",
    "opt = BFGS(slab, logfile='qn.log', trajectory='qn.traj')\n",
    "opt.run(fmax=0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myfile = open(\"qn.log\",'r')\n",
    "txt = myfile.read()\n",
    "print(txt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module 6: Calculate Barriers using nudged elastic band (NEB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "import ase \n",
    "from ase.io import read, write\n",
    "from ase import Atom\n",
    "from ase.build import fcc111\n",
    "from ase.constraints import FixAtoms\n",
    "from ase.neb import NEB\n",
    "from ase.optimize import BFGS\n",
    "from ase.calculators.emt import EMT\n",
    "\n",
    "# Create a bare slab\n",
    "slab = fcc111('Cu', size=(3,3,3), a=3.588, vacuum=10)\n",
    "slab.set_pbc((1,1,0))\n",
    "\n",
    "indices=[atom.index for atom in slab if atom.tag == 3]\n",
    "constraint = FixAtoms(indices)\n",
    "slab.set_constraint(constraint)\n",
    "\n",
    "slab.set_calculator(EMT())\n",
    "\n",
    "opt = BFGS(slab, logfile='bare_slab.log', trajectory='bare_slab.traj')\n",
    "opt.run(fmax=0.01)\n",
    "\n",
    "# create an adsorbate\n",
    "adsorbate = Atom('O')\n",
    "\n",
    "#######################################\n",
    "###### SLAB 1: Cu with O in fcc (position 1) #######\n",
    "#######################################\n",
    "\n",
    "initial = slab + adsorbate\n",
    "initial[-1].position = initial[18].position + (1.3,0.8,1)\n",
    "\n",
    "initial.set_calculator(EMT())\n",
    "opt = BFGS(initial, logfile='initial.log', trajectory='initial.traj')\n",
    "opt.run(fmax=0.01)\n",
    "\n",
    "####################################################\n",
    "###### SLAB 1: Cu with O in fcc (position 2) #######\n",
    "####################################################\n",
    "final = slab + adsorbate\n",
    "final[-1].position = final[19].position + (1.3, 0.8, 1)\n",
    "\n",
    "final.set_calculator(EMT())\n",
    "opt = BFGS(final, logfile='final.log', trajectory='final.traj')\n",
    "opt.run(fmax=0.01)\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}