{
  "cells": [
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "# Introduction to Python and ASE and some useful libraries\n",
        "\n",
        "Python is a programming language that has become a de-facto standard within scientific programming.  It is also very popular for teaching programming and computer science. One of the advantages of the language is that it is relatively easy to read and understand pre-existing code.\n",
        "\n",
        "### Python Notebooks\n",
        "\n",
        "We are working with Python notebooks.  A notebook is split into cells (this is a text cell, most cells contain a block of Python code).  You can execute a cell by pressing Shift+Enter.  You can edit a text cell by doubleclicking it.  You can change the type of a cell on the menu bar above, use \"Markup\" for text and \"Code\" for Python code.  There are also buttons for adding and removing cells.  Cells can be split into two (or merged) on the Edit menu.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Importing the necesary modules\n",
        "\n",
        "A lot of helper packages for e.g. numeric calculations, plotting, atomic-scale simulations etc are available in Python, but you need to \"import\" the module to get access to it (you also often need to install it, but that has been taken care of here).\n",
        "\n",
        "A few packages support special integration with Notebooks, this includes the plotting package ``matplotlib``, where the integration is enabled with a notebook-specific command starting with the ``%`` character.\n",
        "\n",
        "You can add comments to Python code by using the ``#`` character.\n",
        "\n",
        "Some links for reference:\n",
        "* [Python](https://www.python.org/)\n",
        "* Numerical package [NumPy](http://www.numpy.org/)\n",
        "* Scientific package [Scipy](https://www.scipy.org/)\n",
        "* Plotting package [Matplotlib](https://matplotlib.org/)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "%matplotlib notebook\n",
        "import matplotlib.pyplot as plt    # For nice plotting\n",
        "import numpy as np                 # Mathematical operations"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Let's try some simple stuff.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "2 + 3"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "print('Hello')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "print('Hello ' * 5)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# help is useful for getting help of af method\n",
        "help(print)\n",
        "\n",
        "# you can also use a question mark to get help in a Jupyter notebook,\n",
        "# this opens up a new window (close it if you like)\n",
        "?print"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Variables and data types\n",
        "\n",
        "In python, variables contain data of different types, such as numbers or text strings.  You can print a variable with the ``print()`` function:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "a = 42\n",
        "mypi = 3.14\n",
        "b = \"some text\"\n",
        "print(a)\n",
        "print(\"Pi is approximately\", mypi)\n",
        "print(b)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "The usual mathematical operations are possible.  The operator for exponentiation is the double star.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "print(2 * a)\n",
        "print(16**2)\n",
        "print(a + 2 * mypi)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "### A `list` is an ordered collection of arbitrary objects\n",
        "\n",
        "You will need lists of data. A list of data can be specified at once, or built gradually. The latter is demonstated later.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "primes = [2, 3, 5, 7, 11, 13, 17, 19, 23]\n",
        "print(primes)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Lists are indexed starting with 0, so primes[1] is the *second* prime.  You can also access a list from the end, using negative numbers.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "print(primes[1])\n",
        "print(primes[-1])"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "A `list` can contain arbitrary objects\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# a list\n",
        "l = [1, ('gg', 7), 'hmm', 1.2]\n",
        "print(l)\n",
        "print(l[1])   # Python counts from zero, so this is the second element\n",
        "print(l[-2])  # indexing with negative numbers counts from the end"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "### A `dict`  is a mapping from keys to values\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "d = {'s': 0, 'p': 1}\n",
        "print(d)\n",
        "print(d['p'])\n",
        "print('Removing an element from the dictionary')\n",
        "del d['s']\n",
        "print(d)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "### A `tuple`  is an ordered collection like a list but is *immutable*\n",
        "useful for keywords in `dict`\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# with a list we can reassign values\n",
        "x = [2, 3]\n",
        "x[0] = 100\n",
        "print(x)\n",
        "# this it not possible with a tuple\n",
        "y = (2, 3)\n",
        "print('y = ', y)\n",
        "try:\n",
        "    y[0] = 100\n",
        "except Exception as x:\n",
        "    print(x)\n",
        "print('y = ', y)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "### A ``namedtuple`` is a tuple where the elements also have names.\n",
        "\n",
        "They can be accessed by index or by name.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from collections import namedtuple\n",
        "SKN = namedtuple('IndexSKN', ('spin', 'kpt', 'band'))\n",
        "a = SKN(0, 3, 4)\n",
        "print(a)\n",
        "# The elements may be accessed by index (like a normal tuple) or by name:\n",
        "print(a[1])\n",
        "print(a.kpt)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# lets try and use a namedtuple as keys for a dict\n",
        "d = {}\n",
        "d[SKN(0, 10, 5)] = 3.14\n",
        "d[SKN(0, 1, 3)] = 2.72\n",
        "print(d)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "key = SKN(spin=0, kpt=1, band=3)\n",
        "print(d[key])\n",
        "print(d[(0, 1, 3)])  # one can also use a normal tuple as key"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "# NumPy\n",
        "#### NumPy arrays are heavely used in [ASE](https://wiki.fysik.dtu.dk/ase/)\n",
        "ASE makes heavy use of an extension to Python called NumPy. The NumPy module\n",
        "defines an `ndarray` type that can hold large arrays of uniform\n",
        "multidimensional numeric data. An array is similar to a `list` or a `tuple`,\n",
        "but it is a lot more powerful and efficient.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import numpy as np\n",
        "x = np.array([1, 2, 3])\n",
        "print(x)\n",
        "print(x.mean())"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Multidimensional array\n",
        "a = np.zeros((3, 2))\n",
        "a[:, 1] = 1.0\n",
        "a[1, :] = 2.0\n",
        "print(\"The shape of the array:\", a.shape)\n",
        "print(\"Number of dimensions:\", a.ndim)\n",
        "print(\"Data type:\", a.dtype)\n",
        "print(\"And one can print the entire array (if it is not too big):\")\n",
        "print(a)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Matrix muliplication\n",
        "print('shape of a', a.shape)\n",
        "print('shape of a.T', a.T.shape)  # .T transpose a matrix\n",
        "b = np.dot(a, a.T)\n",
        "print(b)\n",
        "# in a more READABLE way one can use @ to dot matrices together\n",
        "c = a @ a.T\n",
        "print('is c equal to b:', (c == b).all())"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Elementwise multiplication\n",
        "d = a * a\n",
        "print(d)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Random Hermitian matrix\n",
        "seed = 12345678\n",
        "rng = np.random.default_rng(seed)\n",
        "rand = rng.random\n",
        "H = rand((6, 6)) + 1j * rand((6, 6))  # 1j = sqrt(-1)\n",
        "H = H + H.T.conj()\n",
        "\n",
        "# Eigenvalues and eigenvectors\n",
        "eps, U = np.linalg.eig(H)\n",
        "\n",
        "#  Make print of numpy arrays less messy:\n",
        "np.set_printoptions(precision=3, suppress=True)\n",
        "print('The eigenvalues are:', eps.real)\n",
        "\n",
        "# lets try and sort them\n",
        "sorted_indices = eps.real.argsort()\n",
        "eps = eps[sorted_indices]\n",
        "U = U[:, sorted_indices]\n",
        "print('after sorting: ', eps.real)\n",
        "\n",
        "# Check that U diagonalizes H\n",
        "D1 = np.diag(eps)  # Diagonal matrix\n",
        "D2 = U.T.conj() @ H @ U  # Diagonalized H matrix\n",
        "print('Diagonal matrix (from eigenvalues):')\n",
        "print(D1)\n",
        "print('Diagonal matrix (tranforming H):')\n",
        "print(D2)\n",
        "print('Are the numbers in the two matrices close to each other?')\n",
        "print(np.allclose(D2, D1))"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "\n",
        "# Plotting with matplotlib\n",
        " see here for more details [Matplotlib](https://matplotlib.org/)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# In a Jupyter Notebook, this magic line gives nice inline figures,\n",
        "# with interactive possibilities.\n",
        "# This line MUST appear before you import matplotlib or a package\n",
        "# using matplotlib (e.g. ase)\n",
        "%matplotlib notebook\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Some data\n",
        "x = [1.4, 2.0, 3.0, 3.3]\n",
        "y = [10, -1.1, 2.2, 5.0]\n",
        "y2 = [1, 2, 3, 5]"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "plt.figure()   # Start a new figure\n",
        "# Plot y versus x (called \"Series 1\" in the legend).\n",
        "#   b = blue.  o = show points.  - = show line.\n",
        "plt.plot(x, y, 'bo-', label=\"Series 1\")\n",
        "# Plot y2 versus x (called \"Series 2\" in the legend).\n",
        "#   r = red.  x = show points as x'es.   -- = show dashed line\n",
        "plt.plot(x, y2, 'rx--', label=\"Series 2\")\n",
        "plt.title(\"Title of the plot\")\n",
        "plt.xlabel(\"Label of the x-axis\")\n",
        "plt.ylabel(\"Label of the y-axis\")\n",
        "# plt.ylim(0, 10)     # Uncomment to set range of y values to show - similar for x.\n",
        "plt.legend()  # Make a legend, let matplotlib decide where to place it.\n",
        "plt.show()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You can save the plot as a figure by pressing the \"floppy disk\" icon.  Note that in some browsers it does not work, then you can stop the interactive plot by clicking the blue on/off button in the upper right corner, and then save the plot as any other figure in your browser (probably by left-clicking it).\n",
        "\n",
        "Sometimes, you need larger fonts in a plot that you want to include in a report.  Below is the same plot, but with larger fonts.  We also overrule the placement of the legend.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "plt.figure()\n",
        "plt.plot(x, y, 'bo-', label=\"Series 1\")\n",
        "plt.plot(x, y2, 'rx--', label=\"Series 2\")\n",
        "plt.title(\"Title of the plot\", size=24)\n",
        "plt.xlabel(\"Label of the x-axis\", size=16)\n",
        "plt.ylabel(\"Label of the y-axis\", size=16)\n",
        "# Make tick marks larger, and increase the font.\n",
        "plt.tick_params(axis='both', which='major', labelsize=16, size=10)\n",
        "plt.ylim(-0.2, 10.2)     # Set range of y values to show - similar for x.\n",
        "plt.legend(loc=\"upper left\", fontsize=16)  # Make a legend, specify location.\n",
        "plt.tight_layout()   # Fixes that otherwise some of the labels are cropped.\n",
        "plt.show()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "More advanced example with multiple sub-plots.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "fig, axs = plt.subplots(1, 2, sharey=True)\n",
        "x = np.linspace(0, 2 * np.pi, 100)\n",
        "axs[0].plot(x, np.cos(x), label='cos')\n",
        "axs[1].plot(x, np.sin(x), label='sin')\n",
        "axs[0].legend()\n",
        "axs[1].legend()\n",
        "plt.show()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Plotting a countour\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "x = np.linspace(-1, 1, 100)\n",
        "y = np.linspace(-2, 2, 100)\n",
        "X, Y = np.meshgrid(x, y)\n",
        "Z = X**2 + Y**2\n",
        "N = 15\n",
        "\n",
        "fig, ax = plt.subplots(1, 1)\n",
        "ax.contour(X, Y, Z, N)\n",
        "ax.set_aspect('equal')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "# ASE (atomic simulation environment)\n",
        "More details can be found here: https://wiki.fysik.dtu.dk/ase/index.html\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Everything starts with a structure!\n",
        "In ASE the most important ingredients is the `Atom` and `Atoms`\n",
        "objects used to setup an atomic structure.\n",
        "### Setting op a molecule using the `Atoms` object\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from ase import Atoms\n",
        "d = 1.1\n",
        "co = Atoms('CO', positions=[[0, 0, 0], [0, 0, d]])"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# lets try and visualize it using the build in viewer in ase\n",
        "from ase.visualize import view\n",
        "view(co)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "### Setting up a periodic structure\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "d = 2.9\n",
        "L = 10\n",
        "wire = Atoms('Au', positions=[[0, L / 2, L / 2]],\n",
        "             cell=[d, L, L], pbc=[1, 0, 0])\n",
        "# lets try and repeat it and visualize primitive and repeated\n",
        "wire10 = wire * (10, 1, 1)\n",
        "view([wire, wire10])"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Nitrogen on copper\n",
        "### Exercise of the relaxation of a molecule on a surface\n",
        "This section gives a quick (and incomplete) overview of what ASE can do.\n",
        "\n",
        "We will calculate the adsorption energy of a nitrogen molecule on a copper\n",
        "surface. This is done by calculating the total energy for the isolated slab\n",
        "and for the isolated molecule. The adsorbate is then added to the slab and\n",
        "relaxed, and the total energy for this composite system is calculated. The\n",
        "adsorption energy is obtained as the sum of the isolated energies minus the\n",
        "energy of the composite system.\n",
        "\n",
        "You can read more about the optimizers in ASE here:\n",
        "https://wiki.fysik.dtu.dk/ase/ase/optimize.html\n",
        "\n",
        "#### 1. Try to go through the script so you understand what is going on\n",
        "#### 2. Calculate the adsorption energy of N2 on a 4x4x2 fcc111 slab (result= 0.324 eV)\n",
        "#### 3. Try a couple of different optimizers and see which one is the fastest\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from ase import Atoms\n",
        "from ase.calculators.emt import EMT\n",
        "from ase.constraints import FixAtoms\n",
        "from ase.optimize import QuasiNewton, BFGS, FIRE\n",
        "from ase.build import fcc111, add_adsorbate\n",
        "from ase.visualize import view\n",
        "\n",
        "h = 2.85\n",
        "d = 1.10\n",
        "\n",
        "# Set up a slab:\n",
        "slab = fcc111('Cu', size=(4, 4, 2), vacuum=10.0)\n",
        "slab.calc = EMT()\n",
        "e_slab = slab.get_potential_energy()\n",
        "# setup a molecule\n",
        "molecule = Atoms('2N', positions=[(0., 0., 0.), (0., 0., d)])\n",
        "molecule.calc = EMT()\n",
        "e_N2 = molecule.get_potential_energy()\n",
        "\n",
        "add_adsorbate(slab, molecule, h, 'ontop')\n",
        "constraint = FixAtoms(mask=[a.symbol != 'N' for a in slab])\n",
        "slab.set_constraint(constraint)\n",
        "dyn = QuasiNewton(slab, trajectory='N2Cu.traj')\n",
        "dyn.run(fmax=0.025)\n",
        "\n",
        "print('Adsorption energy:', e_slab + e_N2 - slab.get_potential_energy())"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Visualize the trajectory\n",
        "!ase gui N2Cu.traj"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Band structure\n",
        "#### Using ASE to setup band structures for Al using a Freelectron model and DFT\n",
        "#### 1. What is the crystal structure of Al?\n",
        "#### 2. Try and look up the recommeded Brillouin zone path for crystal structure [here](https://wiki.fysik.dtu.dk/ase/ase/dft/kpoints.html)\n",
        "#### 3. Can you figure out what the `nbands=-10` and `convergence={'bands': -5}` parameters means in the GPAW DFT input below ? (Hint try and look at the output file `Al.txt`)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# simple free electron calculator\n",
        "from ase.build import bulk\n",
        "from ase.calculators.test import FreeElectrons\n",
        "\n",
        "a = bulk('Al')\n",
        "kpts = {'path': 'GXWLGK', 'npoints': 100}\n",
        "\n",
        "# Simple FreeElectron model calculator\n",
        "a.calc = FreeElectrons(nvalence=3,\n",
        "                       kpts=kpts)\n",
        "a.get_potential_energy()\n",
        "bs = a.calc.band_structure()\n",
        "bs.plot(emax=10, filename='al-free-electron.png')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# setup a DFT calculation with GPAW and repeat\n",
        "from gpaw import GPAW, PW\n",
        "# calc the self-consistent electron density\n",
        "a.calc = GPAW(kpts=(3, 3, 3), mode=PW(200), txt='Al.txt')\n",
        "a.get_potential_energy()\n",
        "# band-structure calculation for a fixed density\n",
        "calc = a.calc.fixed_density(\n",
        "    kpts=kpts,\n",
        "    symmetry='off',\n",
        "    nbands=-10,\n",
        "    convergence={'bands': -5})"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "bs = calc.band_structure()\n",
        "bs.plot(emax=10, filename='al-dft.png')\n"
      ],
      "outputs": [],
      "execution_count": null
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "CAMD2022",
      "language": "python",
      "name": "camd2022"
    },
    "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.6.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 1
}