{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"editable": true,
"nbsphinx": "hidden",
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"\n",
"import warnings \n",
"warnings.filterwarnings(\"ignore\") #ignore some matplotlib warnings\n",
"\n",
"from triqs.plot.mpl_interface import plt, oplot, oplotr, oploti\n",
"plt.rcParams[\"figure.figsize\"] = (8,5) # set default size for all figures"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"# Analytic continuation of Hubbard dimer"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [],
"vscode": {
"languageId": "latex"
}
},
"source": [
"In this notebook, we illustrate how to conduct analytic continuation using `adapol`. We took the strongly correlated hubbard dimer model as an example. On the one hand, it serves as a nontrivial test case which features several closely spaced peaks that could be very challenging. On the other hand, its modest size allows for exact diagonalization, enabling the calculation of ground truth results for rigorous benchmarking. This example was also used as benchmark tests in [this](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.104.165111) and [this](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.107.075151) paper.\n",
"\n",
"Let us consider the following Hamiltonian:\n",
"$$\n",
"\\hat H = \\hat H_0 + \\hat H_1,\n",
"$$\n",
"where\n",
"$$\n",
"H_0 = -t \\sum_{\\sigma \\in \\{\\uparrow, \\downarrow\\}} \\left( \\hat c_{0\\sigma}^\\dagger \\hat c_{1\\sigma} + \\hat c_{1\\sigma}^\\dagger \\hat c_{0\\sigma} \\right) \n",
"-\\mu \\left( \\hat n_{0\\uparrow} + \\hat n_{0\\downarrow} + \\hat n_{1\\uparrow} + \\hat n_{1\\downarrow} \\right),\n",
"$$\n",
"$$\n",
"\\begin{aligned}\n",
" H_1 =& (U+U_a) \\hat n_{0\\uparrow} \\hat n_{0\\downarrow} + (U-U_a) \\hat n_{1\\uparrow} \\hat n_{1\\downarrow} \\\\\n",
"& - \\frac U 2 \\left( \\hat n_{0\\uparrow} + \\hat n_{0\\downarrow} + \\hat n_{1\\uparrow} + \\hat n_{1\\downarrow} \\right) + h \\left( \\hat n_{0\\uparrow} - \\hat n_{0\\downarrow} + \\hat n_{1\\uparrow} - \\hat n_{1\\downarrow} \\right)\\\\\n",
"&+ \\mu_a \\left( \\hat n_{0\\uparrow} + \\hat n_{0\\downarrow} - \\hat n_{1\\uparrow} - \\hat n_{1\\downarrow} \\right)+ h_a \\left( \\hat n_{0\\uparrow} - \\hat n_{0\\downarrow} - \\hat n_{1\\uparrow} + \\hat n_{1\\downarrow} \\right),\n",
"\\end{aligned}\n",
"$$\n",
"Here $\\hat n_{i\\sigma} = \\hat c_{i\\sigma}^\\dagger \\hat c_{i\\sigma}$ is the number operator.\n",
"\n",
"The parameters are taken as follows:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"from triqs.operators import c, c_dag,n\n",
"from triqs.atom_diag import AtomDiag, atomic_g_w, atomic_g_iw\n",
"from itertools import product\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"\n",
"t = 1.0\n",
"mu = 0.7\n",
"h = 0.3\n",
"U_a = 0.5\n",
"mu_a = 0.2\n",
"h_a = 0.03\n",
"beta = 25.0\n",
"U = 5.0"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"Using `triqs`, we can construct the Hamiltonian and obtain the Green's function and spectral function through exact diagonalization:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"f_ops = [(sn,on) for sn, on in product(('up','down'),range(2))]\n",
"H_0 = - t * (c_dag('up', 0) * c('up', 1) + c_dag('up', 1) * c('up', 0) )\n",
"H_0 += - t * (c_dag('down', 0) * c('down', 1) + c_dag('down', 1) * c('down', 0) )\n",
"H_0 += - mu * (n('up', 0) + n('down', 0)) - mu * (n('up', 1) + n('down', 1))\n",
"\n",
"H_1 = (U + U_a) * (n('up', 0) * n('down', 0) ) + (U - U_a) * (n('up', 1) * n('down', 1) )\n",
"H_1 -= (U/2.0) * (n('up', 0) + n('down', 0) + n('up', 1) + n('down', 1) )\n",
"H_1 += h * (n('up', 0) - n('down', 0) + n('up', 1) - n('down', 1) )\n",
"H_1 += mu_a * (n('up', 0) + n('down', 0) - n('up', 1) - n('down', 1) )\n",
"H_1 += h_a * (n('up', 0) - n('down', 0) - n('up', 1) + n('down', 1) )\n",
"\n",
"H = H_0 + H_1\n",
"ad = AtomDiag(H, f_ops)\n",
"\n",
"gf_struct = [('down',2),('up',2)] # fix the bug in the old version by changing orb_names to len(orb_names).\n",
"G_w = atomic_g_w(ad, beta, gf_struct, (-7, 7), 4000, 0.01)\n",
"Tr_G_w = G_w['up'][0,0] + G_w['down'][0,0] + G_w['up'][1,1] + G_w['down'][1,1]\n",
"wmesh = np.array([w.value for w in Tr_G_w.mesh])\n",
"plt.plot(wmesh, -np.imag(Tr_G_w.data)/np.pi)\n",
"plt.title(\"Spectral function\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"Let us also construct the Matsubara Green's function data:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"G_iw = atomic_g_iw(ad, beta, gf_struct, 500)\n",
"TrG_iw = G_iw['up'][0,0] + G_iw['down'][0,0] + G_iw['up'][1,1] + G_iw['down'][1,1]\n",
"iwmesh = np.array([iw.value.value for iw in TrG_iw.mesh])"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"Let us first show how to conduct analytic continuation through the `anacont` function:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from adapol import approx_freq_aaa\n",
"from adapol.sop import SumOfSimplePoles\n",
"\n",
"poles, residues, err = approx_freq_aaa(TrG_iw.data, iwmesh, max_n_poles=12)\n",
"sop = SumOfSimplePoles(poles=poles, residues=residues)\n",
"Spec_approx = -np.imag(sop(wmesh+0.01*1j))/np.pi\n",
"\n",
"plt.plot(wmesh, Spec_approx) \n",
"plt.plot(wmesh, -np.imag(Tr_G_w.data)/np.pi, \"--\",2)\n",
"plt.legend([\"Analytic continuation\", \"Ground truth\"])\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"Next let's see how to conduct analytic continuation through the triqs interface:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from adapol.triqs import approx_gf_imfreq_aaa\n",
"\n",
"poles, residues, err = approx_gf_imfreq_aaa(TrG_iw, max_n_poles=12)\n",
"sop2 = SumOfSimplePoles(poles=poles, residues=residues)\n",
"Spec_approx2 = -np.imag(sop2(wmesh+0.01*1j))/np.pi\n",
"\n",
"plt.plot(wmesh, Spec_approx2) \n",
"plt.plot(wmesh, -np.imag(Tr_G_w.data)/np.pi, \"--\",2)\n",
"plt.legend([\"Analytic continuation\", \"Ground truth\"])\n",
"plt.show()"
]
}
],
"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.12"
}
},
"nbformat": 4,
"nbformat_minor": 4
}