{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "matplotlib.use(\"Pdf\")\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'svg'\n",
    "\n",
    "# set matplotlib parameters\n",
    "params = {'backend': 'ps',\n",
    "          'axes.labelsize': 13,\n",
    "          'font.size': 13,\n",
    "          'legend.fontsize': 13,\n",
    "          'xtick.labelsize': 13,\n",
    "          'ytick.labelsize': 13,\n",
    "          'text.usetex': False,\n",
    "          'text.latex.preamble': \"\\\\usepackage{mathpazo}, \\\\usepackage{amsmath}\",\n",
    "          'font.family': 'sans-serif',\n",
    "          'font.sans-serif': ['Computer Modern Sans serif']}\n",
    "plt.rcParams.update(params)\n",
    "\n",
    "import warnings \n",
    "warnings.filterwarnings(\"ignore\") #ignore some matplotlib warnings\n",
    "\n",
    "# numpy\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this basic example we will perform a VASP calculation for SrVO$_3$, build PLOs for the Vanadium t$_{2g}$ orbitals, and load them as SumK object, which can be then used to perform a DMFT calculation.\n",
    "\n",
    "## VASP setup\n",
    "\n",
    "First we setup the VASP [INCAR link](INCAR) file by specifing the LOCPROJ, EMIN, EMAX and LORBIT flags:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "SYSTEM = SrVO3\n",
    "NCORE = 4\n",
    "LMAXMIX=6\n",
    "EDIFF = 1.E-10\n",
    "\n",
    "# DOS energy window\n",
    "NEDOS = 2001\n",
    "\n",
    "! switch off symmetries\n",
    "ISYM=-1\n",
    "\n",
    "# Smearing procedure\n",
    "ISMEAR = -5\n",
    "\n",
    "# the energy window to optimize projector channels\n",
    "EMIN = 3.9\n",
    "EMAX = 7.1\n",
    "\n",
    "# use the PAW channel optimization\n",
    "LORBIT=14\n",
    "\n",
    "# project to V d\n",
    "LOCPROJ = 2 : d : Pr\n",
    "```\n",
    "Moreover we prepare a [KPOINTS link](KPOINTS), [POSCAR link](POSCAR), and a POTCAR file. For the POTCAR file please use the VASP provided PBE pseudopotentials: `Sr_sv`, `V`, and `O`. \n",
    "\n",
    "Now VASP is executed, which should converge in roughly 27 iterations. Afterwards you should find the files LOCPROJ and PROJCAR in you directory. \n",
    "\n",
    "## PLOVASP\n",
    "\n",
    "First import the PLOVASP module of DFTTools:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Starting run with 1 MPI rank(s) at : 2019-12-05 16:12:52.689539\n"
     ]
    }
   ],
   "source": [
    "# import plovasp converter\n",
    "import triqs_dft_tools.converters.plovasp.converter as plo_converter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, create a configuration file for plovasp [plo.cfg link](plo.cfg):\n",
    "\n",
    "```\n",
    "[General]\n",
    "DOSMESH = -3.0 3.0 2001\n",
    "\n",
    "[Shell 1]\n",
    "LSHELL = 2\n",
    "IONS = 2\n",
    "EWINDOW = -1.4 2.0\n",
    "\n",
    "TRANSFORM = 1.0  0.0  0.0  0.0  0.0\n",
    "            0.0  1.0  0.0  0.0  0.0\n",
    "            0.0  0.0  0.0  1.0  0.0\n",
    "\n",
    "```\n",
    "where the energy window of the t$_{2g}$ bands is specified by `EWINDOW` and the `TRANSFORM` flag picks the correct three orbitals out of the five Vanadium $d$ orbitals [see the guide for the ordering of orbitals](../../guide/conv_vasp.html). Before running PLOVASP, make sure that the Fermi energy is written in the first line of the `LOCPROJ` file, or copy it there (see the VASP [interface guide](../../guide/conv_vasp.html) for more information). The first line should look like\n",
    "```\n",
    "1   729    21     5   5.3834262  # of spin, # of k-points, # of bands, # of proj, Efermi\n",
    "```\n",
    "\n",
    "Now run PLOVASP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Read parameters:\n",
      "0  ->  {'m': 0, 'l': 2, 'isite': 2, 'label': 'dxy'}\n",
      "1  ->  {'m': 1, 'l': 2, 'isite': 2, 'label': 'dyz'}\n",
      "2  ->  {'m': 2, 'l': 2, 'isite': 2, 'label': 'dz2'}\n",
      "3  ->  {'m': 3, 'l': 2, 'isite': 2, 'label': 'dxz'}\n",
      "4  ->  {'m': 4, 'l': 2, 'isite': 2, 'label': 'dx2-y2'}\n",
      "  Found POSCAR, title line: SrVO3\n",
      "  Total number of ions: 5\n",
      "  Number of types: 3\n",
      "  Number of ions for each type: [1, 1, 3]\n",
      "\n",
      "    Total number of k-points: 729\n",
      "   Total number of tetrahedra: 4374\n",
      "eigvals from LOCPROJ\n",
      "\n",
      "  Unorthonormalized density matrices and overlaps:\n",
      "  Spin: 1\n",
      "  Site: 2\n",
      "  Density matrix                                                  Overlap\n",
      "   0.5875772   0.0015679  -0.0003707   0.0015674   0.0000000       0.9294791  -0.0000080  -0.0000078  -0.0000080  -0.0000001\n",
      "   0.0015679   0.5876177  -0.0001854  -0.0016078   0.0003240      -0.0000080   0.9294790  -0.0000042   0.0000080   0.0000070\n",
      "  -0.0003707  -0.0001854   0.5815486  -0.0001854  -0.0000000      -0.0000078  -0.0000042   0.9715751  -0.0000038   0.0000003\n",
      "   0.0015674  -0.0016078  -0.0001854   0.5876172  -0.0003240      -0.0000080   0.0000080  -0.0000038   0.9294791  -0.0000066\n",
      "   0.0000000   0.0003240  -0.0000000  -0.0003240   0.5815487      -0.0000001   0.0000070   0.0000003  -0.0000066   0.9715748\n",
      "\n",
      "  Generating 1 shell...\n",
      "\n",
      "    Shell         : 1\n",
      "    Orbital l     : 2\n",
      "    Number of ions: 1\n",
      "    Dimension     : 3\n",
      "    Correlated    : True\n",
      "    Ion sort      : [1]\n",
      "Density matrix:\n",
      "  Shell 1\n",
      "    Site 1\n",
      "     0.3332630     0.0021719     0.0021714\n",
      "     0.0021719     0.3333128    -0.0022211\n",
      "     0.0021714    -0.0022211     0.3333123\n",
      "      trace:  0.9998880790966638\n",
      "\n",
      "  Impurity density: 0.9998880790966638\n",
      "\n",
      "Overlap:\n",
      "  Site 1\n",
      "[[ 1. -0.  0.]\n",
      " [-0.  1.  0.]\n",
      " [ 0.  0.  1.]]\n",
      "\n",
      "Local Hamiltonian:\n",
      "  Shell 1\n",
      "    Site 1\n",
      "     0.5633806     0.0007563     0.0007563\n",
      "     0.0007563     0.5633801    -0.0007559\n",
      "     0.0007563    -0.0007559     0.5633801\n",
      "\n",
      "Evaluating DOS...\n",
      "  Shell 1\n",
      "    Total number of states: [[[7.33737319 7.48285647 7.28002405]]]\n",
      "  Storing ctrl-file...\n",
      "  Storing PLO-group file 'vasp.pg1'...\n",
      "  Density within window: 0.9999741659673522\n"
     ]
    }
   ],
   "source": [
    "# Generate and store PLOs\n",
    "plo_converter.generate_and_output_as_text('plo.cfg', vasp_dir='./')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PLOVASP created one shell with three orbitals, which are equally filled by 1/3, one electron in total. Additionally we calculated the density of states. Both in VASP and PLOVASP. The later stores the data in the file pdos_x.dat, which can be simply plotted with [matplotlib](https://matplotlib.org/). The result should look similar to:\n",
    "\n",
    "![](svo_dos.svg)\n",
    "\n",
    "Here the gray area highlights the energy window for the PLOs. The total DOS of VASP (blue) coincides with the PLO DOS in the window, as we re-orthonormalized the projector functions in the given window, picking up also Oxygen weight. This setting is closed to the result of maximally localized Wannier functions created with [wannier90](http://www.wannier.org/) without running the actual minimization of the spread. Note, for a proper comparison one can use the hydrogen projector in VASP by using the the line `LOCPROJ=  2 : d : Hy`, instead of `Pr`. \n",
    "\n",
    "\n",
    "## Converting to hdf5 file\n",
    "\n",
    "Finally we can run the VASP converter to create a h5 file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading input from vasp.ctrl...\n",
      "{\n",
      "    \"ngroups\": 1,\n",
      "    \"nk\": 729,\n",
      "    \"ns\": 1,\n",
      "    \"nc_flag\": 0\n",
      "}\n",
      "\n",
      "  No. of inequivalent shells: 1\n"
     ]
    }
   ],
   "source": [
    "# import VASPconverter\n",
    "from triqs_dft_tools.converters.vasp import *\n",
    "\n",
    "\n",
    "# create Converter\n",
    "Converter = VaspConverter(filename = 'vasp')\n",
    "# run the converter\n",
    "Converter.convert_dft_input()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting h5 file `vasp.h5` can now be loaded as sumk object via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SumK\n",
    "from triqs_dft_tools.sumk_dft_tools import SumkDFTTools\n",
    "\n",
    "SK = SumkDFTTools(hdf_file='vasp.h5', use_dft_blocks = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here one should carefully determine the block structure manually. This is important to find degenerate orbitals and spin-channels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "found 1 blocks of degenerate orbitals in shell 0\n",
      "block 0 consists of orbitals:\n",
      "  up_2\n",
      "  up_0\n",
      "  up_1\n",
      "  down_2\n",
      "  down_1\n",
      "  down_0\n"
     ]
    }
   ],
   "source": [
    "Sigma = SK.block_structure.create_gf(beta=40)\n",
    "SK.put_Sigma([Sigma])\n",
    "G = SK.extract_G_loc()\n",
    "SK.analyse_block_structure_from_gf(G, threshold = 1e-3)\n",
    "for i_sh in range(len(SK.deg_shells)):\n",
    "    num_block_deg_orbs = len(SK.deg_shells[i_sh])\n",
    "    mpi.report('found {0:d} blocks of degenerate orbitals in shell {1:d}'.format(num_block_deg_orbs, i_sh))\n",
    "    for iblock in range(num_block_deg_orbs):\n",
    "        mpi.report('block {0:d} consists of orbitals:'.format(iblock))\n",
    "        for keys in list(SK.deg_shells[i_sh][iblock].keys()):\n",
    "            mpi.report('  '+keys)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This minimal example extracts the block structure by calculating once the local Green's functions and then finds degenerate orbitals with a certain threshold in `Gloc`. Afterwards the result is reported, where 1 block is found with size 6 (3x2 orbitals for spin), where a all 6 orbitals are marked as degenerate. This is indeed correct in cubic SrVO$_3$, as all 3 t$_{2g}$ orbitals are degenerate. Note: for a magnetic calculation one has to break the symmetry between up and down at this point manually. Moreover, one can reduce the threshold for example to `1e-5` and finds that then the degeneracy of the 3 t$_{2g}$ orbitals is not found anymore, resulting in only two degenerate blocks for spin up and down, each with size 3x3.\n",
    "\n",
    "Afterwards the exact same DMFT script as in the [Wien2k tutorial](../srvo3.html) can be used. For a more elaborate example including charge self-consistency take a look at the [VASP NiO example](../nio_csc.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}