10_DFPT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5bd44597-d660-405a-a4cf-ebae97b257aa",
   "metadata": {},
   "source": [
    "## Density Functional Perturbation Theory (DFPT)\n",
    "---\n",
    "\n",
    "In this notebook we will move to the perturbation theory calculation of phonon frequencies. This type of calculation is another step up on the scale of computational complexity and demands. Thus, the warnings about *not* following blindly the presented scheme in the real calculations applies here even more. The calculation is modeled after one of the AbInit tutorials.\n",
    "\n",
    "We will calculate for AsAl:\n",
    "- Dynamical matrix at $\\Gamma$ and X\n",
    "- Phonon frequencies at $\\Gamma$ and X\n",
    "- Born effective charges\n",
    "- Dielectric tensor\n",
    "- LO-TO splitting\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68bd7d6d-df37-49a8-b148-2a6df04e3f95",
   "metadata": {},
   "source": [
    "### Calculation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8a0b04f9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "%pylab is deprecated, use %matplotlib inline and import the required libraries.\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import ase\n",
    "import ase.io\n",
    "from ase.spacegroup import crystal\n",
    "from ase.calculators.abinit import Abinit\n",
    "from ase import units as un\n",
    "from elastic import get_pressure\n",
    "\n",
    "# Utility functions\n",
    "from tqdm.auto import tqdm\n",
    "import spglib\n",
    "from abilib import scan_param, set_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "814f5d1e-44df-4c70-88c9-0cec68c6f917",
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import os\n",
    "from ase.calculators.calculator import CalculationFailed\n",
    "\n",
    "def run_calc(cryst):\n",
    "    '''\n",
    "    Run the calculation without reading back any results into cryst.\n",
    "    '''\n",
    "    cryst.calc.write_inputfiles(cryst, ['energy', 'positions'])\n",
    "    cryst.calc.template.execute(cryst.calc.directory, cryst.calc.profile)\n",
    "    # cmd = cryst.calc.command.replace('PREFIX', cryst.calc.prefix)\n",
    "    # proc = subprocess.Popen(cmd, shell=True, cwd=cryst.calc.directory)\n",
    "    # errorcode = proc.wait()\n",
    "    # if errorcode:\n",
    "    #     path = os.path.abspath(cryst.calc.directory)\n",
    "    #     msg = ('Calculator \"{}\" failed with command \"{}\" failed in '\n",
    "    #            '{} with error code {}'.format(cryst.calc.name, cmd,\n",
    "    #                                           path, errorcode))\n",
    "    #     raise CalculationFailed(msg)\n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4493dee0-5fab-426b-a9ba-e61a244344f0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 3.9701,  3.9701,  3.9701, 60.    , 60.    , 60.    ])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = 5.61456927\n",
    "prim = crystal(('As','Al'), basis=((0,0,0), (1/4,1/4,1/4)), \n",
    "        spacegroup=216, cellpar=(a, a, a, 90, 90, 90),\n",
    "        primitive_cell=True)\n",
    "\n",
    "prim.cell.cellpar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f011cfad-83d1-4035-8753-5d45fc78c1d7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.  , 0.  , 0.  ],\n",
       "       [0.25, 0.25, 0.25]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prim.get_scaled_positions()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ce34ce99-b7d9-49e8-948f-8cd09bb0aea4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'F-43m (216)'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spglib.get_spacegroup((prim.cell, prim.get_scaled_positions(), prim.get_atomic_numbers()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "92e92ec1-c2e5-4892-97f8-0a0fc303df60",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FCC(a=5.61457)\n",
      "  Variant name: FCC\n",
      "  Special point names: GKLUWX\n",
      "  Default path: GXWKGLUWLK,UX\n",
      "\n",
      "  Special point coordinates:\n",
      "    G   0.0000  0.0000  0.0000\n",
      "    K   0.3750  0.3750  0.7500\n",
      "    L   0.5000  0.5000  0.5000\n",
      "    U   0.6250  0.2500  0.6250\n",
      "    W   0.5000  0.2500  0.7500\n",
      "    X   0.5000  0.0000  0.5000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "lat = prim.cell.get_bravais_lattice()\n",
    "print(lat.description())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "79484e15-4c80-40f9-a0b7-a11bd29de8fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "work_dir = 'work_10'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0f6020a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_calculator(directory=work_dir, ecut=420):\n",
    "    return Abinit(directory=directory,             \n",
    "              ecut= ecut,\n",
    "              #pps = 'pspnc', tolvrs=1.0e-18, nband=4,\n",
    "              pps = 'pawxml', pawecutdg = 750, xc = 'LDA', toldfe = 1.0e-6,\n",
    "              diemac=9.0,\n",
    "              nstep=15,\n",
    "              kptrlatt=[[-4,  4,  4],\n",
    "                        [ 4, -4,  4],\n",
    "                        [ 4,  4, -4]]\n",
    "              )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ae9d4731",
   "metadata": {},
   "outputs": [],
   "source": [
    "prim.calc = create_calculator()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "56ee0fcb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Energy:  -3164.646 eV\n",
      "  Stress: [   -0.44  -0.44  -0.44   0.00   0.00   0.00 ] GPa\n",
      "Pressure:     0.444 GPa\n"
     ]
    }
   ],
   "source": [
    "print(f'  Energy:  {prim.get_potential_energy():.3f} eV')\n",
    "print( '  Stress: [', 6*' %6.2f' % tuple(prim.get_stress()/un.GPa), '] GPa')\n",
    "print(f'Pressure:     {get_pressure(prim.get_stress())/un.GPa:.3f} GPa')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a6020342",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0., -0.,  0.],\n",
       "       [ 0., -0., -0.]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prim.get_forces()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7f4025d-04a8-4c98-b348-30e38c3abada",
   "metadata": {},
   "source": [
    "### Convergence scan\n",
    "---\n",
    "This is a reminder that you should *always* perform convergence analysis of your calculations. We can skip this step by setting `Scan_Corvengence = False` , since we have already performed this analysis in previous part."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "75bb908b",
   "metadata": {},
   "outputs": [],
   "source": [
    "Scan_Corvengence = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c976ad45",
   "metadata": {},
   "outputs": [],
   "source": [
    "if Scan_Corvengence :\n",
    "    ecut = prim.calc.parameters['ecut']\n",
    "    ecuts, eres = scan_param(prim, 'ecut', 200, 500)\n",
    "\n",
    "    figsize(6,6)\n",
    "    subplot(211)\n",
    "    plot(ecuts, eres[0], 'o')\n",
    "    ylabel('Energy (eV)')\n",
    "    subplot(212)\n",
    "    plot(ecuts, eres[1]/un.GPa, 'o')\n",
    "    ylabel('Pressure (GPa)')\n",
    "    xlabel('ecut (eV)');\n",
    "    prim.calc.set(ecut=ecut)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f516d8b7-af7b-4c9e-98d5-b0650e29423c",
   "metadata": {},
   "source": [
    "### DFPT calculations\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "151b42f9-164b-4f00-852e-10f01eaa1567",
   "metadata": {},
   "outputs": [],
   "source": [
    "cr = ase.Atoms(prim)\n",
    "cr.calc = set_params(create_calculator(f'{work_dir}/AsAl_3'),\n",
    "            ndtset=4,\n",
    "\n",
    "            #Ground state calculation\n",
    "            kptopt1=1,              # Automatic generation of k points, taking\n",
    "                                    # into account the symmetry\n",
    "            tolvrs1=1.0e-18,        # SCF stopping criterion\n",
    "            \n",
    "            #Response Function calculation : d/dk\n",
    "            rfelfd2=2,              # Activate the calculation of the d/dk perturbation\n",
    "            rfdir2=[1, 0, 0],       # Need to consider the perturbation in the x-direction only\n",
    "                                    # This is rather specific, due to the high symmetry of the AlAs crystal\n",
    "                                    # In general, just use rfdir 1 1 1\n",
    "                                    # In the present version of ABINIT (v4.6), symmetry cannot be used\n",
    "                                    # to reduce the number of ddk perturbations\n",
    "            nqpt2=1,\n",
    "            qpt2=[0.0, 0.0, 0.0],   # This is a calculation at the Gamma point\n",
    "\n",
    "            getwfk2=-1,             # Uses as input the output wf of the previous dataset\n",
    "\n",
    "            kptopt2=2,              # Automatic generation of k points,\n",
    "                                    # using only the time-reversal symmetry to decrease\n",
    "                                    # the size of the k point set.\n",
    "\n",
    "            iscf2=-3,               # The d/dk perturbation must be treated \n",
    "                                    # in a non-self-consistent way\n",
    "            tolwfr2=1.0e-22,        # Must use tolwfr for non-self-consistent calculations\n",
    "                                    # Here, the value of tolwfr is very low.          \n",
    "            \n",
    "            #Response Function calculation : electric field perturbation and phonons\n",
    "            rfphon3=1,              # Activate the calculation of the atomic dispacement perturbations\n",
    "            rfatpol3=[1, 2],        # All the atoms will be displaced\n",
    "            rfelfd3=3,              # Activate the calculation of the electric field perturbation\n",
    "            rfdir3=[1, 1, 1],       # All directions are selected. However, symmetries will be used to decrease\n",
    "                                    # the number of perturbations, so only the x electric field is needed\n",
    "                                    # (and this explains why in the second dataset, rfdir was set to 1 0 0).\n",
    "\n",
    "            nqpt3=1,\n",
    "            qpt3=[0.0, 0.0, 0.0],   # This is a calculation at the Gamma point\n",
    "\n",
    "            getwfk3=-2,             # Uses as input wfs the output wfs of the dataset 1\n",
    "            getddk3=-1,             # Uses as input ddk wfs the output of the dataset 2\n",
    "\n",
    "            kptopt3=2,              # Automatic generation of k points,\n",
    "                                    # using only the time-reversal symmetry to decrease\n",
    "                                    # the size of the k point set.\n",
    "            tolvrs3=1.0e-8,\n",
    "            \n",
    "            #Response Function calculation : phonon at X\n",
    "              rfphon4 = 1,            # Activate the calculation of the atomic dispacement perturbations\n",
    "             rfatpol4 = [1,2],           # All the atoms will be displaced\n",
    "               rfdir4 = [1,1,1],         # Need all directions but ABINIT will use symmetries if possible\n",
    "\n",
    "                nqpt4 = 1,\n",
    "                 qpt4 = [0.0,0.5,0.5],   # This is a calculation at the X point\n",
    "\n",
    "              getwfk4 = 1,            # Uses as input wfs the output wfs of the dataset 1\n",
    "\n",
    "              kptopt4 = 3,            # Automatic generation of k points,\n",
    "                                      # no use of symmetries to decrease \n",
    "                                      # the size of the k point set.\n",
    "              tolvrs4 = 1.0e-8,\n",
    "\n",
    "           )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "81dca93c-9c45-4f9f-b876-d95cfee9c232",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_calc(cr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dacb6349-e8d6-40fd-89e6-b2033a7b0081",
   "metadata": {},
   "source": [
    "### Extract the results\n",
    "---\n",
    "\n",
    "The following code is a simple state machine extracting data from the output file.\n",
    "The code is an ad-hoc type but can be adapted to extracting other data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "610609bd-02e0-4784-bfef-60adbbf3e819",
   "metadata": {},
   "outputs": [],
   "source": [
    "cr.calc.directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "9f823bf3-b8ac-4594-870c-075a111778ce",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Dielectric tensor:\n",
      " [[ 9.75441052  0.          0.        ]\n",
      " [-0.          9.75441052 -0.        ]\n",
      " [-0.         -0.          9.75441052]]\n",
      "\n",
      "Effective charges (from electric field response):\n",
      "[[13.27726459  0.         -0.        ]\n",
      " [ 2.2279416  -0.          0.        ]\n",
      " [ 7.75708322  5.52018137 -7.75708322]\n",
      " [ 0.05350041  2.1744412  -0.05350041]\n",
      " [ 7.75708322 -7.75708322  5.52018137]\n",
      " [ 0.05350041 -0.05350041  2.1744412 ]]\n",
      "\n",
      "Effective charges (from phonon response):\n",
      "[[13.27720027  7.75708322  7.75708322]\n",
      " [ 0.          5.52011705 -7.75708322]\n",
      " [-0.         -7.75708322  5.52011705]\n",
      " [ 2.22800672  0.05350041  0.05350041]\n",
      " [-0.          2.17450631 -0.05350041]\n",
      " [ 0.         -0.05350041  2.17450631]]\n",
      "\n",
      "k-point: [0.0, 0.0, 0.0]\n",
      "Frequencies (cm^-1):\n",
      " [31.54545, 31.54545, 31.54545, 361.4994, 361.4994, 361.4994]\n",
      "\n",
      "Direction: [1.0, 0.0, 0.0]\n",
      "Frequencies (cm^-1):\n",
      " [31.54545, 31.54545, 339.4295, 361.4994, 361.4994, 552.8918]\n",
      "\n",
      "Direction: [0.0, 1.0, 0.0]\n",
      "Frequencies (cm^-1):\n",
      " [31.54545, 31.54545, 303.1788, 361.4994, 361.4994, 537.8272]\n",
      "\n",
      "Direction: [0.0, 0.0, 1.0]\n",
      "Frequencies (cm^-1):\n",
      " [31.54545, 31.54545, 303.1788, 361.4994, 361.4994, 537.8272]\n",
      "\n",
      "\n",
      "k-point: [0.0, 0.5, 0.5]\n",
      "Frequencies (cm^-1):\n",
      " [96.79644, 96.79644, 216.2615, 336.1106, 336.1106, 391.4759]\n"
     ]
    }
   ],
   "source": [
    "with open(f'{cr.calc.directory}/abinit.abo') as of:\n",
    "    st = 0\n",
    "    for l in of:\n",
    "        if st==0 and 'Phonon wavevector' in l :\n",
    "            st=1\n",
    "            print(f'\\nk-point: {[float(v) for v in l.split()[-3:]]}')\n",
    "\n",
    "        if 0<st:\n",
    "            #print('>>>', l.strip())\n",
    "            pass\n",
    "            \n",
    "        tok = [] \n",
    "        if st > 0:\n",
    "            tok = l[1:].split()\n",
    "\n",
    "        if st==2 and len(tok)==6:\n",
    "            dt.append(float(l.split()[4]))\n",
    "            if len(dt)==9:\n",
    "                dt=array(dt).reshape((3,3))\n",
    "                print('\\nDielectric tensor:\\n', dt)\n",
    "                st=0\n",
    "            \n",
    "        if st==0 and 'Dielectric tensor' in l:\n",
    "            dt=[]\n",
    "            st=2\n",
    "            continue\n",
    "\n",
    "            \n",
    "        if st==31 and len(tok)==6:\n",
    "            ec.append(float(l.split()[4]))\n",
    "            if len(ec)==18:\n",
    "                ec=array(ec).reshape((-1,3))\n",
    "                print(ec)\n",
    "                st=0                  \n",
    "\n",
    "        if st==3 and '(from':\n",
    "            print(f'\\nEffective charges {l.strip()}:')\n",
    "            st=31\n",
    "            continue\n",
    "            \n",
    "        if st==0 and 'Effective charges' in l:\n",
    "            st=3\n",
    "            ec=[]\n",
    "            continue\n",
    "            \n",
    "            \n",
    "        if st==11 :\n",
    "            if tok:\n",
    "                frqs += [float(v) for v in tok]\n",
    "            else :\n",
    "                print('Frequencies (cm^-1):\\n', frqs)\n",
    "                st=1\n",
    "            \n",
    "        if st==1 :\n",
    "            if 'cm-1' in l:\n",
    "                st=11\n",
    "                frqs = []\n",
    "            if 'direction' in tok[:1]:\n",
    "                print(f'\\nDirection: {[float(v) for v in l.split()[-3:]]}')\n",
    "        \n",
    "        if st==1 and '== DATASET' in l:\n",
    "            print()\n",
    "            st=0\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ea94c08-6900-482b-898c-ce1027140a31",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "878303e5-99c3-4e21-b783-35e175a0865b",
   "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.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}