AbinitioDMFT/01-solid_dmft.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd835fc6-b6d3-4a2d-b742-a5fa925835df",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "import sys\n",
    "from itertools import product\n",
    "import numpy as np\n",
    "np.set_printoptions(precision=4,suppress=True, linewidth=140)\n",
    "import scipy\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from h5 import HDFArchive\n",
    "from triqs.gf import *\n",
    "from triqs.gf.descriptors import Fourier, Omega\n",
    "from triqs.plot.mpl_interface import plt,oplot\n",
    "from triqs.lattice.utils import TB_from_wannier90, k_space_path\n",
    "\n",
    "import solid_dmft.postprocessing.plot_correlated_bands as pcb\n",
    "import solid_dmft.postprocessing.maxent_sigma as sigma_maxent"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f959ccd-d82a-4caa-b869-91b3685ecc61",
   "metadata": {
    "tags": []
   },
   "source": [
    "# TRIQS / solid_dmft tutorial \n",
    "\n",
    "La2CuO4 is one of the high Tc superconductors and can be relatively well described by a 2D Hubbard model on a square lattice which connects the Cu ions. Further below, we will discuss why this material is well suited for DMFT and that the mapping to a low-energy description is straightforward. For a theoretical overview over DFT+DMFT for La2CuO4 you may be interested in this paper by P. Werner: [doi.org/10.1103/PhysRevB.91.125142](https://doi.org/10.1103/PhysRevB.91.125142). \n",
    "\n",
    "This tutorial should give you an introduction to ab-initio DMFT calculations. It will introduce the important observables / results that are obtained, and how to understand them. Furthermore, you will learn the basic workflows of using TRIQS as a library. One of the key aspects of `solid_dmft` is that everything is written in modular functions in the triqs spirit, but can be run as a grayboxed code with mpirun.\n",
    "\n",
    "During the tutorial we marked specific tasks with the icon <i class=\"fa fa-gear fa-2x\" style=\"color: #186391\"></i> that we recommend you do to deepen your understanding of DMFT, the Green's function formalism, and the usage of TRIQS. All tasks can be done directly in this jupyter notebook.\n",
    "\n",
    "more information about the DFT+DMFT formalism and the utilized codes can be found here: \n",
    "* TRIQS documentation: [triqs.github.io](https://triqs.github.io/triqs/latest)\n",
    "* solid_dmft documentation: [triqs.github.io/solid_dmft](https://triqs.github.io/solid_dmft)\n",
    "* solid_dmft input parameter reference: [triqs.github.io/solid_dmft/input_output/DMFT_input/input](https://triqs.github.io/solid_dmft/input_output/DMFT_input/input.html)\n",
    "* more TRIQS tutorials: [github.com/TRIQS/tutorials](https://github.com/TRIQS/tutorials/tree/3.0.x/TRIQSTutorialsPython)\n",
    "* DFT+DMFT implementation using Quantum Espresso + Wannier90 + TRIQS: [doi.org/10.1088/1361-648X/ac5d1c](https://doi.org/10.1088/1361-648X/ac5d1c)"
   ]
  },
  {
   "attachments": {
    "805b530f-f633-422f-a268-c2ed48b9513b.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc0AAAD5CAIAAADKu/lFAAAACXBIWXMAABJ0AAASdAHeZh94AAAc3UlEQVR4Xu3dfahmR30H8O/ZzWaTCtElL5sQY5SkLTYhG7f1j2izvpSqpCYsUrSYUEgjWQ2FtlZbtEEkmNq0BCohhlUsWlAUJeSfohtoi2KlhQayim1SaSmBlhp8IYbamOzm9I9z79zfMzNnnjlzZs6Zl++HoLvneTsz8/t9n7nnPntv1/c9iIgomX36ASIiioo5S0SUFnOWiCgt5iwRUVrn6AeIaLqu67Qj/A4zKR2rgWgSM1Ld2GLEnCWym5qnPthubWLOUuti5alqpa7r+r53Py37rim8PkutcAefP8+IlLG7ecvGQc9no6IxZ6k21lwLECsBGbjEnKWCWZNrqsUCjoHbLOYsFcAaTFPlk1/yTMyhMXDrw5yljJihE6aghHJscuWRgkZEJuYsrcDMlDA1pY8jcOXBmobcDuYsJWSNjABNhQsDtz78/CzFYQ2Fqaqpxq6L2VnuuY34QpQI97M0jbvnPTEaJuEOt3TMWbKztnQANn9EcjLNBWLgZos52zqzXcOwtxfm2OTKI1yXHDBnW2F2Yxj2bW4cgSsPcuFWxJytjbXZArAti8PAzRZztmDWdpqKXVcfBm5umLMFsHbLVGyqBjFwM8GczYi1GQKwbUgjS8IsMwZuaszZFZiFHoZdQQEcm1x5hNUVEXM2IbOOw7DiKQVH4MqDLL/5mLNxWCt1KhY0rYKBmxpzdhprIU7FeqU8MXATYc7aWessACuSSsTAjav1nLWWUQAWHFVJFrbZLAxcT63krFkiYVhP1CzHJlceYY+Yqs1ZsxQmYa0QjXEErjzIJlKqzVlPLAWiYAxcTzF/6ns+rKsuVTlqykoX9VcqlMLdeg1OyKDy/Wzf99aF5zstUQrc4VpV+JY7dkneuvBKffNA62pzP2vlaL1GpqjCUtj6nulYdYw/imgS5qyp2dar/LqBlXuTO7YdJqKZmr2qUNtbbvBSWRdemfpslLMOJ+Rfe5yUf42F+1kfjfRdbaUQnLNKIwvfLC1kBymiljk7Sd19V1UpzA9Zqe6Fb5M1ZAfRo5Y5G6bKvmvx+qwntaLWhY+b6UQ0kA1ltl6hfVdPzppLEovnwqO0tSfKnGOvU1bf1ZOzStJJdyy8PJj0HIhaU3rfVZKz1tlPqvSFp6S6ruv7fvhf/TaaodC+q6QOcphf68IrK54YSdZvhSX6JhhzdgFF9F0NdSAnOofhuBceeZxky1J/flbGK6N2Me6+W3cVaiiCHDazY3Je+5apqI0esmDOZsDRd6ssR/FFkNtmdoxj4ZH3mdcn6X5WC1bm7Loy6btKvg+GBacsjNzgbN6ycTDzUdBUPa/SriqTvtunH6DE+l36DQCAbpd+AyVj/c5YGC5cttbtu7L3s4u9HaWgztm6uvJgiaPLVsRUtbIuVsctbTY8+y7uepWds3WQK7rk2tMy+pFf6kHrWjJwC75uEHciMtHv0m8AkP6rm2ZF2eRy01oQ1UrLdBP3s5la8s22HVHyNMCwpeViLc/aPpNEWbhS97Pzp68U3OEuY2YER+lGCqYaQaPfz4O748IUv5+NOx05kyM1C0geaWdOytIxi2czKz/MwgtR5MIzUxR32TU+OQ5x/z3YsApbZ5s5689d2P7mTHjEnCl7Pztz8BVQM2CtS3WQE5UaZziAtWjDpJj/Pt4HRcrL2VgjrwwDN39D3za4BNaaDFPo7JWXs0qhM56aZ+CCExhJm9FpZS25MJlMqdrSzlzlwnI24kJWT5aFdd64yaUw1nIK00jtFZazSiPLE4vnJpezOtXUbc6wP5r0kBVZqyVAKeN1m7NwpeYshWHgkslaDAGqLBt16WCOknKWKRARA3e+sA3OWlta60KHWf7kS1dSzlIKnoELdlchrOsYhis+UFva4DfIYnKWO6zU5MRae5VLIFmnaEkRT4ALmloxOUtL8tzkNt6fwcOfdOnAugRhPF+RHPwXTiojZyOWGk3CwF2MdYbDcDmiU5cOwpSRswoLaC0MXClsUzNQc2WdSX/BJ0DLKyBnZ5YjxSXb21waeaTxIDAnJ0zj05gPtaUNeJeNlrOTimriSe6YOjZKrc1Nrmwz68DDVDZLJEXL2aFI4lUdlcSdOxUErjYu6zA9aV8QlDsnLZu6cNFydtD3jNqmlR641tMOk+0YKZi6dDDVtFT253Mynq+cf3OSg7su11pT91lNMmcIw2nMeQZaWFgcpcpZ+EWtNHYiYQOjDLnTLcX6ul/Rn9ynxz3P6E9IScmK8l+4yNcNpOEctDqXJ6bdJP/qff5UEplWm7dsHPQvX8X6hAHcL81MpD7o0sHu77t9EZ/8NZz7Gnzn7MbN0lfehXOP4PQZ/biboyz7fuM/qev2/gPGn4LK1O/SbwAAdLvGjpu0e27X3zn8p87EcT5J9UF9S2WZ8HvFL78CBy7DZfv141uZMWo1lrnDjUC/G7tUD3fAzQ1T8fwa/X4eOm5mCcC2L8us9Jw1H/f8KVxzLd79fnz9aRz4Gb50D97xOhz7GCbua63pOUoErn5GcqvrPUzK3ZwERNQ8JYpuM2cP4B/vw/WX4dzzcPWb8NnvAMCZS3HrcVz0LD73MJ75Bj7/BK48huNHjfzzEFr2HdCNbXWZuZVxp6S8VdLvF1vSzWzPSwfF8ly4je+DnXkcf3YYH/k0LvwxPnk33nMzLj2Ntx3Bh4/gyQfxxUdxy2F89yX40MdxideTx6QKXVa7Nkb512RNQStIl3Funl1ErZn61rj5eYPz8dHP4rYLAeCGA7jqVvzVKbztnXjuMdz2x3jV3fjsUVz/drzn9XjkduOKwxq07pMD1yZhpT6lMoz1jMr3sTvEMvTtWm8nlNpGzu6/Ar/00p0/v/Q1+Pn9+Pfv4Sxw3rW4/xPoj+PQITz0AA7fnEXImrjVJbe9uOzvlEe2BtwyObjAS1Asakvrs2qb+9kXxUeoXkSP3YuwB3Hsjp3DN51Q98gat7otG9t+7u1PsVPHWztkGSru1ZlncmIUxcbG9OxTOP3DnT//6J/xb2dx1S9g+ue4ctSPf2is43fSitWN6Efoj8/JcHrqVOVw9LtSZrau0eZ+9nnc89vo34uLn8GDf4L/fQXueOvG7XWQ7cbLC0UYq+Mlo3PIvsVeUb2QjNrFXp18DCWhH7XZyNlX3YVPXowPnMATP8GVN+AzD+AtF8jbK6TVLS8vrGusauPmi7pokAl3fMubGLiFci1wDmTjLXyqIy0PMHNjsEbqMpmyd3EWJzdv2c6diWECnjPd5NBUPmuR8OfIlI5b3SiseQpnUTYlIGQhZs+nyWl1zFlfsoy16JB/bbbamaerMANXHqQlOd4ymbMhPDMXNcZu43naR/1uWMSnks/DTe6ShpLQj25izs6lVbIjdsuq+bHSYevmz9zkctXWlXXOjrV6zjy3uvmU/dgkszNTi7iZHVN94FqLd8Uhjq1p1jmr+OzMM+TY6i5/eWFsAq1lQQ5xLx0swwxcebAyI5WesMu2BlQZOVuHxba6Y0tea1+VaGyNUpOBW8Emt/f+BdvrDpE5u44oW92xXi23bQoyf0s78+EzyVcvN3BHOmBDDmNizmbBc6uLzZ+uXlxXVGPr14lucx6bgnlVIc/Smjpt+QyCOZsRW/tplbL313xqiALkGWRZBa6lGwz9yHWD9c7arpicHXYQK656RLY8BTxqOuzyAqVQU0GazMCVB6MbaQid5+t73m1Jef7A7koM32ow9SP0xxt6/nTH8g0FoB/NlSxOVcD6nSaSteqoWFnt1po3ffmdOHAdTk/9HbGzbZ0T5mwcqgSlvQTdpD84iKMEfYo4c12h510vWcCqwvU72aRLVWl4yOVX4MBluCy/n5ldzHWDTIzVVqz0DCZfXztH+de1T3Oy1SfWYUicSWc49f55UkNQ7SCODP/vEncC+h7Pn8I11+LIW3DdD3DgZ/jSPfj7U/jBTfi7u3MJOO5n7dQ7tka8qW/QH78qtSMY/pN8dhZEProOQD/856gorRrn9Mr/fRu/dxzXXImXHMRFV+OOk/hJDwBnLsWtx3HRs/jcw3jmG/j8E7jyGI4f3fx0zqoyifvVdNbSyHsbNZUcSk1b3eIM79P60UKMNIpp737RB3vwZXjZ1fij23HVhTj9EP7gd3Hpr+DeX8bPHcGHj+DJB/HFR3HLYXz3JfjQx3GJ7wkvIeuFVyE4nOScMm0hT6camRIgg8zVlj6uOT/n2+Rflv73XJ2jNiTHaGTHpRj1mcdx9LV4+Un8ze+gA557DDe+AfvuxteO4vq348hJPHL7cl+tbx1sbftZ5qk/bUrkzGmzyMmbL+eQHWmaDVPPXQ421rvm09/ER+7FqcfwPz9B3+OFszj8ws5N512L+z+B/jgOHcJDD+DwzcuFrI9Sc5Z5Gp2cOV5emCTnDNWM9I0u7mjU5MwJ3Bf/E7f9Bh57Db7wLbzpKuw7jaOvFTcfxLE7dv54U16//g3IOWetS2I9SNFxq+uv9/g3uFvvkI7PKy+2iGbgyoNuL3wb//Qsfv19eOvVAPDT/8b3X8Rh/V6ZyjdnrTyXhOLiVne+ZUo3q1R1kIHruXnafyVeeQ4e/yr+5Y04+D386YfxI9fd85LVRQwqQD/+AZ2u1Q+Nube0jpvm2zrncr2sq7auXhgyd2y6zrkOJ+/DBadw5OV4/V244mP4/V/U75OtfK8rme9yXTlXwdo00iDA9D2UufoRxf28wcBRnI6bpnLMsBLppdaUdPWj094YrOdc2HUDypksMC0R5F9tdUh2W4O1ysmUWyvtSP6sb6vMWUpCqzSZF1p2GDVZpOHLXrPBrAfHtJmqDmbgyoMFKSZnJ9Ur5YZbXSsGqycZuCVucgvIWTWtjNo6+G91S6dV7NbRsbq32pzPYgK3gJyVCppZ8mTb6qpA6lU25b/gtksHe+dvlf+gcmZeVcg2FgrLWSX/maUAw2Jasyn/ywtdB3ew5nnaFcg/cPP9SlybMnkh3CrbgdBU1m5xrP+klVcf6hrM+WiX45QGk06MIpJxkToZzGgyX7GknDXftayyHRF5suas5Fj/kUfs0EJ24Bm1jhcdDC89FKpxAYHWsUDgmnFkvlC+1eDIWckcpGR9CGVua85qHCUgn8AasgMzah3PqVjPjjmbral15cmMIPP5860Gz5xVzNFK7sdSVub0w2gV9KMhu6PTo1bjfy7D+QecPC1jToFp5FM5nrbU74OZ1NisgeuYAqqJtry2WvA1p1JYZjkzsyL1etWTs4o5idJiM0s5UItsKYVNLIcGWbMiRTJUmLOKdRIVBm6erIu1gBQ/XIZKIUMgRTKUkbPdtouzbtZJtB6Z8yoU3cx1D8PAbZy5P5tfhGXkbETmJEoRZ5YCWBdlLQzcxplZERwLK+wXPMmxJd3XuHs73euSyboWM5dg7+Nc3Unrpw5kjDo++wXvwE1arrQi82tfLabkTVK+BbFYzirWJlcWOAHyrFp/GyE72Izaseh0By7GHwjmbAPMrHBXbL4FsXzOKuYkapY8mXZoZeqoWh8bQSlC1pGPY9yZaz7hwuVKK7Jmhbn6+RbEijkrWedRWeusqjSWs+qIJz0WN3eyZiz605950/DMK9YqLU+LKfXnvXsAYM76Y+AmZa1R60EHew4OOTs7ZCX7C+2K9SqUP+ZsKgzcFKw1aj24leWKwe412egJyMBtmVafjnLNKL80WpxleJ4M3FjkTGrz5qhdh534i3fFYCt34CLxq9Mq/HM238/PDttY/WhO1Gxaz9ORHTQm5kSNfMAgERWjXdehv3PzRkAEMQO3QfnuZ+F8f8iWNXOVggayGPcbUkANdDhhDdllAq4T17jcm9xlzofS8d/PMmdTYeB62rrKW++g2Zv4/kSPk0PYLRZqMmf3DjJwK2K2NnN2feaqSEUPbT73ZnYwqQz2NrPio7IdQj42G8aaswoDtyzu5lWYsxlxr1kdY5zKZ4l9snhgDdmFuXNWcQcumLnLcvfmVr3zQ11gzq7Fva71jddqQoB6VIK8LDt+r+Q8c1ZyZy4DNyJ33/kzN7DuEp1cE0tyn3od3Atf8cAxZX233jOTkEVQzioM3FjcbeXPvZTM2fK4K6PKGfBfXzk51jur2203LmpOzioMXE/urvEXtmTM2YK5S6eaqZi6uI7753BZVomSs8qcwF34gxZJuZvCX8SlAXO2Dltrq+hpmbq4Y/fPKmQRO2cVd+DClqfqIeZNOdta9j5SLIGJOVsbd/EVNz9hK2s+KsOQxcRBBXBnrvw0m3kwH+6S9pR6qt2Ys9VyV2cpExW2stqj8vnel5JoMzvGHbimVdLWXbH+lpxYT8zZ+rnLN+cZk2c+6Tw3htzfmVvIYvGcVfwDN13UugvS3yoTGIY52xB3fWc4dXOWVQx297F5XDEYrJWzik/gWi8syONu7nrzt+5ExWLOhnVcK5eF25yGbJC55JocpjF4MzvQczankEUGOav4BK5JTubWcvKUyYQkYs6Sdbz79ANUrF7QbwMAdLv0G9YwdpLeutxCtgIdTgQXiSw/Sb9fk5izdXJXeXAvzRThFcWPdmXIjhndzHafQvcp/eBEIkI36PcjId+f801RqAawZpw6uHCfhL3c7gcM5iZFu2w/gHxM2BqRVRk522Vz2atcqweu9XWn2fmAQQ90aLsqRnesszU7pUllnbN99r+6pkSykczplUcStVzY04pvgJ2Afta1SRejAPrdn31O0Y3Vdu47ggX2WQRb4EpRJn/OUmr/7mvOUyXlv8VOmnQ+V66tJ+DzQJK0xhlbfd+yWMsC2yuSEgXunGTUQnbn4IwnTEfmrDXIouk+hdljD/v8LEmV5Cxy7ajqxQ3c4EXcC1lAPjT4CSNaIEkHK46R3JizFMf8wA1ewbGQxeZZTX1aT4vFqJJoIJROPTmLGY1KcYVlbvDy7b3a5i9V3L15J6qmPq2yfJIOgk+YcuOZs1l/3oByo8rIGrjWPLXec4wRfBuXZY1btwt4yATjSaqMNR41pbCc7by/n0tJBQTu1oWzZGJ/Atj+OaSu6yZ9An8CjyQd9H3P4qQxxVSGtXUpH9bAldwL507SUSoHw3LWO0YVxyiYsw3idQNalHuHO1Z/yU1P0sFqJ0w1KukdmFva4sjMda9a4H4W4Uk6cJ/VJNzPNoj7WVpfn82/nB5rAKIFFJmz3Dg0xfKJLsW2mWVt0GI8dxKFBRYvHRRn5qWDsX8MunPn2R+hncnssbXOhNbiU+HMWUrOZ9V27tJP+xf3Ps8chZmnA/N1FzslyoTPiheWs/AbFWVl65KpEBu53W7r0wbwz9OtfLY5VAGfOizy+iyVpd+9htXZLqyHhexMEfN0jHqqrut8WpEqVnDOWpuWyhIcsmNBqRm725KVI1+Lgdum8nJWbY6oRNZ3R+PABMOzjZWE+Vrrkptc7QhVrLycpRJZ3x2HA/45Yz6DOlhcWjFwm1J2zlo3R1QEW2busObpQNu9VrD6ZuDKg1SHInPWujmizKlV67oOUDliWck2U0aOuqY3EkKhOUtlMaJ0Jzv6fu/P/lrIIHOTW/Fgi2bUtl2pOSs3RyzBTIzV3GZqyJClLRi4pXCvS6k5Syvamqdj+r7ffai8dDBBs3FjBq48SKsYawRTDTnLLW0iY2UUPNsiZCmQDNxm33Vys3X+C85ZdemA5nDM4dbqmUR7nYB3R8epNkjOHgM3cwXnLE0yFlLLdKZ68b7XAzfAMudcEPOqAqcotUlTXXbOqi1twOaoVuvmqZUM2c3jE1ZtbFwkMXDzVHbOtmwsd7LtK3Ve6t0xTLYDzIoZuPIgzTT1bayenJ20OSrIWCSVMtjh9As52QrJwJ2aDhRL8Tk7c3OUj7FRFN0SI2Pa4/nuyICYT04d53OOsVZ1KD5ny+JYofoqXo3VHFk1746FMq8q1Fd+C/CftBpyVjWt5+ZoAY4QyeQMU3OE7FSOyaSZGLhThVVjDTm7orFJb7xYfUI24N3R824UwAxceZBMkyantpz1b9pJmKcBoszN2MxTIjJwucnVBFdjJTmrNkczjT0J68zfMIWTJszn3XHrHSguOeEMXM3UeagkZ6dinkbUdXupOjKvdrHeHSk186pCg50yp1YrzFm5ORqbmgarJKkhatVkx5rdlrs6T9bOam11Asa7/eu1gpipWtPo8mRM+eSQdYSp4ybKRDuBO6caK9zPSu0UQSYizvGcsqbFyNWpeMlmDq2qnLV+UaPMnCkymdM8HAmbYHnBh0pkNiAXdFB5ZVsDV6l77AtwzO6kqTXbkl+I1MFc2RLNH0VV+1mTtW+V+dNHVhGnk0tTNGsDNrimle9nTdbAVVqbjTnGJjJsCuV7XuM9WbfiFjfKCVe+nzXJmTIzlzvcOaLPGVehPtYGLGKh55xkc/tZKzNwJU6RyZywmZNkXQLOfCOyDdwom1kwZzXWblc4V4qcp1izok0+Z7tBuQVurPNp7rqBm5pNa+DGmvSacCYoImsDrtVu1hAIw/3sFu65bnP2himJPnTuZ8m0YuBG3FcxZ30xcJVO/OCYuCJWNtVnyfKIm++8buBLzrWZuUtWwOoaGCLlyLyqsEC7RXkJ5mwIc72VuG+DRGQyGzDzXuN1gzjMwJU4yZ5KaRvKTdz9TfQ65H42DvMNVoq+bEQkWS/r5dNu3M+mYg1chdNulWGHULnCyinsUW7cz6bCHS7RusweXKvd9ukHKLZe0G8Dul36Dev5yrtw7hGcPqMfJyqUakDVbmMdN3Z8Ju5nF6Wi1lxOecSayIu5/AoceAaX7dePE5VONqDqOGu7WQ8G4352HeoNVr8BwEqb3OdP4Zpr8e734+tP48DP8KV78I7X4djHwH0t1Uc1YO+xyZ2PObsytdj6DQCSBu4LePQ+vOnVuOB8HLocxz6A7x/Grcdx0bP43MN45hv4/BO48hiOH0WC1x6l5iHJkKl5e4Eq6HeKvZkF+HmD/FgXXom1Xl//Q7z1L3H9HXjfm/H8U3jqFfjob2E/8OSDuPHPccNhfPc6fOvTuMR1Lkmo4ccaKTXI3URbRa895my+3LUyZ+H6H+M3X4m/fSOefASHxYs89xhufAP23Y2vHcX1b8eRk3jk9qW/5GHOkj93j/hLXWz8Pli+5Nqb9TQnj87+B574KV79q7ho81nPuxb3fwL9cRw6hIcewOGblw5ZIiuz/sMENEsUzNkyOC5cyiOTyqjbZ1x7PYhjd+z88aYTG7cQLcAs7zCTGmEB3K8Upt+l3wBAXObXb9i0/5W46iD+9R/wQ/vTECWnalXS77SNageNfr+18fps8dzVaV/fHl99L275DG64C3e9GXga33wGf/FBnK/fbwVqOPYzp9K469NT6cXAnK2Hu6D1hX4Oj9yLe/8a3/kvnHMx3vxBfPn9OLhxj3UwZwvlLj9/Va47c7ZC7orPfMWZs5lzV5e/ptaX3werkKxgsysYZOTDrJwwLDMwZ6unqtxsG3mEzdAyszYCsIQcmLOtcASuPMhuqZh16adihQRgzjaHgVs968oGYA3EwpxtFwO3dNaFC8AlTo05Swzc3FnXJQBXcC3MWdoj+9DsbQZuauach+EC5YY5S3aOTa48wpYOY85qAE5+KZiztIUjcOVB9ryVddKm4tyWjv8ejCZzZ8f8iioxu91z4q+gIZM/7mdpspZ3uNYhB6hycmgMc5bCVRy41hEFKHHsFB1zliKQaWImVOaBa55wgDyHRplgzlJkjk2uPLJ8MJnnE2D506YKMGcpFUfgyoPRk8v6cgGinxg1izlLySUKXOuzBZj6ukRT8XNdtALPiByK0/POW7HUaS3MWVpTrAyVWNKUG+Ys5WJq5rJ0qRTMWcqOGbisUioac5aIKK19+gEiIoqKOUtElBZzlogoLeYsEVFazFkiorT+H1fEVgHfUk9KAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "id": "2ab1e52e-0fe3-4f99-9fb0-3a03209df8e0",
   "metadata": {},
   "source": [
    "## 0. Performing a DFT Quantum Espresso calculation and creating a Wannier Hamiltonian\n",
    "\n",
    "(This step is a recap of running DFT and creating Wannier function. Due to the lack of time we provide all the output here, but feel free to repeat the steps.)\n",
    "\n",
    "The aim is to perform a SCF DFT calculation with Quantum Espresso (QE), to identify a correlated subspace (here: Cu $d_{x^2-y^2}$ orbital), and to create a Wannier Hamiltonian for this orbital. The results will be converted to a TRIQS-readable h5 format, via the [triqs.github.io/dft_tools](https://triqs.github.io/dft_tools/latest/) package. The input for the QE SCF calculation can be found in the folder `data/scf/lco.scf.in` , and consecutively a NSCF calculation has to be performed to obtain the input for Wannier90 (W90) `data/mlwf`. We also perform a bandstructure calculation `data/bands/lco.bnd.in` to compare our Wannier fit with the DFT bands. Last, we construct a Wannier Hamiltonian using W90 (`data/mlwf`). In this specific case we construct one Wannier function for the Cu $d_{x^2-y^2}$ orbital with disentanglement. The input for W90 looks like (`data/mlwf/lco.win`):\n",
    "\n",
    "The W90 input contains the following important parameters:\n",
    "\n",
    "```\n",
    "num_wann = 1\n",
    "num_bands = 5\n",
    "mp_grid 7 7 7\n",
    "\n",
    "exclude_bands : 1-31,34-40\n",
    "\n",
    "dis_win_min = 10.0\n",
    "dis_win_max = 16.0\n",
    "dis_froz_min = 12.3\n",
    "dis_froz_max = 16.0\n",
    "```\n",
    "\n",
    "The latter is necessary to properly disentangle the target WF from the other bands. \n",
    "\n",
    "The bands in QE and W90 will be plotted along the following path (green line) in the BZ:\n",
    "\n",
    "* G:  [ 0.00,  0.00,  0.00]\n",
    "* S': [ 0.25, -0.25,  0.25]\n",
    "* X:  [ 0.00,  0.00,  0.50]\n",
    "* Z:  [ 0.25,  0.25, -0.25]\n",
    "\n",
    "![Screenshot_20220503_125745.png](attachment:805b530f-f633-422f-a268-c2ed48b9513b.png)\n",
    "\n",
    "\n",
    "# 1. Creating a low-energy description of La$_2$CuO$_4$\n",
    "\n",
    "In this section we will discuss the process of constructing a low-energy model of the material.\n",
    "\n",
    "### DFT density of states\n",
    "\n",
    "Let us first take a look at the density of states to identify the character of the states close to the Fermi level:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32e3bc5f-8cdd-4398-805e-aa68fcc3eda1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read partial DOS from QE\n",
    "loc = './data/scf/dos/'\n",
    "seedname = 'lco'\n",
    "fermi = 12.7367\n",
    "\n",
    "dos_data = {}\n",
    "elems = {'La': range(1,3), 'Cu': range(3,4), 'O': range(4,8)}\n",
    "char_table = {'La': {1: 's', 2: 'p', 3: 'd', 4: 's', 5: 'p', 6: 'f'},\n",
    "              'Cu': {1: 's', 2: 'p', 3: 'd', 4: 's', 5: 'p'},\n",
    "              'O': {1: 's', 2: 'p'}}\n",
    "\n",
    "with open(loc + f'{seedname}.pdos_tot', 'r') as f:\n",
    "    dos_data['tot'] = np.loadtxt(f)\n",
    "\n",
    "for ct, (el, n_atom) in enumerate(elems.items()):\n",
    "    n_char = char_table[el].keys()\n",
    "    for atom, char in product(n_atom, n_char):\n",
    "        with open(loc + f'{seedname}.pdos_atm#{atom}({el})_wfc#{char}({char_table[el][char]})') as f:\n",
    "            dos_data[f'{el}-{atom}-{char}{char_table[el][char]}'] = np.loadtxt(f)\n",
    "\n",
    "# read the w90 DOS for comparison\n",
    "with open('./data/mlwf/lco-dos.dat', 'r') as f:\n",
    "    w90_dos = np.loadtxt(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22ec8879-9937-48d6-b114-7487efad7131",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, dpi=150, figsize=(8,3))\n",
    "\n",
    "# plot total density of states\n",
    "ax.fill_between(dos_data['tot'][:,0] - fermi, dos_data['tot'][:,1], zorder=2, color='gray', lw=0., alpha=0.3, label='total')\n",
    "\n",
    "# plot partial DOS. The entry number corresponds to the orbital order in data:\n",
    "# 0 col is the energy, then the summed partial DOS, the columns are dz2, dxz, dyz, dx2-y2, dxy\n",
    "Cu_dx2y2 = np.sum([dos_data[f'Cu-3-3d'][:,5]], axis=0)\n",
    "ax.plot(dos_data['tot'][:,0] - fermi, Cu_dx2y2, label=r'Cu-$d_{x^2-y^2}$', zorder=10)\n",
    "\n",
    "Cu_dz2 = np.sum([dos_data[f'Cu-3-3d'][:,2]], axis=0)\n",
    "ax.plot(dos_data['tot'][:,0] - fermi, Cu_dz2, label=r'Cu-$d_{z^2}$')\n",
    "\n",
    "O = np.sum([dos_data[f'O-{atom}-2p'][:,1] for atom in elems['O']], axis=0)\n",
    "ax.plot(dos_data['tot'][:,0] - fermi, O, label=r'O-$p$')\n",
    "\n",
    "# plot the w90 DOS\n",
    "ax.plot(w90_dos[:,0]-fermi, w90_dos[:,1], label='w90 DOS')\n",
    "\n",
    "ax.set_xlim(-8,2.8)\n",
    "ax.set_ylim(0,4)\n",
    "ax.set_ylabel('DOS')\n",
    "ax.set_xlabel(r'$\\omega$ (eV)')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f443cbe2-29d9-445a-9a02-d010cc570e55",
   "metadata": {},
   "source": [
    "It can be observed that the states right at the Fermi level are dominantly of Cu $d_{x^2-y^2}$ (blue) character hybridized with O $p$ states (green). The W90 DOS matches exactly the total DOS in the window given by the Cu $d_{x^2-y^2}$ states, capturing the nature of the bonding orbitals as well. This means that this is the best compromise of having a localized Cu $d_{x^2-y^2}$-character-like state, but still describing properly the electronic states corresponding to the band close to the Fermi level.\n",
    "\n",
    "<i class=\"fa fa-gear fa-2x\" style=\"color: #186391\"></i>&nbsp; What is the reason for the contributions to the DOS of the Cu $d_{x^2-y^2}$ (blue) states at low energies between -8 and -5 eV?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c7ccc1d-85d2-41bd-bf2f-875eefca7bdc",
   "metadata": {},
   "source": [
    "### Plot QE bands together with W90 bands to check quality of Wannier fit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32e3e4b0-759b-4145-9ba2-1b017bdfd817",
   "metadata": {},
   "source": [
    "To check more carefully whether the tight-binding Hamiltonian describes the DFT bands close to the Fermi level correctly, we now plot the QE and the W90 bands on top of each other. To plot the W90 bands we will make use of the solid_dmft postprocessing module [plot_correlated_bands](https://triqs.github.io/solid_dmft/_ref/postprocessing.plot_correlated_bands.html) (`pcb`). This utility reads the real space Wannier Hamiltonian, and then leverages Wannier interpolation to construct $H(k)$ on any given $k$ point mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2b75b37-c21a-4fd4-9b4b-8c5a0ebea5ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "# w90 information\n",
    "w90_dict = {'w90_path': './data/mlwf/',\n",
    "            'w90_seed': 'lco',\n",
    "            'n_orb': 1,\n",
    "            'mu_tb': fermi}\n",
    "\n",
    "# plotting options\n",
    "plot_dict = {'colorscheme_bands': 'coolwarm'}\n",
    "\n",
    "# tight-binding-like bands we want to plot from the given W90 Hamiltonian\n",
    "tb_dict = {'bands_path': [('G', 'S\\''), ('S\\'', 'X'), ('X', 'G'), ('G', 'Z')], 'n_k': 50,\n",
    "            'G': [ 0.0,  0.0,  0.0], 'S\\'': [0.25, -0.25, 0.25], 'X': [0.0,  0.0,  0.5 ], 'Z': [0.25,  0.25, -0.25] }\n",
    "\n",
    "# get_dmft_bands extracts H(k) from the above specified path\n",
    "tb_bands, alatt_k_w_path, freq_dict = pcb.get_dmft_bands(**w90_dict, **tb_dict, eta=0.01)\n",
    "\n",
    "# also load the bands.dat.gnu file for a comparison\n",
    "band_data = np.loadtxt('./data/bands/lco_bands.dat.gnu') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "203e6b10-cfa5-414b-9262-f02e7653f27d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,dpi=150,figsize=(5,6))\n",
    "\n",
    "# since every code uses a different scale for a k-path we have to scale the QE bands\n",
    "factor = tb_bands['k_mesh'][-1]/band_data[-1,0]\n",
    "\n",
    "# plot QE bands\n",
    "ax.scatter(band_data[:,0] * factor, band_data[:,1] - fermi, color='C0', s=4, label='data')\n",
    "\n",
    "# plot W90 bands using the plot_bands function of pcb\n",
    "pcb.plot_bands(fig, ax, alatt_k_w_path, tb_bands, freq_dict, n_orb=w90_dict['n_orb'], tb=True, alatt=False, **plot_dict)\n",
    "\n",
    "ax.set_ylim(-8,2.5)\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "578b9623-5c80-4396-a126-3c3d434dbad9",
   "metadata": {},
   "source": [
    "### Convert Wannier90 to TRIQS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28b31c97-97f7-4a86-94f0-c4b2febf948b",
   "metadata": {},
   "source": [
    "Now, we will convert the W90 output to a TRIQS readable format using the Wannier90Converter of [triqs/dft_tools](https://triqs.github.io/dft_tools/latest/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c62fd2fe-dd30-49ce-bd78-9a7d86a99f4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from triqs_dft_tools.converters.wannier90 import Wannier90Converter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53f4f897-73ca-4efc-9bea-a7321a16c231",
   "metadata": {},
   "source": [
    "The wannier90 converter input `lco.inp` contains the following information:\n",
    "```\n",
    "0   7 7 7              # specification of the k-mesh\n",
    "1.0                    # electron density\n",
    "1                      # number of atoms\n",
    "0  0  2  1  0  0       # atom, sort, l, dim, SO, irep for each atom\n",
    "```\n",
    "\n",
    "Now we run the converter in the according sub-directory producing an h5 archive `lco_wannier.h5`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96acd28c-c476-448f-98bd-6df3fccd790b",
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd \"data/mlwf\"\n",
    "Converter = Wannier90Converter(seedname='lco', hdf_filename='lco_wannier.h5')\n",
    "Converter.convert_dft_input()\n",
    "%cd \"../../\"\n",
    "# and copy the h5 archive in the cwd\n",
    "!cp data/mlwf/lco_wannier.h5 ./"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8d28595-2475-4b0e-88b1-b6d32cf4c0a3",
   "metadata": {},
   "source": [
    "### The HDFArchive (h5)\n",
    "\n",
    "An h5 archive is a standardized storage format similar to a zip archive, but for directly storing arrays, strings etc. and is interfaced to many different programming languanges, which immensely simplifies data transfer across different codes. TRIQS provides its own h5 interface layer build on top of h5py making reading / writing very easy. After opening an h5 archive, you can access the object as if it were a dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cc08af5-f1ec-4c23-b256-c47bc25c1648",
   "metadata": {},
   "outputs": [],
   "source": [
    "with HDFArchive('lco_wannier.h5', 'r') as h5:\n",
    "    print(h5['dft_input'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3b5b682-3cec-4870-a314-1397d1d26427",
   "metadata": {},
   "source": [
    "As we can see, the Wannier90Cconverter calculated $H(k)$ on a predefined k-mesh and stored it with additional information in the archive. The information stored is standardized for all DFT codes in TRIQS (see [triqs.github.io/dft_tools/3.2.x/h5structure](https://triqs.github.io/dft_tools/3.2.x/h5structure.html)) and interfaces to other codes exist. This is very practical, since we can then run solid_dmft starting from any compatible DFT code."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1672719a-329c-43e7-8030-7a30d4788c6f",
   "metadata": {},
   "source": [
    "## 2. Run a one-shot DMFT calculation using solid_dmft"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f389774-f621-4c0d-8ce2-dce67658f9bd",
   "metadata": {},
   "source": [
    "Now that we have finished all preparations, we are able to run the one-shot DMFT calculation. We will use the [triqs/cthyb impurity solver](https://triqs.github.io/cthyb/latest/) to solve the impurity problem within the DMFT loop. This is a Monte Carlo solver, sampling the impurity Green's function $G^{imp}(\\tau)$ in imaginary time $\\tau=i t$ at finite temperature $T$ by expanding the impurity Hamiltonian in the hybridization function, which couples the bath with the impurity. We will perform our calculations at inverse temperature `beta = 10`, which corresponds to $T \\approx 1160$ K, where $\\beta = \\frac{1}{k_B T}$. Higher temperatures are numerically easier. The solver is fully MPI parallelized, and written in C++.\n",
    "\n",
    "The accuracy of the solver is ultimately determined by the total number of Monte Carlo (MC) moves it performs and how many successive measurements of $G^{imp}(\\tau)$ are made. This is controlled in the `[solver]` section of the solid_dmft input file with the following parameters:\n",
    "\n",
    "```\n",
    "n_warmup_cycles = 1e+4\n",
    "n_cycles_tot = 2e+6\n",
    "length_cycle = 60\n",
    "```\n",
    "The first parameter `n_warmup_cycles` is important to thermalize the system before any measurement is performed. This task is performed on every MC \"walker\", i.e. on every MPI rank, and the number strongly depends on the problem at hand. For simple problems at higher temperatures a value between 5000-10000 is usually enough. The next parameter, `n_cycles_tot` determines the number of measurements performed. The number is absolute for all MC walkers, and the workload of performing the measurements is divided among the MPI ranks. The last parameter `length_cycle` is equally crucial for accuracy. This is the number of MC moves performed between two measurements. So the number of total moves will be `n_cycles_tot * length_cycle`. The length of a cycle is important to decorrelate two measurements. If the length is not large enough, measurements will depend on each other and general accuracy will go down.\n",
    "\n",
    "Since we have a limited amount of resources we will first perform 6 DMFT iterations with $2 \\times 10^6$ measurements to pre-converge our results. The cycle length is optimized for this temperature and problem. No need to modify. A detailed tutorial on how to find the optimal parameter set for cthyb can be found in [this tutorial from cthyb](https://triqs.github.io/cthyb/unstable/guide/cthyb_convergence_tests.html).\n",
    "\n",
    "solid_dmft is controlled via a single configuration file. This file is usually called `dmft_config.toml` but can be named differently when the name is passed as argument to solid_dmft. The file is divided into sections: `general`, `solver`, `dft` and `advanced`. It follows the typical [toml](https://toml.io/en/) file specification. Let's have a look at the config file we provide:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b470942-1675-4917-a2ab-587d18c48e6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat dmft_config.ini"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a22faaa-ca3f-4a7d-b26e-dfbe5068b963",
   "metadata": {},
   "source": [
    "Most important parameters not yet introduced are `U` (similar to DFT+U), the choice of solver `type` of solver (cthyb in our case as introduced above), the number of DMFT iterations `n_iter_dmft`, the jobname where the calculation is performed in `jobname`, and the choice of double counting `dc_type` (similar to DFT+U). For now you can leave all parameters as is.\n",
    "\n",
    "The reference for these parameters can be found here [triqs.github.io/solid_dmft/input_output/DMFT_input/input](https://triqs.github.io/solid_dmft/input_output/DMFT_input/input.html). Let us know if you have any questions.\n",
    "\n",
    "<i class=\"fa fa-gear fa-2x\" style=\"color: #186391\"></i> Now go to the terminal, navigate to the directory of the tutorial and start the run either interactively (the calculation will take 3-4 min on 4 cores): \n",
    "\n",
    "`>>> mpirun solid_dmft > lco.dmft.out 2>&1 &`\n",
    "\n",
    "and inspect the output with\n",
    "\n",
    "`>>> less lco.dmft.out`\n",
    "\n",
    "The calculation will create a directory `b10-U3.6`, where everything will be stored into the h5 archive and text files. Take a look around and see if you can identify the important steps of a DMFT calculation. Take special attention to the impurity density matrix printed before and after the solver is running. This is a good indicator for how well the calculation is going. In- and output density matrix should be the same at convergence.\n",
    "\n",
    "During the calculation you can also monitor the `observables_imp0.dat` file, which is stored in the `jobname` directory. This file summarizes important observables in each iteration. Make sure that the chemical potential `mu` converges to a fixed value while the calculation is running. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5d7954d-cd3b-4f47-8597-40f3823b0e83",
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat b10-U3.6/observables_imp0.dat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e5b78aa-39d2-4938-a97e-90f526f4b6eb",
   "metadata": {},
   "source": [
    "## 3. Check convergence"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "683afa7d-7e04-4ad7-a517-e6ae8bd6728e",
   "metadata": {},
   "source": [
    "After the calculation is finished we can take a look at the results (this step also works while the calculation is still running) and plot the results. First, we load the results from the h5 archive in the `jobname` directory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "893a37d2-50b9-42a3-aa10-6b8138d2b714",
   "metadata": {},
   "outputs": [],
   "source": [
    "with HDFArchive('b10-U3.6/lco_wannier.h5','r') as h5:\n",
    "    print(h5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46004291-7c0e-4d59-b95d-ee10aa2f97fb",
   "metadata": {},
   "source": [
    "The h5 archive now has additional `groups` for the in / output of the calculation. You can print each sub-group to view how the data is stored. The organized folder-like structure allows us to later read the calculation again and continue where we left off. Let us now load the results of our calculation at the current stage:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ab1733f-63e9-4844-aeec-7db01d7ac0fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "with HDFArchive('b10-U3.6/lco_wannier.h5','r') as h5:\n",
    "    # the next two lines load all directly measured observables and the convergence metrics per iteration\n",
    "    obs = h5['DMFT_results/observables']\n",
    "    conv_obs = h5['DMFT_results/convergence_obs']\n",
    "    # The next lines load the last iteration impurity Green's function and self-energies\n",
    "    Sigma_imp_iw = h5['DMFT_results/last_iter/Sigma_freq_0']\n",
    "    G_imp_iw = h5['DMFT_results/last_iter/Gimp_freq_0']\n",
    "    G_imp_tau = h5['DMFT_results/last_iter/Gimp_time_0']\n",
    "    # last, the chemical potential and double counting potential\n",
    "    mu = h5['DMFT_results']['last_iter']['chemical_potential_post']\n",
    "    dc = h5['DMFT_results']['last_iter']['DC_pot']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0266f20b-63f8-425d-a0da-b402bc558510",
   "metadata": {},
   "source": [
    "We now plot the data and take a look at the most important observables to gauge if the calculation is converged: \n",
    "* the chemical potential\n",
    "* the impurity occupation\n",
    "* the spectral weight at the Fermi level, which is directly accessible on the imaginary axis, without analytical continuation (see also [ref](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.83.235113)):\n",
    "\n",
    "$$\\bar{A}(\\omega=0)= \\lim_{\\beta \\rightarrow \\infty} -\\frac{\\beta}{\\pi} G(\\tau=\\frac{\\beta}{2})$$\n",
    "\n",
    "* the quasiparticle weight Z, i.e. renormalization of the DFT bands:\n",
    "\n",
    "$$Z = \\left[ 1 - \\frac{\\partial \\text{Im} \\Sigma(i\\omega)}{\\partial \\omega} \\big|_{\\omega \\rightarrow 0} \\right]^{-1}$$\n",
    "\n",
    "* the convergence of the Weiss field $\\mathcal{G}^0(i \\omega_n)$\n",
    "* the convergence of the DMFT self-consistency condidation $|| G^{loc} - G^{imp} ||$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa96e2b9-88c4-4e5f-9a9d-b9ea42c88a50",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=6, dpi=120, figsize=(7,10), sharex=True)\n",
    "\n",
    "# chemical potential\n",
    "ax[0].plot(obs['iteration'], obs['mu'], '-o', color='C0')\n",
    "ax[0].set_ylabel(r'$\\mu$ (eV)')\n",
    "\n",
    "# imp occupation\n",
    "imp_occ = np.array(obs['imp_occ'][0]['up'])+np.array(obs['imp_occ'][0]['down'])\n",
    "ax[1].plot(obs['iteration'], imp_occ, '-o', color='C1')\n",
    "ax[1].set_ylim(0.98,1.02)\n",
    "ax[1].set_ylabel('Imp. occupation')\n",
    "\n",
    "# A(w=0)\n",
    "Aw0 = -1*(np.array(obs['imp_gb2'][0]['up'])+np.array(obs['imp_gb2'][0]['down']))\n",
    "ax[2].plot(obs['iteration'], Aw0 , '-o', color='C2')\n",
    "ax[2].set_ylim(0,)\n",
    "ax[2].set_ylabel(r'$\\bar{A}(\\omega=0$)')\n",
    "\n",
    "# Z\n",
    "Z = 0.5*(np.array(obs['orb_Z'][0]['up'])+np.array(obs['orb_Z'][0]['down']))\n",
    "ax[3].plot(obs['iteration'], Z, '-o', color='C3')\n",
    "ax[3].set_ylim(0,1.0)\n",
    "ax[3].set_ylabel(r'QP weight Z')\n",
    "\n",
    "# convergence of Weiss field\n",
    "ax[4].semilogy(obs['iteration'][1:], conv_obs['d_G0'][0], '-o', color='C4')\n",
    "ax[4].set_ylabel(r'dG$_0$')\n",
    "\n",
    "# convergence of DMFT self-consistency condition Gimp-Gloc\n",
    "ax[5].semilogy(obs['iteration'][1:], conv_obs['d_Gimp'][0], '-o', color='C5')\n",
    "ax[5].set_ylabel(r'|G$_{imp}$-G$_{loc}$|')\n",
    "\n",
    "ax[-1].set_xticks(range(0,len(obs['iteration'])))\n",
    "ax[-1].set_xlabel('Iterations')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "237d2a2c-7705-4d0f-a59a-488d17b1bb56",
   "metadata": {},
   "source": [
    "All metrics seem to indicate good convergence in the first few iterations. However, for the small number of measurements performed some metrics like the impurity occupations show room for improvement. The last two plots show overall accuracy of the calculation. Note that the Weiss field $\\mathcal{G}^0(i \\omega_n)$ generally converges first, and only after a few extra iterations the convergence of $|| G^{loc} - G^{imp} ||$ will follow.\n",
    "\n",
    "<i class=\"fa fa-gear fa-2x\" style=\"color: #186391\"></i>&nbsp; Why is the QP weight larger than one? Understand how it is computed in [solid_dmft](https://github.com/TRIQS/solid_dmft/blob/unstable/python/solid_dmft/dmft_tools/observables.py#L590) and compare with self-energy plot below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "136e027d-961c-4b53-b2a9-7759376d32bf",
   "metadata": {},
   "source": [
    "### Reminder: The TRIQS BlockGf object\n",
    "\n",
    "`Sigma_imp_iw` loaded from the h5 archive is a TRIQS Block Green's function object ([see doc here](https://triqs.github.io/triqs/unstable/documentation/manual/triqs/gfs/py/full.html)). This is a structured container holding not only the raw data of the Green' function like a numpy array, but also information about the mesh, and its block (i.e. spin/orbital) structure. So when accessing a Block Green's function object you only get back information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9968b91-b45b-42e1-87a7-d55d628a436f",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(Sigma_imp_iw)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82570137-7ce7-4ddb-ae7f-243739f53db0",
   "metadata": {},
   "source": [
    "we see that `Sigma_imp_iw` consists of two blocks. Each block has a name and is printed as separate line. The block names itself can be accessed via the `indices` property:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5252efe6-3d1a-4a13-84af-eb81b9234cca",
   "metadata": {},
   "outputs": [],
   "source": [
    "list(Sigma_imp_iw.indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d407e2c-e585-4d58-b9c7-7992b1d8d17f",
   "metadata": {},
   "source": [
    "We have one block for each spin index. Since we are running a paramagnetic calculation both blocks hold the same data. `solid_dmft` averaged both blocks for us to ensure a paramagnetic solution. To get an overview over all accessible attributes and functions of the Block Green's function, take a look [here](https://triqs.github.io/triqs/unstable/documentation/python_api/triqs.gf.block_gf.BlockGf.html).\n",
    "\n",
    "If we access one block we get back a view on the underlying Green's function in memory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31edfa5c-5cf6-482c-ade7-3718d35b46b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(Sigma_imp_iw['up_0'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cdcf20c-0129-4746-b221-654d6d011220",
   "metadata": {},
   "source": [
    "Every Green's function in TRIQS is defined on a mesh which can be accessed via `Sigma_imp_iw.mesh`. Also note that triqs implements many useful linear algebra operations that run directly in C++, like multiplication, addition, inversion etc. This leverages fast pre-compiled C++ functionality. We will later learn more how to use these operations."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d4b43f2-9a90-4e9a-9722-2d2c3ecd701b",
   "metadata": {},
   "source": [
    "### Plot of the self-energy Σ\n",
    "\n",
    "Let us now take a look at the impurity self-energy. It is the central output of DMFT calculation, which determines the corrections to the non-interacting DFT input Hamiltonian. Here we calculated Σ on the imaginary axis. Hence, in general features are not easily transferable to the real frequency axis. However, at $i \\omega_n \\rightarrow 0$ both the real and imaginary axis meet and we can directly interpret physical properties. Note that Σ itself like any other Green's function is complex-valued. Around $i \\omega_n \\approx 0$ we can read off quasiparticle properties. As long as the imaginary part of Σ goes to 0 for small Matsubara frequencies, the material has quasiparticle-like features. If the imaginary part of Σ diverges at 0, the material is Mott insulating. The slope around 0 determines the band renormalization, indicating the strength of correlations of the material. Let us plot the imaginary part of Σ via the TRIQS `oplot` plotting layer on top of matplotlib. This uses directly the `mesh` information to directly plot the Green's function object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69406e44-fca8-4af9-a107-f1f7777785aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot self energy\n",
    "fig, ax = plt.subplots(dpi=150, figsize=(6,3))\n",
    "\n",
    "ax.oplot(Sigma_imp_iw['up_0'].imag, '-o', c='C0', label='d$_{x^2-y^2}$')\n",
    "\n",
    "ax.set_ylabel(r\"Im $\\Sigma (i \\omega)$\")\n",
    "ax.set_xlim(0,20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ac5e6db-d41a-4d8e-aeda-63a11022c92f",
   "metadata": {},
   "source": [
    "We see that the self-energy is quite close to divergence around 0, but the first Matsubara goes back towards 0 indicating a metallic solution, but with strong renormalization. You should always make sure that the high frequency tail of the self-energy convergens smoothly to 0 on the imaginary axis, and to a constant value on the real axis.\n",
    "\n",
    "<i class=\"fa fa-gear fa-2x\" style=\"color: #186391\"></i>&nbsp; check the convergence of the real part of Σ . Does it look converged after increasing the number of MC cycles?\n",
    "\n",
    "<i class=\"fa fa-gear fa-2x\" style=\"color: #186391\"></i>&nbsp; now we want to converge the calculation better, by performing 4 additional iterations with better sampling, that are appended to the previous iterations. To do so, open the `dmft_config.toml` file in the main directory (not in the job folder) and change first the number of dmft iterations to 4, and then increase the number of MC measurements to $10 \\times 10^6$. Then just run solid_dmft again with the same command, but storing the output log file to `lco.dmft-hd.out`. Check that solid_dmft in the beginning notices that there is already a calculation folder with that name and continues the calculation starting from your current self-energy instead of starting a new one. After this is finished you can just rerun the above cells to load the new results, as the new results will be added to the existing h5 archive.\n",
    "\n",
    "<i class=\"fa fa-gear fa-2x\" style=\"color: #186391\"></i>&nbsp; while you wait for this to finish, plot the impurity Green's function in $\\tau$, this time the real part, which we loaded as `G_imp_tau`. You will see that it looks like a very noisy function. However, The true MC error is hard to estimate since the $\\tau$ mesh is very fine. To get a better feeling how accurate your measured Green's function is rebin on less $\\tau$ points by calling the member function `G_imp_tau['up_0'].rebinning_tau(new_n_tau=xx)` on it. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4800cad8-b7de-4a56-8b17-e9ca856dc269",
   "metadata": {},
   "source": [
    "## 4. Analytic continuation of the self-energy\n",
    "\n",
    "Since we are interested in the real frequency properties of the material we have to perform an analytic continuation. We do this using the TRIQS application [triqs.github.io/maxent](https://triqs.github.io/maxent/latest/). This is a stochastic method, which also performs quite well on noisy QMC data. We wrapped the functionality in [solid_dmft/maxent](https://triqs.github.io/solid_dmft/_ref/postprocessing.maxent_sigma.html) to automatically load a result from a DMFT calculation with all necessary input.\n",
    "\n",
    "Disclaimer: The analytic continuation of a self-energy is delicate and should be done with great care. Always check if the results are physically sound and avoid overinterpreting specific features of a MaxEnt-continued Green's function.\n",
    "\n",
    "We use now the `sigma_maxent` module of `solid_dmft`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a824c2f2-bad6-4528-b9dc-5a723b891772",
   "metadata": {},
   "outputs": [],
   "source": [
    "# use pcb maxent script to continue sigma\n",
    "Sigma_real_freq = sigma_maxent.main(external_path='b10-U3.6/lco_wannier.h5', \n",
    "                                     omega_min = -10, omega_max= 10, \n",
    "                                     maxent_error=0.03, iteration= None,\n",
    "                                     n_points_maxent=101,\n",
    "                                     continuator_type='inversion_sigmainf')[0][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23c19940-01dc-4817-b286-85ce51c2b20d",
   "metadata": {},
   "source": [
    "This takes just a few seconds in the current setup. You can play around using another, better converged iteration from the DMFT calculation from the h5 archive, or with the `n_points_maxent` (number of freq. points for the continuation). The function also stores the result consistently in the same h5 archive.\n",
    "\n",
    "Let us plot the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00299dbc-9246-44c5-afca-0b841cc91f6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, sharex=True, dpi=150, figsize=(10,4))\n",
    "\n",
    "# plot real part of the self-energy and subtract constant offset at 0.0 (Hartree shift)\n",
    "# Sigma_real_freq['up_0'](0.0) evaluates the self-energy at w=0.0\n",
    "ax[0].oplot((Sigma_real_freq['up_0'] - Sigma_real_freq['up_0'](0.0)[0,0]).real, '-', c='C0', label='Re')\n",
    "# plot imaginary part of the self-energy\n",
    "ax[1].oplot(-1*Sigma_real_freq['up_0'].imag, '-', c='C0', label='Im')\n",
    "\n",
    "ax[0].set_ylabel(r\"Re $\\Sigma (\\omega)$ (eV)\")\n",
    "ax[1].set_ylabel(r\"-Im $\\Sigma (\\omega)$ (eV)\")\n",
    "ax[-1].set_xlim(-2.5,2.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "227ebbbd-68cc-467a-b4c8-b4e935b6b3b7",
   "metadata": {},
   "source": [
    "<i class=\"fa fa-gear fa-2x\" style=\"color: #186391\"></i>&nbsp; how can the quasiparticle properties be read of the real frequency Σ? "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4ba2381-cafe-461e-81a2-23594105e64c",
   "metadata": {},
   "source": [
    "## 5. Plot Spectral Function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a34f90f6-423b-4ecf-94d8-91710942c942",
   "metadata": {},
   "source": [
    "The spectral function, which can be compared to photoemission experiments, is obtained by constructing first the lattice Green's functions: \n",
    "\n",
    "$$\\hat G(k,\\omega) = \\left[ \\omega + \\mu -\\hat{\\epsilon}(k)-\\left( \\hat\\Sigma(\\omega)^{imp}-\\hat\\Sigma^{dc} \\right) \\right]^{-1} $$\n",
    "\n",
    "using the DMFT analytically continued self-energy from above, the Wannier Hamiltonian $\\epsilon(k)$ and the obtained chemical potential. The spectral function is then calculated as: \n",
    "\n",
    "$$ A(k,\\omega) = - \\frac{1}{\\pi} \\text{Im} \\ G(k, \\omega) $$\n",
    "\n",
    "which is done by the function `get_dmft_bands` from solid_dmft postprocessing routines of `plot_correlated_bands.py`. By passing the option `with_Sigma='calc'` the self-energy from MaxEnt is loaded directly from the h5 file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb0ad5f8-9d6b-41a5-ad68-a94c42af82bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we have to specify the mesh of the spectral function and pass a new keywargs object sigma_dict\n",
    "w_mesh = {'window': [-4.0, 4.0], 'n_w': int(501)}\n",
    "sigma_dict = {'dmft_path': 'b10-U3.6/lco_wannier.h5',\n",
    "              'it': 'last_iter',\n",
    "              'spin': 'up', 'block' : 0,\n",
    "              'w_mesh': w_mesh}\n",
    "# plotting options\n",
    "plot_dict = {'colorscheme_bands': 'Greys',\n",
    "             'colorscheme_alatt': 'Spectral_r',\n",
    "             'colorscheme_kslice': 'Spectral_r'}\n",
    "\n",
    "tb_bands, alatt_k_w_path, freq_dict = pcb.get_dmft_bands(with_sigma='calc', **w90_dict, **tb_dict, **sigma_dict)\n",
    "\n",
    "fig, ax = plt.subplots(1,dpi=150,figsize=(7,4))\n",
    "pcb.plot_bands(fig, ax, alatt_k_w_path, tb_bands, freq_dict, n_orb=w90_dict['n_orb'], tb=True, alatt=True, **plot_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7278531f-58fd-4aa4-abe1-936a6ec03b35",
   "metadata": {},
   "source": [
    "The original W90 band is plotted in black and the spectral function is shown as intensity plot. The quasiparticle band is the strongly renormalized bright \"band\"-like structure around $\\omega=0$. The strong intensities at larger energies are Hubbard sattelite features forming the Hubbard bands. To visualize this effect better we now want to calculate the k-summed spectral function $A(\\omega)$, which we can compare with the DOS from W90.\n",
    "\n",
    "We first construct the full lattice Green's function on a homogeneous mesh:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21ba9b1d-0828-4d74-8e0f-d5ab43a3e6b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# triqs function to load a w90 hamiltonian as tight binding object\n",
    "TB = TB_from_wannier90(seed='lco' ,path='./data/mlwf/')\n",
    "\n",
    "# calculate epsilon(k) on a given mesh\n",
    "# since La2CuO4 is mostly 2D, we use a 2D mesh to speed things up\n",
    "k_dim = 7\n",
    "k_mesh = TB.get_kmesh(n_k = (k_dim, k_dim, 1))\n",
    "n_k = np.prod(k_mesh.dims)\n",
    "e_k = TB.fourier(k_mesh)\n",
    "\n",
    "# create a Mesh product of k and omega mesh:\n",
    "mesh_k_w = MeshProduct(k_mesh, Sigma_real_freq.mesh)\n",
    "\n",
    "# create a triqs Gf on a combined k and omega mesh:\n",
    "G_k_w = Gf(mesh = mesh_k_w, target_shape = [1,1])\n",
    "\n",
    "# create here first a Gf called G_loc_w for the k summation\n",
    "G_loc_w = None\n",
    "\n",
    "# now we loop over k-points\n",
    "for k in G_k_w.mesh.components[0]:\n",
    "    # use equation above to set at each k point the lattice Gf\n",
    "    G_k_w[k,:] << inverse(Omega + mu - e_k[k] - (Sigma_real_freq['up_0']-dc[0]['up']) )\n",
    "\n",
    "    # do the summation for G_loc_w"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c45006ff-fc86-4504-8bee-6af9731d8bfa",
   "metadata": {},
   "source": [
    "<i class=\"fa fa-gear fa-2x\" style=\"color: #186391\"></i>&nbsp; adapt the cell above to obtain the local real frequency Green's function `G_loc_w` defined in line 16. Think about how it should be normalized and also check [triqs.github.io/triqs/unstable/documentation/manual/triqs/gfs/py/block.html#operations](https://triqs.github.io/triqs/unstable/documentation/manual/triqs/gfs/py/block.html#operations) for how to perform basic operations on Gf objects.\n",
    "\n",
    "If everything went fine we can plot the result (hint: after it looks roughly correct increase `k_dim`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0e04dac-110a-410c-b190-85687427b1d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(dpi=150, figsize=(6,3))\n",
    "\n",
    "ax.oplot(-1/np.pi*G_loc_w.imag, label=r'A($\\omega$) U=3.6 eV')\n",
    "\n",
    "ax.plot(w90_dos[:,0] - fermi, w90_dos[:,1], label='w90 DOS')\n",
    "\n",
    "\n",
    "ax.set_ylabel(r'A($\\omega$)')\n",
    "ax.set_xlabel(r'$\\omega$')\n",
    "ax.set_xlim(-5,5)\n",
    "ax.set_ylim(0,1.1)\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa608e22-8d5b-424a-9c7e-19c888c9c801",
   "metadata": {},
   "source": [
    "### Analytically continue the impurity Green's function\n",
    "\n",
    "Instead of the self-energy we can also continue the impurity Green's function to the real frequency axis by running MaxEnt directly on `G_imp_tau` from the solution. Use the according function from solid_dmft (see below) to run MaxEnt on the impurity Green's function and plot the spectral function from it in comparison with that obtained from the lattice Green's function with the continued self-energy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5557dcaf-65dc-4931-92f6-6b50ae2eb1ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "import solid_dmft.postprocessing.maxent_gf_imp as gf_imp_maxent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c8434e5-e305-42d7-90d4-b47ab73c3031",
   "metadata": {},
   "outputs": [],
   "source": [
    "# solution:\n",
    "results = gf_imp_maxent.main(external_path='b10-U3.6/lco_wannier.h5', \n",
    "                             sum_spins=True, maxent_error=0.04,\n",
    "                             omega_min=-10, omega_max=10)[0]\n",
    "mesh = results['mesh']\n",
    "A_imp_w = results['Aimp_w_line_fit']['total_0']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "321254f9-7e9a-4614-b76e-2c17b9a15e48",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,dpi=150,figsize=(6,3))\n",
    "\n",
    "ax.oplot(-1/np.pi*G_loc_w.imag, label=r'A($\\omega$) $\\Sigma$ maxent')\n",
    "ax.plot(mesh, A_imp_w[0,0,:].real , label=r'A($\\omega$) $G_{imp}$ maxent')\n",
    "ax.plot(w90_dos[:,0]-fermi,w90_dos[:,1], label='w90 DOS')\n",
    "\n",
    "ax.set_ylabel(r'A($\\omega$)')\n",
    "ax.set_xlabel(r'$\\omega$')\n",
    "ax.set_xlim(-5,5)\n",
    "ax.set_ylim(0,1.1)\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d47dbe7-9e4e-4fce-84b6-49c3e70a1471",
   "metadata": {},
   "source": [
    "## 6. Plot the Fermi slice at $k_z = 0.0$ using the solid_dmft post-processing tool"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a31b147b-5af2-4bcd-b159-07beadb28f76",
   "metadata": {},
   "source": [
    "We can also look at in-plane Fermi slices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77e4c5f7-f711-45bf-8562-8cbdf003af2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "tb_kslice = {key: tb_dict[key] for key in list(tb_dict.keys()) if key.isupper()}\n",
    "kslice_update = {'bands_path': [('S\\'', 'G'),('G', 'R')], 'R': np.array([-0.25,0.25,0.25]), 'n_k': 50, 'kz': 0.0}\n",
    "tb_kslice.update(kslice_update)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32ded285-50e1-4d64-9a37-43c68ced8a32",
   "metadata": {},
   "outputs": [],
   "source": [
    "w90_dict = {'w90_path': './data/mlwf/',\n",
    "            'w90_seed': 'lco',\n",
    "            'n_orb': 1,\n",
    "            'mu_tb': fermi}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4195f701-57dc-48b3-8513-5d00573122cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we have to specify the mesh of the spectral function and pass a new keywargs object sigma_dict\n",
    "w_mesh = {'window': [-0.5, 0.5], 'n_w': int(1e6)}\n",
    "sigma_dict = {'dmft_path': 'b10-U3.6/lco_wannier.h5',\n",
    "              'it': 'last_iter',\n",
    "              'spin': 'up',\n",
    "              'block' : 0 ,\n",
    "              'w_mesh': w_mesh}\n",
    "\n",
    "tb_bands, alatt_k_w_path, freq_dict = pcb.get_dmft_bands(fermi_slice = True, with_sigma='calc', **w90_dict, **tb_kslice, **sigma_dict)\n",
    "\n",
    "fig, ax = plt.subplots(1,dpi=150,figsize=(7,4))\n",
    "pcb.plot_kslice(fig, ax, alatt_k_w_path, tb_bands, freq_dict, n_orb=w90_dict['n_orb'],\n",
    "                tb_dict=tb_kslice, tb=True, alatt=True, quarter=[0,1,2,3], **plot_dict)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc9b4192-3c86-4b07-8edd-f951b44c4bdd",
   "metadata": {},
   "source": [
    "We will now store our calculated self-energy on the real axis in the appropriate format be read by the FermiSee tool in the following tutorial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5152933d-c98e-41e0-bc56-086e4e0b585f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare self-energy for FermiSee\n",
    "w_mesh = [w.value for w in Sigma_real_freq.mesh]\n",
    "# store to FermiSee readable format\n",
    "with HDFArchive('sigma_fermisee.h5', 'w') as h5:\n",
    "    h5.create_group('self_energy')\n",
    "    h5['self_energy']['Sigma'] = Sigma_real_freq\n",
    "    h5['self_energy']['w_mesh'] = w_mesh\n",
    "    h5['self_energy']['n_w'] = len(w_mesh)\n",
    "    h5['self_energy']['n_orb'] = 1\n",
    "    h5['self_energy']['dc'] = dc[0]['up'][0,0]\n",
    "    h5['self_energy']['dmft_mu'] = mu\n",
    "    h5['self_energy']['orbital_order'] = [0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4310827a-6d4b-4eb4-a6e5-faa1baae4c91",
   "metadata": {},
   "source": [
    "## This concludes the first tutorial. You may now continue with tutorial 02-phase-diagram-exploration."
   ]
  }
 ],
 "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.10.10"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}