05_Running_DFT_calculation_on_PACE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running a simple DFT calculation on PACE Supercomputer Cluster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this lecture, we will learn about running DFT calculations on PACE. We will first learn about PBS submission scripts that are used to submit jobs on PACE. As a member of Medford group, you will have access to two supercomputer clusters here at Georgia Tech: PACE-Phoenix and PACE-Hive. If you are part of the VIP course, you will have access to PACE-ICE cluster. If you are on PACE-Phoenix, you will need to request for resources from your PI's account. This means that you are paying for the jobs you are running therefore, these resources should be used judiciously. An example run.sh script to submit a QE job is given below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/bin/bash\n",
    "#PBS -l nodes=1:ppn=8\n",
    "#PBS -l pmem=3GB\n",
    "#PBS -l walltime=12:00:00\n",
    "#PBS -A GT-amedford6\n",
    "#PBS -N optimizer\n",
    "#PBS -o stdout\n",
    "#PBS -e stderr\n",
    "\n",
    "cd $PBS_O_WORKDIR\n",
    "\n",
    "source ~/.bashrc\n",
    "\n",
    "module load anaconda3 #It is optional to load anaconda. You don't need to load it if it's already been loaded in the environment \n",
    "\n",
    "source /storage/coda1/p-amedford6/0/shared/rich_project_chbe-medford/medford-share/envs/espresso-6.7MaX-beef-ipi\n",
    "\n",
    "python calc.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To submit a job to the PBS system, you have to use the **qsub** command. We usually call this file run.sh for convenience. The command for submiiting the job will be **qsub run.sh**. When you do this, you are essentially handing the file to the PBS queuing system which will scan the file for the blocks containing #PBS that specify the resources you are requesting for the job. Each line in this file starting with #PBS specifies something different. The line with tag *-l* are requesting certain resources (nodes and processors per node). You can request the resources based on requirement of your calculations. In this example script, we are requesting for 1 node and 8 processors per node. The next line (#PBS -l pmem=3GB) specifies the memory you request per processor. Here it is 3 GB. Hence the total memory we request for is 24 GB. For Phoenix cluster, the default queue is *inferno* where the jobs are run. -N is the name you want to give the job, -o and -e are the filenames where the outputs and errors will be printed respectively. You can run *pace-whoami* check the maximum resources available for a particular queue. \n",
    "\n",
    "If you are submitting a job to PACE-HIVE, you will need to specify the queue using **#PBS -q**. HIVE has different queues like hive, hive-himem, hive-interact which offer different resources and walltime. After the PBS block, you can source your environment and run the code. The next line with change directory (cd) command with an alias for current working directory (PBS_O_WORKDIR) changes to the correct folder. More information about PBS script and submission of jobs on PACE can be found in this [link](https://pace.gatech.edu/pace-user-documentation). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running Quantum Espresso on PACE Phoenix\n",
    "\n",
    "For the purpose of DFT calculations, our group primarily uses Quantum Espresso (QE) which is a plane-wave method code.  Various versions for QE are present in the shared folder of Medford group. The environment can be sourced in the bash file. An example command to source QE environment is:\n",
    "\n",
    "**source /storage/coda1/p-amedford6/0/shared/rich_project_chbe-medford/medford-share/envs/espresso-6.7Max-beef-ipi**\n",
    "\n",
    "Version 6.7Max is the environment for latest compiled QE code. We will first go through the details a simple code for running a DFT calculation using Quantum Espresso. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#example code for running Quantum Espresso\n",
    "\n",
    "#import all the necessary modules\n",
    "from ase import Atom, Atoms\n",
    "from ase.optimize import BFGSLineSearch as BFGS\n",
    "from ase.build import molecule\n",
    "from espresso import iEspresso as Espresso\n",
    "from ase.io import read, write\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "t0 = time.time()\n",
    "images = []\n",
    "\n",
    "#building the \"molecule\" object\n",
    "image = molecule('H2O',vacuum=10.) \n",
    "image.set_cell([10, 10, 10]) \n",
    "image.set_pbc([1,1,1])\n",
    "image.center()\n",
    "image.rattle(0.0001)\n",
    "\n",
    "#Use QE as a single-point calculator (just SCF)\n",
    "calc = Espresso(atoms=image,\n",
    "                pw=500.0,\n",
    "                xc='PBE',\n",
    "                kpts=\"gamma\")\n",
    "# retrieve Energy, force, stress\n",
    "\n",
    "print('1st calculation with PBE')\n",
    "print('Energy = {}'.format(image.get_potential_energy()))\n",
    "print('Force = {}'.format(image.get_forces()))\n",
    "print('Stress = {}'.format(image.get_stress()))\n",
    "t1 = time.time()\n",
    "print('Time = {} sec'.format(t1-t0))\n",
    "print('-'*20)\n",
    "\n",
    "# Run optimization (will change the atomic positions as it is relaxation calculation)\n",
    "dyn = BFGS(image,logfile='opt_pbe.log', trajectory='h2o_pbe_optimization.traj')\n",
    "dyn.run(fmax=0.005)\n",
    "image.calc.close()\n",
    "\n",
    "# Attach BEEF-vdw calculator\n",
    "t2 = time.time()\n",
    "calc = Espresso(atoms=image,\n",
    "                pw=350.0,\n",
    "                xc='BEEF-vdw')\n",
    "\n",
    "# retrieve Energy, force, stress after optimization\n",
    "print('2nd calculation with BEEF-vdw')\n",
    "print('Energy = {}'.format(image.get_potential_energy()))\n",
    "print('Force = {}'.format(image.get_forces()))\n",
    "print('Stress = {}'.format(image.get_stress()))\n",
    "\n",
    "E_ensemble = image.get_ensemble_energies()\n",
    "print('Ensemble energies = {} +/- {}'.format(E_ensemble.mean(), E_ensemble.std()))\n",
    "print('Time = {} sec'.format(time.time()-t2))\n",
    "\n",
    "# Run optimization with BEEF\n",
    "dyn = BFGS(image,logfile='opt_beef.log',trajectory='h2o_beef_optimization.traj')\n",
    "dyn.run(fmax=0.005)\n",
    "image.calc.close()\n",
    "\n",
    "print('Totol time = {} sec'.format(time.time()-t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's take a look at this codeblock line by line to understand what is happening as we run this code. For setting up a DFT calculation, we first need to build the chemical system whose properties we want to calculate. In the example given, the object is a water molecule. By building the system, we are essentially specifying its atomic positions. Since it is a water molecule, it has to be put in a isolated box with 10 Angstrom vacuum in all three directions so that electron density of the molecule doesn't overlap with that of its images. The molecule should be located at the center of the box. We can verify this by using ASE's visualize method as we learnt in the previous lectures. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#building the \"molecule\" object\n",
    "image = molecule('H2O',vacuum=10.) \n",
    "image.set_cell([10, 10, 10]) \n",
    "image.set_pbc([1,1,1])\n",
    "image.center()\n",
    "image.rattle(0.0001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After defining the system, we will set up a calculator where we define the parameters for the DFT calculation. These parameters will vary for different DFT codes. The important parameters that you will need to specify for your system in Quantum Espresso calculator are the following:\n",
    "\n",
    "1. atoms: The Atoms object \n",
    "1. pw: Plane-wave cutoff\n",
    "2. kpts: K-points\n",
    "3. xc: Exchange-correlation functional\n",
    "4. spinpol: Spin-polarization\n",
    "5. dipole: Dipole-correction\n",
    "\n",
    "You will also need to decide if you want to run a single-point calculation or relaxation calculation. A single-point calculation might result in very high forces is most cases as the atoms in the initial guess might not be relaxed. The first step in the example script is to run a single-point calculation (also known as SCF calculation). SCF stands for self-consistent field. Note that the calculator also takes in the atoms object so that it knows where to get the data from. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calc = Espresso(atoms=image,\n",
    "                pw=500.0,\n",
    "                xc='PBE',\n",
    "                kpts=\"gamma\")\n",
    "\n",
    "print('1st calculation with PBE')\n",
    "print('Energy = {}'.format(image.get_potential_energy()))\n",
    "print('Force = {}'.format(image.get_forces()))\n",
    "print('Stress = {}'.format(image.get_stress()))\n",
    "t1 = time.time()\n",
    "print('Time = {} sec'.format(t1-t0))\n",
    "print('-'*20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some of the properties you want to obtain from a DFT calculation can be:\n",
    "\n",
    "1. Energy and forces\n",
    "2. Density of states \n",
    "3. Electron density\n",
    "4. Electrostatic potential\n",
    "\n",
    "The next step in the example script is to run an optimization that minimizes the forces on atoms in the system. This is also known as geometric optimization because the geometry of system changes such that forces on atoms are below a certain threshold which is user defined. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run optimization (will change the atomic positions as it is relaxation calculation)\n",
    "dyn = BFGS(image,logfile='opt_pbe.log', trajectory='h2o_pbe_optimization.traj')\n",
    "dyn.run(fmax=0.005)\n",
    "image.calc.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For optimization, we can use the optimizers available within ASE's implementation. The force threshold in the example is 0.005 eV/Angstrom. The file \"opt_pbe.log\" will save the steps of optimization. You can find out the force and total ground state energy in this file. The optimization trajectory is saved in the \"h2o_pbe_optimization.traj\" file. You can visualize it to observe how the atoms move to obtain their relaxed atomic positions. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next few lines in the example script show how to calculate energy, forces and stresses using the BEEF-vdw exchange correlation functional which is known to account for Van der waals interactions in chemical systems. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Strategy on selecting pseudopotentials for DFT calculations:\n",
    "\n",
    "For running DFT calculations, we will need to select pseudopotentials. We have a large pool of pseudopotentials in our shared folder on PACE-HIVE and PACE-PHOENIX. We also have a module for selecting the pseudopotentials which can be imported using the following command:\n",
    "\n",
    "from getpseudo import getpseudo\n",
    "\n",
    "By default, we use the standard solid-state pseudopotential (SSSP) for the calculations. In future, we will default to a set of pseudopotentials which is compatible with both SPARC and Quantum Espresso. To check the available set of pseudopotentials, source the environment QE or SPARC environment (we will learn about SPARC in the future lessons). The commands for sourcing the environments are the following:\n",
    "\n",
    "* For QE:\n",
    "\n",
    "**source /storage/coda1/p-amedford6/0/shared/rich_project_chbe-medford/medford-share/envs/espresso-6.7MaX-beef-ipi**\n",
    "\n",
    "* For SPARC:\n",
    "\n",
    "**source /storage/coda1/p-amedford6/0/shared/rich_project_chbe-medford/medford-share/envs/SPARC_latest_2021**\n",
    "\n",
    "After source the environment, you can use the command \"setpseudos\" on terminal which will give a list of available pseudopotentials in the following format:\n",
    "\n",
    "*esp#pseudos*\n",
    "\n",
    "*test_mostafa#210128#v1*\n",
    "\n",
    "*test_mostafa#psp#GGA#ultrasoft#210217#v0*\n",
    "\n",
    "*test_oncv#psp#pbe*\n",
    "\n",
    "*test_mostafa#210128#v1#UPF*\n",
    "\n",
    "*test_mostafa#210224#oncv#psps_soft*\n",
    "\n",
    "*test_mostafa#210224#oncv#psps_standard*\n",
    "\n",
    "*espresso_sg15_pp#brown#brenda*\n",
    "\n",
    "*espresso_sg15_sparc#qimen_SG15#pseudopotentials#psp8#upf_LDA#upf#cut*\n",
    "\n",
    "*espresso_sg15_sparc#qimen_SG15#pseudopotentials#psp8#upf_PBE#upf#cut*\n",
    "\n",
    "*espresso_sg15_sparc#qimen_SG15#pseudopotentials#psp8#upf_PBE#psp#cut*\n",
    "\n",
    "*espresso_sg15_sparc#qimen_SG15#pseudopotentials#psp8#upf_LDA#psp#cut*\n",
    "\n",
    "*espresso_sg15_rel*\n",
    "\n",
    "For setting the path to pseudopotential directory in a DFT calculation, you can import the **getpseudo** module and set the path as an input argument to the calculator. For example, in SPARC, you can set the calculator as given in the following example:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calc = SPARC(atoms = sim, pseudo_dir=getpseudo(\"espresso_sg15_sparc#qimen_SG15#pseudopotentials#psp8#upf_PBE#psp#cut\"), **parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where, **pseudo_dir** is used to set the path to pseudopotential. \"parameters\" refers to all other input arguments to the SPARC calculator. \n",
    "\n",
    "Similar to SPARC, for setting pseudopotential path in QE, you can refer to the following example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calc = Espresso(atoms=sim,\n",
    "                pw = 500,\n",
    "                psppath = getpseudo(\"espresso_sg15_mcsh#related#psps_mcsh#related#psps_QE\"),\n",
    "                convergence = {\"energy\": 0.0005, \\\n",
    "                               \"mixing\": 0.1, \\\n",
    "                               \"maxsteps\": 750, \\\n",
    "                               \"diag\": \"david\"},\n",
    "                smearing = \"gauss\",\n",
    "                sigma = 0.05,\n",
    "                nbands = -10,\n",
    "                xc = 'PBE',\n",
    "                kpts = (4,4,1),\n",
    "                spinpol = False,\n",
    "                dipole = {\"status\" : True},\n",
    "                nosym=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise:\n",
    "    \n",
    "1. Run a simple DFT calculation using QE on PACE-PHOENIX by following the steps given below:\n",
    "\n",
    "* Use ASE to generate a CO_${2}$ molecule. \n",
    "* Run a single point calculation by sourcing the QE environment using a PBS script and report the energy. Use PBE as the exchange correlation functional. \n",
    "* Run an optimization calculation for CO_$2$ molecule by refering to the example given in this lecture for H_${2}$O molecule and report the energy of the molecule, forces and stresses on each atom. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": 4
}