IntroductionPotentialFitting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f2adfe85-18c8-4430-b370-1b4d6f1b05b3",
   "metadata": {},
   "source": [
    "# Interatomic potentials\n",
    "\n",
    "In general interatomic potentials can be written as a sum of functional terms depending on the positions of the atoms in a structure. Then the energy $U$ of the system is\n",
    "\n",
    "$U = \\sum_{i=1}^N U_1(\\vec{r_i}) + \\sum_{i,j=1}^N U_2(\\vec{r_i}, \\vec{r_j}) + \\sum_{i,j,k=1}^N U_2(\\vec{r_i}, \\vec{r_j}, \\vec{r_k}) +...$\n",
    "\n",
    "The one body term only matter when the system is within an external field. Most classic interatomic potentials don't use 4th and higher order body terms, pair potentials only use the 2 body term. As a general rule potentials that include higher order body terms can be more accurate but are slower.\n",
    "\n",
    "There are many different forms for $U_i$, which cover a wide range of different run time and accuracy necessities.\n",
    "\n",
    "Simple pair potentials (f.e Lennard-Jones, Morse, Buckingham) only contain very few parameters, typically less than 10. In many body potentials (f.e. EAM, MEAM, Tersoff) the typical number of parameters is 10-50 and for machine learning potentials the number of parameters can reach several thousands.\n",
    "\n",
    "# Fitting\n",
    "\n",
    "In the fit process the parameters of the chosen functions for $U_i$ are optimized. For this purpose an objective or cost function is defined and minimized. In general the objective function is defined as\n",
    "\n",
    "$\\chi^2 = \\sum_i w_i r_i$\n",
    "\n",
    "where $w_i$ is a weight and $r_i$ is a residual that describes the difference to target values. This residual can be defined in different ways, so it is not possible to simply compare the residual for different fitting processes or codes. A more in depth explanation and some examples can be found on https://atomicrex.org/overview.html#objective-function.\n",
    "\n",
    "The minimization can be done with local or global optimization algorithms.\n",
    "Generally local optimization algorithms should all be able to find the local minimum coming from some initial parameter set, so the \"best\" local algorithm is the one finding the minimum in the shortest time. Typically used local algorithms are f.e. (L)BFGS or Nelder-Mead.\n",
    "Examples for global algorithms are evolutionary algorithms or simulated annealing. For most problems it is impossible to tell a priori which global algorithm will give the best results, so using global algorithms typically involves testing many of them.\n",
    "\n",
    "# EAM potentials\n",
    "\n",
    "EAM potentials are pair functionals. \n",
    "In a generalised form they are equal to Finnis-Sinclair, effective medium theory or glue potentials. Their total energy can be written as\n",
    "\n",
    "$E = \\frac{1}{2}\\sum_{ij}V(r_{ij}) + \\sum_i F(\\rho_i)$\n",
    "\n",
    "with\n",
    "\n",
    "$\\rho_i = \\sum_j \\rho(r_{ij})$\n",
    "\n",
    "The original functions for V, $\\rho$ and F were derived from different theories, but they can be chosen freely.\n",
    "\n",
    "# Fitting code\n",
    "\n",
    "Fitting is done using the pyiron interface to the atomicrex code https://atomicrex.org. It can be used to fit different types of classic interatomic potentials:\n",
    "- pair potentials\n",
    "- EAM\n",
    "- MEAM\n",
    "- Tersoff\n",
    "- ABOP\n",
    "- ADP (in development)\n",
    "\n",
    "It allows to fit different properties (energies, forces, lattice constants, elastic properties, etc.) and implements the LBFGS minimizer. Additionally it offers an interface to the nlopt library which implements several global and local optimization algorithms and the ability to apply arbitrary constraints to parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f4075857-362a-400f-b043-2fc5ad54a407",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyiron import Project, ase_to_pyiron\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d70cf5ae-8e25-444c-b24f-2d676ef8d2f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "pr = Project(\".\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d83afb5-b99d-4081-9883-ce67c3557f5d",
   "metadata": {},
   "source": [
    "### Get the training data\n",
    "Load a job that contains structures with energies and forces from DFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ea3fcb75-1d69-4f9e-901d-ba001c6dab2e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The job TrainData was saved and received the ID: 738901\n"
     ]
    }
   ],
   "source": [
    "tc = pr.create.job.TrainingContainer(\"TrainData\", delete_existing_job=True)\n",
    "df = pd.read_pickle(\"Al_3PreFinal.pckl.gz\")\n",
    "for r in df.itertuples():\n",
    "    s = ase_to_pyiron(r.ase_atoms)\n",
    "    if r.energy_corrected > 0:\n",
    "        continue\n",
    "    vol = s.get_volume()\n",
    "    natoms = len(s)\n",
    "    if vol/natoms > 120:\n",
    "        continue\n",
    "    tc.add_structure(structure=s, energy=r.energy_corrected, forces=r.forces, name=f\"ID{r.Index}\")\n",
    "tc.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57cc1f71-1255-4d15-9c31-e598098fe1e4",
   "metadata": {},
   "source": [
    "### Have a look at the training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "25b99c5a-e8a2-42de-aaf0-2794d675385b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>V</th>\n",
       "      <th>E</th>\n",
       "      <th>space_group</th>\n",
       "      <th>crystal_system</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>16.484415</td>\n",
       "      <td>-3.482751</td>\n",
       "      <td>225</td>\n",
       "      <td>cubic</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>14.835973</td>\n",
       "      <td>-3.433909</td>\n",
       "      <td>225</td>\n",
       "      <td>cubic</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>15.165661</td>\n",
       "      <td>-3.452735</td>\n",
       "      <td>225</td>\n",
       "      <td>cubic</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>15.495350</td>\n",
       "      <td>-3.466526</td>\n",
       "      <td>225</td>\n",
       "      <td>cubic</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>15.825038</td>\n",
       "      <td>-3.475815</td>\n",
       "      <td>225</td>\n",
       "      <td>cubic</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",