Workshop_2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<script>\n",
       "code_show=true; \n",
       "function code_toggle() {\n",
       "    if (code_show){\n",
       "        $('div.cell.code_cell.rendered.selected div.input').hide();\n",
       "    } else {\n",
       "        $('div.cell.code_cell.rendered.selected div.input').show();\n",
       "    }\n",
       "    code_show = !code_show\n",
       "} \n",
       "$( document ).ready(code_toggle);\n",
       "</script>\n",
       "Some hidden code... to show/hide this cell's raw code input, click (or double click if stubborn!) <a href=\"javascript:code_toggle()\">here</a>."
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<style>\n",
       "div.warn {    \n",
       "    background-color: #b3b3ff;\n",
       "    border-color: #0000b3;\n",
       "    border-left: 5px solid #0000b3;\n",
       "    padding: 0.5em;\n",
       "    }\n",
       " </style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "from IPython.display import display\n",
    "\n",
    "\n",
    "# Taken from https://stackoverflow.com/questions/31517194/how-to-hide-one-specific-cell-input-or-output-in-ipython-notebook\n",
    "tag = HTML('''<script>\n",
    "code_show=true; \n",
    "function code_toggle() {\n",
    "    if (code_show){\n",
    "        $('div.cell.code_cell.rendered.selected div.input').hide();\n",
    "    } else {\n",
    "        $('div.cell.code_cell.rendered.selected div.input').show();\n",
    "    }\n",
    "    code_show = !code_show\n",
    "} \n",
    "$( document ).ready(code_toggle);\n",
    "</script>\n",
    "Some hidden code... to show/hide this cell's raw code input, click (or double click if stubborn!) <a href=\"javascript:code_toggle()\">here</a>.''')\n",
    "display(tag)\n",
    "\n",
    "############### Write code below ##################\n",
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./ES_ML_SOSSO/colors.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CH413: Computational Workshop 2 - Assignment\n",
    "## Global optimisation: clusters of 12, 13 and 14 atoms with the Morse potential\n",
    "\n",
    "Your task is to perform global optimisation - as on the workshop -, to find the global minimum structure of 12, 13 and 14 atom clusters, described by the Morse potential, and compare the per atom energy of the three global minimum structures. The notebook we used during the workshop is a good starting point! The Morse potential is also included in the atomic simulation environment (ase), so we can use it as easily as the Lennard-Jones.\n",
    "\n",
    "For small cluster sizes the behaviour of the Morse potential is very similar to that of the Lennard-Jones potential. By using the default potential parameters in ase, the energy and distance scale will be very similar to the one in the Lennard-Jones system we used at the workshop. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=warn>\n",
    "<h2>Assessment</h2>\n",
    "    \n",
    "* <b>Step 1 [5%]</b>: Create a markdown cell, and with using latex provide the equation for the Morse potential. Define all symbols involved.  <br>\n",
    "* <b>Step 2 [5%]</b>: Import the necessary libraries and initialise a calculator for the Morse potential. (Use the default parameters, thus no need to assign or define any arguments for the potential.)<br>\n",
    "* <b>Step 3 [10%]</b>: Create one example initial structure for all three cluster sizes using a cubic simulation cell with boxlength of 15 Angstrom and enabling periodic boundary conditions. (Think about how and where you place the atoms and visualise the structure to check the initial structure is sensible.)<br>\n",
    "* <b>Step 4 [10%]</b>: Perform a test minimisation on the largest cluster to determine the parameters (``fmax`` and ``steps``) for an appropriate level of convergence.<br>\n",
    "* <b>Step 5 [25%]</b>: Perform 40 minimisation cycles for all three cluster sizes and save the minimum energies and structures in separate arrays. <br>\n",
    "* <b>Step 6 [15%]</b>: Create a scatter plot showing the energy of the minimised configurations for all three cluster sizes. (Adjust the energy axis so data for all three clusters are visible, but do not show the minimisations that failed and got stuck at really high energies.) Save the plot in png format.<br>\n",
    "* <b>Step 7 [10%]</b>: Choose the lowest energy configuration of all three cluster sizes and visualize them to check the structure. Create a markdown cell and briefly comment on your findings (describe briefly the structures and how does these compare to known global minimum structures of Morse clusters)<br>\n",
    "* <b>Step 8 [10%]</b>: Once you are satisfied that the global minima has been found for all three cluster sizes, create a plot showing the energy of the global minimum per atom for the three global minima as a function of the number of atoms in the cluster. Save the plot in png format.<br>\n",
    "* <b>Step 9 [10%]</b>: What is the trend for the energy/atom for clusters of 12, 13 and 14 atoms? Give a short explanation of your findings in a markdown cell.<br>   \n",
    "\n",
    "<div/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the equation for the Morse potential: $$V'(r)=D_e(1-e^{-\\alpha(r-r_e)})^2$$\n",
    "\n",
    "where $r$ is the distance between the atoms, $r_{e}$ is the equilibrium bond distance, $D_e$ is well depth, and $\\alpha$ controls the width of the well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code imports all the packages needed to execeute the rest of the workbook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import numpy as np # numpy = numeric python. Useful when dealing with e.g. arrays...\n",
    "# import matplotlib.pyplot as plt # useful to plot data\n",
    "# import ase, ase.io # ase=atomic simulation environment. A MUST for every computational chemist...\n",
    "# from ase.visualize import view # to allow us to creat a 3D snapshot of the atomic structures\n",
    "# from ase.optimize import BFGS # Broyden–Fletcher–Goldfarb–Shanno minimisation algorithm\n",
    "# from ase.calculators.morse import MorsePotential # built-in Morse calculator in ase\n",
    "# calc=MorsePotential() # initialising calculator for Morse Potential\n",
    "# import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cells visualise an exmaple initial structure for 12, 13 and 14 atom clusters using a cubic box of 15 Angstrom, with the positions randomly distributed between 4.5 Angstroms and 9.75 Angstroms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster of 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# N_atoms=12\n",
    "# cell=np.eye(3)*15.0\n",
    "# at_12Ne=ase.Atoms(pbc=[(True,True,True)],cell=cell)\n",
    "\n",
    "# for i in range(N_atoms): \n",
    "#     pos = np.random.rand(3)*(0.3*cell[1,1])+0.35*cell[1,1] \n",
    "#     at_12Ne.append(ase.Atom(\"Ne\",position=pos)) \n",
    "    \n",
    "# view(at_12Ne,viewer='x3d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster of 13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# N_atoms=13\n",
    "# cell=np.eye(3)*15.0\n",
    "# at_13Ne=ase.Atoms(pbc=[(True,True,True)],cell=cell)\n",
    "\n",
    "# for i in range(N_atoms): \n",
    "#     pos = np.random.rand(3)*(0.3*cell[1,1])+0.35*cell[1,1] \n",
    "#     at_13Ne.append(ase.Atom(\"Ne\",position=pos)) \n",
    "    \n",
    "# view(at_13Ne,viewer='x3d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster of 14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# N_atoms=14\n",
    "# cell=np.eye(3)*15.0\n",
    "# at_14Ne=ase.Atoms(pbc=[(True,True,True)],cell=cell)\n",
    "\n",
    "# for i in range(N_atoms): \n",
    "#     pos = np.random.rand(3)*(0.3*cell[1,1])+0.35*cell[1,1] \n",
    "#     at_14Ne.append(ase.Atom(\"Ne\",position=pos)) \n",
    "    \n",
    "# view(at_14Ne,viewer='x3d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These lines of code are to run a minimisation on the 14 atom cluster to assess the best values of fmax, the number of steps, and to view the resulting structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# at_14Ne.set_calculator(calc)\n",
    "# dyn = BFGS(atoms=at_14Ne, trajectory='Ne14e_bfgs.traj')\n",
    "# dyn.run(fmax=0.001,steps=1000)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view(at_14Ne,viewer='x3d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 5+6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following lines of code run 40 minimisations for the 3 atom clusters. The array of minimum energies and the lowest of the array is then printed. A plot showing the minimum energies of all 3 atom clusters is then produced."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster of 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# minima_energy_M12=[] \n",
    "# minima_structure_M12=[] \n",
    "# for j in range(40):\n",
    "    \n",
    "#     N_atoms=12\n",
    "#     cell=np.eye(3)*15.0\n",
    "#     at_12Ne=ase.Atoms(pbc=[(True,True,True)],cell=cell)\n",
    "\n",
    "#     for i in range(N_atoms): \n",
    "#         pos = np.random.rand(3)*(0.3*cell[1,1])+0.35*cell[1,1] \n",
    "#         at_12Ne.append(ase.Atom(\"Ne\",position=pos)) \n",
    "    \n",
    "#     at_12Ne.set_calculator(calc)\n",
    "#     dyn = BFGS(atoms=at_12Ne, trajectory='Ne12_bfgs.traj')\n",
    "#     dyn.run(fmax=0.001,steps=1000)\n",
    "    \n",
    "#     e = at_12Ne.get_potential_energy() \n",
    "#     minima_energy_M12.append(e) \n",
    "#     minima_structure_M12.append(at_12Ne) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(minima_energy_M12)\n",
    "# print(min(minima_energy_M12))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster of 13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# minima_energy_M13=[] \n",
    "# minima_structure_M13=[] \n",
    "\n",
    "# for j in range(40):\n",
    "    \n",
    "#     N_atoms=13\n",
    "#     cell=np.eye(3)*15.0\n",
    "#     at_13Ne=ase.Atoms(pbc=[(True,True,True)],cell=cell)\n",
    "\n",
    "#     for i in range(N_atoms): \n",
    "#         pos = np.random.rand(3)*(0.3*cell[1,1])+0.35*cell[1,1] \n",
    "#         at_13Ne.append(ase.Atom(\"Ne\",position=pos)) \n",
    "    \n",
    "#     at_13Ne.set_calculator(calc)\n",
    "#     dyn = BFGS(atoms=at_13Ne, trajectory='Ne13_bfgs.traj')\n",
    "#     dyn.run(fmax=0.001,steps=1000)\n",
    "    \n",
    "#     e = at_13Ne.get_potential_energy() \n",
    "#     minima_energy_M13.append(e) \n",
    "#     minima_structure_M13.append(at_13Ne) \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(minima_energy_M13)\n",
    "# print(min(minima_energy_M13))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster of 14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# minima_energy_M14=[] \n",
    "# minima_structure_M14=[] \n",
    "# for j in range(40):\n",
    "    \n",
    "#     N_atoms=14\n",
    "#     cell=np.eye(3)*15.0\n",
    "#     at_14Ne=ase.Atoms(pbc=[(True,True,True)],cell=cell)\n",
    "\n",
    "#     for i in range(N_atoms): \n",
    "#         pos = np.random.rand(3)*(0.3*cell[1,1])+0.35*cell[1,1] \n",
    "#         at_14Ne.append(ase.Atom(\"Ne\",position=pos)) \n",
    "    \n",
    "#     at_14Ne.set_calculator(calc)\n",
    "#     dyn = BFGS(atoms=at_14Ne, trajectory='Ne14_bfgs.traj')\n",
    "#     dyn.run(fmax=0.001,steps=1000)\n",
    "    \n",
    "#     e = at_14Ne.get_potential_energy() \n",
    "#     minima_energy_M14.append(e) \n",
    "#     minima_structure_M14.append(at_14Ne) \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(minima_energy_M14)\n",
    "# print(min(minima_energy_M14))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# fig=plt.figure(num=None,figsize=(4,2),dpi=500,facecolor='w',edgecolor='k') \n",
    "# plt.tick_params(axis='both', which='major', labelsize=4) \n",
    "# plt.xlabel('Number of minimised structures',fontsize=6) \n",
    "# plt.ylabel('Energy of the minimised structure / eV',fontsize=6) \n",
    "# plt.ylim(-50.0, -10.0) \n",
    "\n",
    "# plot1=plt.scatter(np.arange(len(minima_energy_M12)),minima_energy_M12,color='rebeccapurple',s=1,label='12 atom cluster') # data(x), data(y), colour, size\n",
    "# plot2=plt.scatter(np.arange(len(minima_energy_M13)),minima_energy_M13,color='red',s=1,label='13 atom cluster')\n",
    "# plot3=plt.scatter(np.arange(len(minima_energy_M14)),minima_energy_M14,color='yellow',s=1,label='14 atom cluster')\n",
    "# plt.legend(prop={'size': 3})\n",
    "# plt.savefig('Morse_min.png'); "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 7\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The minimum energy structures for the 3 cluster sizes are visualised, and a markdown cell shows the energy values of the global minima, which are then visualised/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view(minima_structure_M12[minima_energy_M12.index(min(minima_energy_M12))]),viewer='x3d')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view(minima_structure_M13[minima_energy_M13.index(min(minima_energy_M13))],viewer='x3d')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view(minima_structure_M14[minima_energy_M14.index(min(minima_energy_M14))],viewer='x3d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The global minimum energies found via the minimisations have very similar, almost the same energies as the global minima. It would be sensible to say that if the energies are similar, they would have similar, if not the same, structures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The global minima for clusters of 12, 13 and 14 atoms, with a cutoff of 6.0 are as follows:\n",
    "\n",
    "Cluster of 12: 12B, with a global minimum of -36.400278\n",
    "\n",
    "Cluster of 13: 13A, with a global minimum of -42.439863\n",
    "\n",
    "Cluster of 14: 14B with a global minimum of -45.619277\n",
    "\n",
    "Values taken from: http://doye.chem.ox.ac.uk/jon/structures/Morse/tables.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CELD_12=pd.read_table(\"http://doye.chem.ox.ac.uk/jon/structures/Morse/points/12B\", sep = \"\\s+\", skiprows = 1, header = None)\n",
    "# at_12_minimum=ase.Atoms('Ne12',\n",
    "#                        positions=CELD_12.values,\n",
    "#                        cell=[(10,0,0),(0,10,0),(0,0,10)],\n",
    "#                        pbc=[(True, True, True)])\n",
    "# at_12_minimum.center()\n",
    "# view(at_12_minimum, viewer='x3d')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CELD_13=pd.read_table(\"http://doye.chem.ox.ac.uk/jon/structures/Morse/points/13A\", sep = \"\\s+\", skiprows = 1, header = None)\n",
    "# at_13_minimum=ase.Atoms('Ne13',\n",
    "#                        positions=CELD_13.values,\n",
    "#                        cell=[(10,0,0),(0,10,0),(0,0,10)],\n",
    "#                        pbc=[(True, True, True)])\n",
    "# at_13_minimum.center()\n",
    "# view(at_13_minimum, viewer='x3d')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CELD_14=pd.read_table(\"http://doye.chem.ox.ac.uk/jon/structures/Morse/points/14B\", sep = \"\\s+\", skiprows = 1, header = None)\n",
    "# at_14_minimum=ase.Atoms('Ne14',\n",
    "#                        positions=CELD_14.values,\n",
    "#                        cell=[(10,0,0),(0,10,0),(0,0,10)],\n",
    "#                        pbc=[(True, True, True)])\n",
    "# at_14_minimum.center()\n",
    "# view(at_14_minimum, viewer='x3d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 8+9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two arrays are generated, one for the minimum energies for the 3 cluster sizes, the other for the number of atoms in the cluster. A plot of energy/atom is then produced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# minimum_energy_per_atom = [(min(minima_energy_M12))/12,\n",
    "#                            (min(minima_energy_M13))/13,\n",
    "#                            (min(minima_energy_M14))/14]\n",
    "# no_of_atoms = [12,13,14]\n",
    "\n",
    "# print(minimum_energy_per_atom)\n",
    "# print(no_of_atoms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# xdata=no_of_atoms\n",
    "# ydata=minimum_energy_per_atom\n",
    "# plt.plot(xdata,ydata)\n",
    "# plt.ylabel('Energy per Atom (eV)')\n",
    "# plt.xlabel('Number of Atoms in Cluster')\n",
    "# plt.savefig(\"Energy_per_atom.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comment on the trend of energy/atom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The structure of a 12 atom cluster is reasonably stable, which leads to a low value of energy/atom. There is an empty space which is filled for a 13 atom cluster, causing a sharp decrease in the energy/atom value and the most stable strcture. For the 14 atom cluster, the energy/atom value increases a little because of the extra atom causing a less stable structure."
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}