Workshop_2_2020-2021.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CH413: Computational Workshop 2\n",
    "## Global optimisation of clusters\n",
    "\n",
    "In this workshop, we are going to use simple minimisation to find local minima of the 13 atom Lennard-Jones cluster, starting from random initial configurations. This is like the AIRSS (ab initio random structure search) technique we were learning about, but with using the Lennard-Jones potential to describe the interactions instead of ab inito energy calculations, as that would be way too expensive (= would take too long).\n",
    "\n",
    "We will use the same environment as in Workshop 1: Jupyter notebook and python. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='background:lightblue'>\n",
    "As before, these blue boxes contain instructions for you to follow, or stuff for you to do\n",
    "<h2>How to access this Jupyter notebook</h2>\n",
    "\n",
    "* <b>Step 1</b>: Open a web browser and go to [this page](https://mnf145.csc.warwick.ac.uk:8987/module/CH413). <br>\n",
    "* <b>Step 2</b>: The Notebook Launcher pops up: select the CH413 module and fill in the boxes using your SCRTP username and password <br>\n",
    "* <b>Step 3</b>: Open the Jupyter notebook you are interested in, e.g. Workshop_2_Class.ipynb <br>\n",
    "* <b>Step 4</b>: Make a copy of the orginal notebook (which is read only). In the toolbar on top of the notebook: File - Make a Copy (you can also rename it)<br>\n",
    "* <b>Step 5</b>: You're all set! <br><br>\n",
    "<b> Remember: </b> You can access this notebook at home at any time by going through the same steps on e.g. your laptop - all the changes you have made will be saved and synced! <br>\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='background:lightblue'>\n",
    "Let's warm up by writing a bit in a new Markdown cell. Since we are going to use the Lennard-Jones potential to describe the interactions between atoms, let's remind ourselves of the functional form: use LaTeX to write down the expression for the potential. (Remember: LaTeX is fun and Google is your friend!)\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the Lennard-Jones pair potential: $$U_\\mathrm{LJ}(r)=abc?$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import libraries\n",
    "Before we are able to do anything else, we need to import some libraries to help us define and manipulate coordinates, plot data, perform minimisations, visualise atomic structures...etc. Central to this workshop is the ase package. It stands for Atomic Simulation Environment (ASE), and it is a set of tools and Python modules for setting up, manipulating, running, visualizing and analyzing atomistic simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
    "import nglview as nv # visualisation\n",
    "import pytraj as pt # allowing handling and visualisation of trajectories\n",
    "from ase.optimize import BFGS # Broyden–Fletcher–Goldfarb–Shanno minimisation algorithm\n",
    "from ase.calculators.lj import LennardJones # built-in LJ calculator in ase"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  The Atoms Object\n",
    "In order to describe the atomic system (the cluster in our case), we need to define what kind of atoms build it up, how many of them, what are their coordintes and how big is the cell. Ase allows us to store this information in a very convenient way: the Atoms object! The Atoms object has a lot of possible keyword arguments, e.g.:\n",
    "\n",
    "- ``positions`` stores the coordinates\n",
    "- ``cell`` stores the simulation cell parameters in a 3*3 matrix \n",
    "- ``forces`` stores the forces on every single atoms\n",
    "- ``pbc`` determines if periodic boundary conditions are applied\n",
    "\n",
    "Let's look at an example of 3 nitrogen atoms positioned along a line, each at 1.5 Angstrom distance, in a cubic box with length of 15 Angstrom. The Atoms object for this system looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "at_3N = ase.Atoms('O3',  # three nitrogen atoms\n",
    "                  positions=[(5,5,5),(3.5,5,5),(6.5,5,5)], # xyz coordinates for all three atoms\n",
    "                  cell=[(15,0,0),(0,15,0),(0,0,15)], # three cell vectors in a matrix\n",
    "                  pbc=[(True,True,True)]) # applying periodic boundary conditions for all directions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualisation\n",
    "It's always a good idea to visualise the atomic structures we are working with. Not just because they look really cool, but a quick glance can give us a lot of information on what's going on! E.g. to see whether there're some silly unphysical arrangements, quickly find some obvious symmetries in our system...etc.\n",
    "\n",
    "We will use the nglview visualisation package, with which you can \n",
    "* rotate the structure, \n",
    "* zoom in, \n",
    "\n",
    "The colour and the size of the sphere the atom is represented with depends on the type of atom set in the Atoms object, e.g. nitrogen is dark blue, oxygen is red...etc., and bonds are automatically reresented if two atoms are closer than the typical bonding distance of the given atom type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = nv.show_ase(at_3N)\n",
    "view.add_unitcell()\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4> Create an Atoms object from scratch</h4>\n",
    "    \n",
    "<p style='background:lightblue'> In the cell below, create and visualise an Atoms object which contains 4 neon atoms in a cubic box with length of 10 Angstrom. Place the atoms such that they form a square, nearest neighbours placed at about 2 Angstrom distance. Set periodic boundary conditions True.     \n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create your Atoms object here\n",
    "at_4Ne=\n",
    "\n",
    "view = nv.show_ase(at_4Ne)\n",
    "view.add_unitcell()\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Minimisation\n",
    "It's time to do some minimisations, so we will take the 4 Ne atom configurations we just created and see what is the nearest energy minimum. \n",
    "\n",
    "We are going to use the Lennard-Jones potential to describe the interactions, so we will need to set the Lennard-Jones parameters and initialise the calculator within ase, which will take care of all the energy and force calculations for us. \n",
    "(A quick note: it doesn't really matter whether we call our atoms neon atoms or else. I suggest using Ne, as it will be shown with a sphere size of about the right size to see the structures the best, but if you'd like to use something else, that's fine, it won't make any difference to numerical results, as long as the LJ paramteres are set as below.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='background:lightblue'> Set the LJ parameters $\\sigma$ and $\\epsilon$ to 1.0, the potential truncation cutoff (rc) to 5.0.\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma=  \n",
    "epsilon=  \n",
    "rc=  \n",
    "calc = LennardJones(sigma=sigma, epsilon=epsilon, rc=rc) # We need to set this once, then we can use this ``calc`` again\n",
    "at_4Ne.set_calculator(calc) # we need to set the calculator for every Atoms object we want to work with"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<h4> Set up and perform the minimisation</h4>\n",
    "\n",
    "We will use the BFGS (Broyden-Fletcher-Goldfarb-Shanno) algorithm to perform the minimisation (this is built in in ase). ``fmax`` and ``steps`` control when the minimisation is stopped: either when the largest force on any of the atoms becomes smaller than ``fmax``, or when the number of iteration reaches ``steps``, whichever is sooner. The minimisation will print out the energy and the largest force at every iteration step, and the Atoms object will be updated with the final coordinates. <br>\n",
    "<p style='background:lightblue'> \n",
    "* How does the energy and the force changes during the minimisation? <br></p>\n",
    "<p style='background:lightblue'>\n",
    "* Try decreasing the ``fmax`` and increasing the ``steps``, and rerun the minimisation. (It will continue from the last previous step). What happens? <br></p>\n",
    "<p style='background:lightblue'>\n",
    "* Visualise the structure, how did it change from the initial configuration? <br></p>\n",
    "<p style='background:lightblue'>\n",
    "* Go back to the cell where the initial atoms object was created, and modify the structure such that the atoms are not perfectly in plane. Repeat the minimisation. How does the final structure and energy compares to the previous minimisation? Can you explain the difference <br></p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dyn = BFGS(atoms=at_4Ne, trajectory='Ne4_square_bfgs.traj') # set up BFGS minimiser for the desired Atoms object\n",
    "dyn.run(fmax=0.01,steps=25) # perform minimisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualising the final structure\n",
    "view = nv.show_ase(at_4Ne)\n",
    "view.add_unitcell()\n",
    "view"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ne4_square_traj=ase.io.trajectory.TrajectoryReader(\"Ne4_square_bfgs.traj\") # read the trajectory to an atoms object\n",
    "ase.io.write(\"Ne4_square_bfgs.traj.pdb\",Ne4_square_traj) # write a pdb formatted file which can be handled by pytraj\n",
    "p_traj = pt.load(\"Ne4_square_bfgs.traj.pdb\")\n",
    "p_view = nv.show_pytraj(p_traj)\n",
    "p_view.add_unitcell()\n",
    "p_view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding the global minimum of LJ$_{13}$\n",
    "\n",
    "Now we are going to try to find the global minimum structure of 13 Lennard-Jones atoms starting from random configurations. We will use the same cell as before, but this time placing 13 atoms randomly inside, doing this with a cycle and a random number generator, instead of doing it by hand. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_atoms=13 # number of atoms to be included\n",
    "cell=np.eye(3)*10.0 # creating a 3*3 matrix with 10.0 in the diagonal, zeros elsewhere. this is our simulation box.\n",
    "at_13Ne=ase.Atoms(pbc=[(True,True,True)],cell=cell) # create the Atoms object, but do not add any atoms yet\n",
    "\n",
    "for i in range(N_atoms): # this is a for cycle repeated N_atoms times \n",
    "    pos = np.random.rand(3)*cell[1,1] # generate three uniform random real number in the range [0,cell[1,1]=10]\n",
    "    at_13Ne.append(ase.Atom(\"Ne\",position=pos)) # add these as the positions of the atom added to the Atoms Object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='background:lightblue'> \n",
    "Check the coordinate values that were generated and visualise the random structure. \n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(at_13Ne.positions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = nv.show_ase(at_13Ne)\n",
    "view.add_unitcell()\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='background:lightblue'>    \n",
    "* Perform the minimisation of this 13 atom random structure! Don't forget that you will need to set the calculator for this Atoms object as well. Adjust the ``fmax`` and ``steps`` parameters to achieve convergence.\n",
    "</p>    \n",
    "<p style='background:lightblue'>\n",
    "* Visualise the structure. What has happened? Has one single cluster been formed? \n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up and perform the minimisation\n",
    "at_13Ne.set_calculator(calc)\n",
    "dyn = BFGS(atoms=at_13Ne, trajectory='Ne13_bfgs.traj') # set up BFGS minimiser for the desired Atoms object\n",
    "dyn.run(fmax=0.01,steps=520) # perform minimisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = nv.show_ase(at_13Ne)\n",
    "view.add_unitcell()\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a good chance that multiple smaller clusters were formed during the minimisation, instead of one single cluster of all 13 atoms. This demonstrates that starting from a completely random structure is not always fruitful. E.g. if initially some of the atoms are further away than the cutoff distance of the potential, they won't 'feel' each other. To overcome this problem, we need to start from a more 'sensible' initial structure, thus making sure that none of the atoms are too far away."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='background:lightblue'>   \n",
    "Create another Atoms object, in which all 13 Ne atoms are placed randomly, but not within the entire simulation cell, just within a 5 Angstrom range in the middle of the box (thus all coordinates need to be generated within the range $[2.5,7.5]$). Type in a new formula to achieve this. (It's a good idea to visualise this new structure as well, to be sure all is fine.)\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_atoms=13 \n",
    "cell=np.eye(3)*10.0 \n",
    "at_13Ne_B=ase.Atoms(pbc=[(True,True,True)],cell=cell) \n",
    "\n",
    "for i in range(N_atoms): \n",
    "    pos = # fill in with new formula for coordinates HERE!\n",
    "    at_13Ne_B.append(ase.Atom(\"Ne\",position=pos)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = nv.show_ase(at_13Ne_B)\n",
    "view.add_unitcell()\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='background:lightblue'>   \n",
    "* Now perform the minimisation of this 13 atom random structure! Don't forget that you will need to set the calculator for this Atoms object as well. Adjust the ``fmax`` and ``steps`` parameters if necessary!!!</p>\n",
    "<p style='background:lightblue'>\n",
    "* Visualise the structure. Has a single cluster been formed?</p>\n",
    "<p style='background:lightblue'>\n",
    "* Is this the global minimum for LJ$_{13}$? (Use Google to find out the global minimum looks like.)\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up and perform the minimisation\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = nv.show_ase(at_13Ne_B)\n",
    "view.add_unitcell()\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform a series of minimisations and collect data ###\n",
    "Now we have set up everything to search for minima of the LJ$_{13}$ cluster and hopefully find the global minimum as well. We will repeat the minimisation (as we've done before) multiple times, using a ``for`` cycle,  always starting from a different initial random structure. We will collect the minimum energies and the minimum configurations to be able to examine the result. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='background:lightblue'>\n",
    "Use the cell below to set up a series of minimisations - first we will start with 10 cycles. Fill in the gaps by copying the necessary bits from previous cells. Don't forget that indentation does matter in python!  \n",
    "Run the cell. (It might take a minute or two to have this finished.)\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "minima_energy_LJ13=[] # the minimum energies will be collected in this array\n",
    "minima_structure_LJ13=[] # the minimised structures will be collected in this array (every element will be an Atoms object)\n",
    "\n",
    "# if you want to run this cell again to perform more minimisations, comment out the two lines above, otherwise\n",
    "# the arrays will be initialised again and data from the first round of calculation will be lost.\n",
    "\n",
    "for j in range(10):\n",
    "       \n",
    "    # create the initial structure here\n",
    "    \n",
    "    \n",
    "    # set calculator here\n",
    "    \n",
    "    \n",
    "    # set up and run the minimisation here\n",
    "    \n",
    "\n",
    "    \n",
    "    e = at_13Ne_B.get_potential_energy() # get the potential energy of the Atoms object (this is the final minimum energy)\n",
    "    minima_energy_LJ13.append(e) # add the energy value to the array\n",
    "    minima_structure_LJ13.append(at_13Ne_B) # add the Atoms obejct (structure) to the array\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4> Plot the minimum energies</h4>\n",
    "\n",
    "<p style='background:lightblue'>\n",
    "Let's plot the energies of the minimised configurations, using matplotlib. There might be a few random structures that could not be properly minimised (e.g. the two atoms were overlapping too much in the initial structure, so the forces were too high and the algorithm got stuck), if so, you will need to adjust the range of the $x$ axis. To find out what is a good range, print the ``minima_energy_LJ13`` array first in a separate cell.\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig=plt.figure(num=None,figsize=(4,2),dpi=500,facecolor='w',edgecolor='k') # create the figure\n",
    "plt.tick_params(axis='both', which='major', labelsize=4) # define font size for the tick labels\n",
    "plt.xlabel('number of minimised structures',fontsize=6) # set label of the x axis\n",
    "plt.ylabel('energy of the minimised structure',fontsize=6) # set label of the y axis\n",
    "#plt.ylim(-100.0, +100.0) # set limits for the y axis (bottom,top) - adjust the values!!!\n",
    "plt.scatter(np.arange(len(minima_energy_LJ13)),minima_energy_LJ13,color='rebeccapurple',s=1) # data(x), data(y), colour, size\n",
    "plt.savefig('LJ13_min.png'); # save the plot as a figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=warn>\n",
    "<h4> Have you found the global minimum?</h4>\n",
    "\n",
    "Check out the Cambridge Cluster Database (CCD) to find out what is the global minimum energy.\n",
    "* What is the lowest energy value you have got? \n",
    "* What is the corresponding structure? (Visualise the corresponding element of the ``minima_structure_LJ13`` array.)  Is this the global minimum?  (If the energy value differs only in the third or fourth digit from the value in CCD, do not worry, such small differences can arise if the cutoff distance is different.) \n",
    "\n",
    "\n",
    "If none of the structures were the global minimum, continue the minimisations.\n",
    "<div/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The lowest energy value I've got:\n",
    "\n",
    "Is this the known global minimum? Yes/No"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=warn>\n",
    "\n",
    "<h4> If you're done with the tasks above...</h4>\n",
    "If you have done all the tasks in the notebook so far, and would like to experiment a bit further, here are a couple of further tasks/questions you can work one.\n",
    "\n",
    "* What is the second lowest energy structure? \n",
    "* Look at some of the lowest energy local minima, how different are they? \n",
    "\n",
    "If you are more familiar with python:\n",
    "\n",
    "* some minimisations gets stuck, because the initial structure contains atoms which overlap. Modify the code such that these configurtions are discarded and a new initial structure is generated to replace them (hint: e.g. use the initial energy value to determine which ones should be discarded.) \n",
    "* if you have at least 50-60 minimisations, produce a histogram showing how many times did you find the different local minima\n",
    "  \n",
    "<div/>"
   ]
  }
 ],
 "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": 2
}