{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geometry optimization\n",
    "\n",
    "MyBigDFT comes with some classes implementing particular workflows of calculations. These workflows define a queue of jobs, that can easily be run sequentially, without having to worry about the Job context manager. They also generally define a particular post-processing procedure, run after all the BigDFT calculations in order to extract some meaningful imformation.\n",
    "\n",
    "The example provided here shows how to obtain the geometry optimization of the N$_2$ molecule by using the [Geopt](https://mmoriniere.gitlab.io/MyBigDFT/geopt.html) workflow class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialization\n",
    "\n",
    "You first need to import some useful classes to define a ground state calculation as well as the `Geopt` class, that is a workflow:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from mybigdft import Posinp, Atom, InputParams, Job\n",
    "from mybigdft.workflows import Geopt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The [Geopt](https://mmoriniere.gitlab.io/MyBigDFT/geopt.html) class\n",
    "\n",
    "This class allows to relax the input geometry of a given system in order to find the structure that minimizes the forces. It is meant to ease the creation of this type of calculation by automatically setting the main input parameters under the ``geopt`` key from the BigDFT input parameters.\n",
    "\n",
    "Of course, you can specify the input parameters you want to use for all these calculations, via the ``inputparams`` argument. You may also want to give a specific name to the runs and a specific folder where to run them, as usual.\n",
    "\n",
    "The base job has a given set of input parameters, and default values for the main parameters of the ``geopt`` key are automatically updated. The extra arguments of the ``geopt`` input parameters key can also be passed as keyword arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms = [Atom('N', [0, 0, 0]), Atom('N', [0, 0, 1.0975])]  # experimental geometry\n",
    "pos = Posinp(atoms, units=\"angstroem\", boundary_conditions=\"free\")\n",
    "inp = InputParams()\n",
    "base_job = Job(posinp=pos, inputparams=inp, name=\"N2\", run_dir=\"N2/geopt/from_experimental_structure\")\n",
    "geopt = Geopt(base_job)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The workflow is run as usual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/geopt/from_experimental_structure\n",
      "Logfile log-N2.yaml already exists!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "geopt.run(nmpi=6, nomp=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting from the experimental geometry, the converged atomic structure found using default BigDFT parameters is the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2   angstroem\n",
      "free\n",
      "N  -2.81097166026014e-21  -6.30580455586082e-22   0.00200502001148875\n",
      "N   2.83565311682474e-21   7.23934260935048e-22   1.09549494660991\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(geopt.final_posinp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.093489926598422\n"
     ]
    }
   ],
   "source": [
    "# Distance between both nitrogen atoms after geometry optimization\n",
    "distance = geopt.final_posinp.distance(0, 1)\n",
    "print(distance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Starting from another starting point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms_2 = [Atom('N', [0, 0, 0]), Atom('N', [0, 0, 1.05])]\n",
    "pos_2 = Posinp(atoms_2, units=\"angstroem\", boundary_conditions=\"free\")\n",
    "base_job_2 = Job(posinp=pos_2, inputparams=inp, name=\"N2\", run_dir=\"N2/geopt/from_worse_staring_point\")\n",
    "geopt_2 = Geopt(base_job_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can run this new workflow. Note that a warning was raised, telling us that, at some point, the input wavefunctions had to be recalculated: this is because the forces found for the initial geometry were high. The geometry optimization procedure then had to move the atoms by a large amount, which was sufficient for the input wavefunctions to be recalculated for the second geometry. For that reason, it is better to start from an initial position that is already close to a minimum. You might want to reach that \"rather good\" initial geometry with low quality input parameters (such as a small grid extension and a large grid step) before performing a more expensive one which gives converged results with respect to those input parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/geopt/from_worse_staring_point\n",
      "Logfile log-N2.yaml already exists!\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:265: UserWarning: The norm of the residue is too large, need to recalculate input wavefunctions\n",
      "  warnings.warn(warning, UserWarning)\n"
     ]
    }
   ],
   "source": [
    "geopt_2.run(nmpi=6, nomp=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same warning was actually found after the second step of the geometry optimization procedure. After that, the difference in initial geometry between two steps were low enough for the wavefunctions found at the end of the former step to be used as input guess for the next step, thus reducing the computational time. This is an indicator that your input geoemetry might not be close enough to the minimal geometry. The worst case scenario would be that the changes in the first steps of the geometry optimization procedure could leave you far away from the relaxed structure, and lead to a very long convergence (if any).\n",
    "\n",
    "You can see readily see that the second geopt workflow did need more steps than the first one to reach convergence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First geopt took 8 steps while second geopt took 27 steps.\n"
     ]
    }
   ],
   "source": [
    "n_1 = len(geopt.queue[0].logfile)\n",
    "n_2 = len(geopt_2.queue[0].logfile)\n",
    "print(f\"First geopt took {n_1} steps while second geopt took {n_2} steps.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It however reached a similar minimum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.093490323564331\n"
     ]
    }
   ],
   "source": [
    "distance = geopt_2.final_posinp.distance(0, 1)\n",
    "print(distance)"
   ]
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}