SOAP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Smooth Overlap of Atomic Positions\n",
    "\n",
    "SOAP is a **local** descriptor, that maps the local environment around a point very accurately. It eliminates rotational, and permutation redundancies by integrating the overlap of smoothed out atomic positions, by gaussian smearing, and mapping them into coefficients of orthornormal basis functions.\n",
    "\n",
    "This is done by the following steps:\n",
    "\n",
    "* Smooth out the atomic positions:\n",
    "\n",
    "    The atomic positions are point objects in space. Integrating them would need a lot of basis functions. Thus, the atoms' positions are smeared as gaussian functions. \n",
    "    $$ \\rho(r) = \\sum_i e^{-(r-r_i)^2}$$\n",
    "    However, this also makes all the elements indistinguishable. Thus, SOAP for individual elements, in molecule/unit-cell, is calculated, and then the values are concantenated at end. \n",
    "    ![soap depiction](./images/soap_depiction.png) *Image courtesy Jäger Marc*\n",
    "\n",
    "\n",
    "* Generate orthonormal basis set:\n",
    "\n",
    "    The obtained smeared atomic position, or atomic density, if you will, is decomposed using Laplace Spherical Harmonics -- spherical harmonics in real space -- and orthogonal basis set: $\\Upsilon_{lm}(\\theta, \\phi)$ and  $g_n(r) $. \n",
    "    Basis function for s orbital of hydrogen: ![basis function](./images/basis_set.jpg) \n",
    "    Laplace spherical harmonics $\\Upsilon_{ℓm}$ for l = 0, …, 4 (top to bottom) and m = 0, …, l (left to right). The negative order harmonics $\\Upsilon_{ℓ-m}$ would be shown rotated about the z axis by $90^o$ with respect to the positive order ones.\n",
    "    ![rotating spherical harminics](./images/Rotating_spherical_harmonics.gif)  *Image courtsey wikipedia.org/wiki/User:Cyp*  \n",
    "    \n",
    "\n",
    "* Integrate for all coefficients:\n",
    "    $$c_{nlm} = \\left< \\rho | g_n(r)\\Upsilon_{lm} \\right>  = \\int_V g_n(r)\\Upsilon_{lm}(\\theta, \\phi)\\rho(r, \\theta, \\phi)dV$$\n",
    "    Further, a power spectrum, or a density matrix, per se, is made out of these parameters and summed for all m's for rotational invarience.\n",
    "    $$P_{nn'l} = \\sum_m c_{nlm}c^*_{n'lm}$$\n",
    "    \n",
    "\n",
    "For more info see:\n",
    "[Bartók, Albert P., Risi Kondor, and Gábor Csányi. <i>Physical Review B 87.18</i> (2013): <b>184115</b>](https://arxiv.org/pdf/1209.3140.pdf;)\n",
    "\n",
    "For calculating SOAP, we use the [DScribe package](https://singroup.github.io/dscribe) as developed by [Surfaces and Interfaces at the Nanoscale, Aalto](http://physics.aalto.fi/en/groups/sin/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "\n",
    "We are going to see SOAP in action for a simple NaCl system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --- INITIAL DEFINITIONS ---\n",
    "import numpy, math, random\n",
    "from visualise import view\n",
    "from ase import Atoms\n",
    "import sys\n",
    "sys.path.insert(0, './data/descriptor_codes/')\n",
    "sys.path.insert(0, './data/descriptor_codes/src')\n",
    "from dscribe.descriptors import SOAP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Atom description\n",
    "\n",
    "We'll make an ase.Atoms class for NaCl:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the system under study: NaCl in a conventional cell.\n",
    "NaCl_conv = Atoms(\n",
    "    cell=[\n",
    "        [5.6402, 0.0, 0.0],\n",
    "        [0.0, 5.6402, 0.0],\n",
    "        [0.0, 0.0, 5.6402]\n",
    "    ],\n",
    "    scaled_positions=[\n",
    "        [0.0, 0.5, 0.0],\n",
    "        [0.0, 0.5, 0.5],\n",
    "        [0.0, 0.0, 0.5],\n",
    "        [0.0, 0.0, 0.0],\n",
    "        [0.5, 0.5, 0.5],\n",
    "        [0.5, 0.5, 0.0],\n",
    "        [0.5, 0.0, 0.0],\n",
    "        [0.5, 0.0, 0.5]\n",
    "    ],\n",
    "    symbols=[\"Na\", \"Cl\", \"Na\", \"Cl\", \"Na\", \"Cl\", \"Na\", \"Cl\"],\n",
    ")\n",
    "view(NaCl_conv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting SOAP hyper-parameters\n",
    "Next we set the hyper-parameters to SOAP. \n",
    "1. calcpos, center of SOAP calculation\n",
    "1. rcut, sets the cutoff for atoms whose gaussian densities will be included in the integral.\n",
    "2. nmax, sets the number of orthogonal radial basis functions to use.\n",
    "3. lmax, sets the number of angular momentum terms, so l = 0, 1, ...., lmax\n",
    "\n",
    "    **Note: even when giving one SOAP calculation position, it should be wrapped in a list, as shown in example below**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computing SOAP\n",
    "calcpos = [0, 0, 0]\n",
    "\n",
    "soaper = SOAP(\n",
    "    rcut=8,\n",
    "    nmax=5,\n",
    "    lmax=5,\n",
    "    species=['Na', 'Cl'],\n",
    "    sparse=False\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculation\n",
    "Now we call the soap function, and pass all the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculation\n",
    "soap1 = soaper.create(NaCl_conv, positions=[calcpos])\n",
    "print(\"Size of descriptor: {}\".format(soap1.shape[1]))\n",
    "print(\"First five values, for position {}: {}\".format(calcpos, soap1[0,:5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rotational invariance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Rotation of positions\n",
    "\n",
    "print(\"Original positions:\\n {}\".format(NaCl_conv.positions))\n",
    "\n",
    "NaCl_conv.rotate(90, [0,1,1], center=calcpos)\n",
    "\n",
    "print(\"Rotated positions:\\n {}\".format(NaCl_conv.positions))\n",
    "\n",
    "view(NaCl_conv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recompute SOAP for the same atom, after rotation and show the difference in descriptors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soap2 = soaper.create(NaCl_conv, positions=[calcpos])\n",
    "print(numpy.linalg.norm(soap1 - soap2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remark\n",
    "\n",
    "The power spectrum at a desired position x is the fingerprint of the local chemical environment at this specific position. Thus, it can be used to:\n",
    "1. Compare the similarity of two local chemical environments by comparing their SOAP descriptors.\n",
    "2. Machine learn local properties, like charges, adsorption energies, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Smoothness\n",
    "\n",
    "Verify that the SOAP is smooth under translations of point of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DIY..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Construct a global SOAP\n",
    "Use the atomic environments to construct an average SOAP descriptor for molecules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# atomic positions as matrix\n",
    "molxyz = numpy.load(\"./data/molecule.coords.npy\")\n",
    "# atom types\n",
    "moltyp = numpy.load(\"./data/molecule.types.npy\")\n",
    "\n",
    "atoms_sys = Atoms(positions=molxyz, numbers=moltyp)\n",
    "view(atoms_sys)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# build SOAP at each atom location\n",
    "# ...\n",
    "# compute average soap for each specie\n",
    "# ...\n",
    "# concatenate the soaps to the the overall global one\n",
    "# ..."
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}