{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute the electronic polarizability tensor\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 electronic polarizability tensor of the N$_2$ molecule by using the [PolTensor](https://mmoriniere.gitlab.io/MyBigDFT/elecpoltensor.html#mybigdft.workflows.poltensor.PolTensor) class. \n",
    "\n",
    "The electronic polarizability tensor represents the response of the charges of a system (its dipole) to the application of an external electric field.\n",
    "To compute this tensor, some BigDFT calculations are performed, where the system is subject to an external electric field along each direction of space ($x$, $y$ and $z$). The elements of the electronic polarizability tensor are then defined by the ratio of the delta of the dipole in one direction and the delta of the electric field amplitudes:\n",
    "$\\alpha_{ij} = \\frac{\\Delta D_i}{\\Delta E_j}$\n",
    "where $i, j \\in \\{x, y, z\\}$, $\\Delta D_i$ is the variation of the dipole along the $i$ direction and $\\Delta E_j$ is the variation of the electric field amplitude along the $j$ direction"
   ]
  },
  {
   "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 `PolTensor` class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mybigdft import Job, Posinp, Atom\n",
    "from mybigdft.workflows import PolTensor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You are then able to define an instance of the `PolTensor` class. To that end, you only need to provide a Job defining the ground state calculation. It simply amounts to proving initial positions, some non-default input parameters (if needed) and initializing the job. Default values are used to define the electric field amplitudes applied along each direction, but you may specify them via the `ef_amplitudes` optional argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms = [Atom('N', [0, 0, 0]), Atom('N', [0, 0, 1.095])]\n",
    "pos = Posinp(atoms, units='angstroem', boundary_conditions='free')\n",
    "gs = Job(posinp=pos, name='N2', run_dir='../../../tests/pol_tensor_N2')\n",
    "pt = PolTensor(gs)\n",
    "# The line above actually amounts to:\n",
    "# pt = PolTensor(gs, ef_amplitudes=[1.e-4]*3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can then run the workflow without having to care about the context manager to run each calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/pol_tensor_N2\n",
      "Logfile log-N2.yaml already exists!\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/pol_tensor_N2/EF_along_x+\n",
      "Logfile log-N2.yaml already exists!\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/pol_tensor_N2/EF_along_y+\n",
      "Logfile log-N2.yaml already exists!\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/pol_tensor_N2/EF_along_z+\n",
      "Logfile log-N2.yaml already exists!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pt.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the calculations were run with different names in the same folder as the one specified by the ground state job: this was taken care of by the PolTensor class.\n",
    "\n",
    "You can loop over the jobs in the queue to see which were the input parameters used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N2: {}\n",
      "N2: {'dft': {'elecfield': [0.0001, 0.0, 0.0]}}\n",
      "N2: {'dft': {'elecfield': [0.0, 0.0001, 0.0]}}\n",
      "N2: {'dft': {'elecfield': [0.0, 0.0, 0.0001]}}\n"
     ]
    }
   ],
   "source": [
    "for job in pt.queue:\n",
    "    print(f\"{job.name}: {job.inputparams}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can finally get the value of the electronic polarizability tensor (in atomic units) via the `pol_tensor` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.05558e+01, -2.00000e-04,  0.00000e+00],\n",
       "       [-2.00000e-04,  1.05558e+01,  0.00000e+00],\n",
       "       [-2.00000e-04, -2.00000e-04,  1.50535e+01]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pt.pol_tensor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the electronic polarizability tensor is diagonal and $\\alpha_{xx} = \\alpha_{yy} < \\alpha_{zz}$. The values (in atomic units) are consistent with what is reported in the litterature and the experiment (see Table II of [George Maroulis and Ajit J. Thakkar, *J. Chem. Phys.* **88**, 7623 (1988)](http://dx.doi.org/10.1063/1.454327)). Keep in mind that default input parameters were used, so that the the present calculations cannot be considered as converged.\n",
    "\n",
    "**Be carefull**: if the electric fields applied are high enough for some electrons to escape the system, then you cannot trust the results anymore. In such a case, the system you are trying to simulate is not the one you expected, because an occupied orbital has a positive energy under the perturbation of the elctric field. There is actually a rule of thumb to know if the electric field applied along a direction is not too high: the difference of potential along the electric field (*i.e.* the amplitude of the electric field times the length of the BigDFT grid along the direction of the electric field) has to be roughly two times smaller than the absolute value of the HOMO when there is no perturbation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean polarizability and the polarizability tensor can easily be converted to $\\unicode[serif]{xc5}^3$. Another units that are often considered when it comes to polarizability are the SI units ($1 \\unicode[serif]{xc5}^3 = 1.1126501~10^{-40}$ C$^2$ m$^2$ J$^{-1}$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.7863720221055357"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from mybigdft.globals import B_TO_ANG\n",
    "pt.mean_polarizability*B_TO_ANG**3"
   ]
  }
 ],
 "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
}