{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 3: Density functional theory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Authors:  <u>Jan Janssen</u>, Tilmann Hickel, Jörg Neugebauer ([Max-Planck-Institut für Eisenforschung](https://www.mpie.de))**\n",
    "\n",
    "The scope of this first exercise is to become familar with:\n",
    "* Density functional theory calculation, \n",
    "* Calculate the equilibrium lattice constant \n",
    "* Compare the results to interatomic potentials."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='red'> Update Notice </font> \n",
    "Since the workshop in 2020, there have been some updates to pyiron commands, and here in the exercises new commands are used. \n",
    "\n",
    "The changes include:\n",
    "- Now pyiron has various modules, like pyiron_atomistics, pyiron_continuum, pyiron_experimental. So we can import `Project` from `pyiron_atomistics` directly.\n",
    "- using `pr.create` to create instances of various objects, such as structures and jobs. For example:\n",
    "```python\n",
    "pr.create.structure.ase.bulk() ## this replaces pr.create_ase_bulk()\n",
    "```\n",
    ", or\n",
    "```python\n",
    "pr.create.job.Lammps('job_name') ## this replaces pr.pr.create_job(job_type=pr.job_type.Lammps,'job_name')\n",
    "```\n",
    "\n",
    "Similary, one can create pyiron tables via `pr.create.table()`.\n",
    "\n",
    "By this approach, the user easily gets the available options via autocompeletion.\n",
    "\n",
    "**Please note that in the video tutorials, the old commands are still presented.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reminder\n",
    "In the first session we learned how to create a pyiron project object and then use this pyiron project object to create atomistic structure objects. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the Project object\n",
    "from pyiron_atomistics import ____"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a Project object instance for a project named dft\n",
    "pr = ____(\"dft\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a cubic aluminum fcc structure\n",
    "al_fcc = pr.create.structure.ase.bulk(____, _____=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Confirm the final structure has 4 atoms by calculating \n",
    "# the length of the structure object\n",
    "____(al_fcc) == 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Density Functional Theory \n",
    "To provide a brief introduction to density functional theory (DFT) we calculate the equilibrium lattice constant for a cubic fcc aluminium structure. Besides the pseudo potential which defines the electron electron interaction, the DFT precision is dominated by the convergence parameters, namely the plane wave energy cutoff and the kpoint mesh. Both can be set in pyiron using the corresponding properties `encut` and `kpoint_mesh`.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a DFT job with the S/PHI/nX quantum engine named spx\n",
    "job_dft_template = pr.create.job.Sphinx(\n",
    "    job_name=________\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the default DFT convergence parameters energy cutoff\n",
    "print(\n",
    "    job_dft_template.encut, \n",
    "    job_dft_template.kpoint_mesh\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Increase the energy cutoff to 400eV and increase the kpoint mesh to 5x5x5\n",
    "job_dft_template.set_encut(_______)\n",
    "job_dft_template.set_kpoints([__, __, __])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Energy Volume Curve \n",
    "The discretisation of the plane waves on the fourier mesh and the kpoint mesh in the brillouin zone cause fluctuations in the energy surface. Therefore a minimization like we used it for the interatomic potentials in the previous section is insufficient. Instead we calculate the energy for various volumes around the equlibrium volume and then determine the equilibrium volume by interpolating the minimum between these volumes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the numerical library numpy\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create 7 strains ranging from -5% (-0.05) to +5% (0.05)\n",
    "strain_lst = np.linspace(_____, ______, _____)\n",
    "strain_lst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copy the cubic fcc aluminium supercell\n",
    "al_fcc_copy = _________.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply of -5% to the copy of the aluminium\n",
    "# supercell and compare the volume \n",
    "al_fcc_copy.apply_strain(_____)\n",
    "al_fcc.get_volume(), al_fcc_copy.get_volume()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Iterate over the list of strains \n",
    "for _______ in strain_lst: \n",
    "    job_strain = job_dft_template.copy_to(\n",
    "        # Define the job name based on the current strain \n",
    "        new_job_name=\"spx_\" + str(1 + _____).replace(\".\", \"_\"), \n",
    "        new_database_entry=False\n",
    "    )\n",
    "    \n",
    "    # Copy the cubic fcc aluminium supercell\n",
    "    al_fcc_copy = ________.copy()\n",
    "    \n",
    "    # Apply the strain of using the for loop\n",
    "    al_fcc_copy.apply_strain(_____)\n",
    "    \n",
    "    # Set the strained structure to the job\n",
    "    job_strain.structure = al_fcc_copy \n",
    "    \n",
    "    # Execute the job\n",
    "    job_strain.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the status of the calculation in the job_table\n",
    "pr.job_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reminder\n",
    "We again use the pyiron table object to collect the simulation results and aggregate them in a pandas DataFrame. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a pyiron table object\n",
    "table = pr.create.table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Implement a filter function, which returns true \n",
    "# for finished jobs and jobs with \"spx_\" in the job_name\n",
    "def filter_jobs(job):\n",
    "    return job.status == _______ and _____ in job.job_name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Many commonly used functions are already available using tab completion\n",
    "# We select the get_volume function and the get_energy_tot function\n",
    "table.add.________\n",
    "table.add._________"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assign the filter function defined above\n",
    "table.filter_function = _____"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute the pyiron table just like a pyiron job object\n",
    "table._____"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Return a pandas DataFrame with the collected results \n",
    "df_res = table.get_dataframe()\n",
    "df_res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualise the energy volume curve \n",
    "We again use the matplotlib library to visualise the calculated energy volume curve and calculate the equilibrium volume by fitting a second order polynomial and calculate the roots of the derivative. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit a second order polynomial to the energy volume curve \n",
    "fit = np.polyfit(_____.volume, _____.energy_tot, ___)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the polynomial of the fit \n",
    "fit_poly = np.poly1d(fit)\n",
    "\n",
    "# Calculate the roots of the derivate of the polynomial\n",
    "vol_roots = np.polyder(fit_poly).roots\n",
    "\n",
    "# Select the root within the volume range which is\n",
    "# smaller (<) than maximum volume and larger (>) than \n",
    "# the minimum volume. \n",
    "vol_eq = vol_roots[\n",
    "    (vol_roots ____ df_res[\"volume\"].max()) & \n",
    "    (vol_roots ____ df_res[\"volume\"].min())\n",
    "][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the matplotlib library for plotting. \n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Plot the volume and the total energy from the DataFrame \n",
    "# with the collected results\n",
    "plt.plot(_____.volume, _____.energy_tot)\n",
    "\n",
    "# Plot the fitted equilibrium volume\n",
    "plt.axvline(______, linestyle=\"--\", color=\"red\")\n",
    "plt.xlabel(\"Volume\")\n",
    "plt.ylabel(\"total Energy\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## pyiron Master Jobs \n",
    "While managing calculations with with for loops and aggregating calculation results in pandas DataFrame is already a very scaleable concept, we have implemented master jobs which can execute multiple calculations in series or in parallel to automate common simulation tasks. The calculation of the energy volume curve is one of those examples. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a DFT job with the S/PHI/nX quantum engine named \"spxjob\"\n",
    "job_master_template = pr.create.job.Sphinx(\n",
    "    job_name=______\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assign the cubic fcc aluminium structure to the template job \n",
    "job_master_template.structure = _____"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the job object to create the Murnaghan object, named \"murn\" \n",
    "murn = pr.create.job.Murnaghan( \n",
    "    job_name=_______\n",
    ")\n",
    "murn.ref = job_master_template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute the Murnaghan object \n",
    "murn.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the status of the calculation in the job_table\n",
    "pr.job_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Finally we can use the internal functionality of the master job \n",
    "# to visualise the energy volume curve\n",
    "murn.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary \n",
    "In this section you learned: \n",
    "* how to calculate equilibrium bulk material properties using density functional theory.\n",
    "* how the finite plane wave energy cutoff and the finite kpoint mesh limit the precision of DFT calculation. \n",
    "* and how to use master jobs in pyiron to automate common tasks like calculating energy volume curves. \n",
    "\n",
    "Suggestions: \n",
    "* To learn more about convergence you can plot the convergence of an individual energy over the change in plane wave energy cutoff or kpoint mesh. \n",
    "* To visualise the discretisation of the energy volume curve you can calculate the energy volume curve at a small volume range of +/- 1% for a low kpoint mesh of 3x3x3 and an energy cutoff of 300eV with 21 points.\n",
    "* To validate the interatomic potentials from the previous section we can calculate energy volume curves for those and compare the energy differences to the DFT results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}