Merge pull request #15 from SpM-lab/fix/author-issues-py

Fix the Python notebooks reported in the author issues

Author: shinaoka

pyproject.toml

@@ -7,7 +7,7 @@ dependencies = [
     "matplotlib",
     "ghp-import",
     "jupyterlab",
-    "sparse-ir>=2.0.0a9",
+    "sparse-ir>=2.1.5",
     "admmsolver",
-    "pylibsparseir>=0.7.5",
+    "pylibsparseir>=0.10.0",
 ]

src/FLEX_py.ipynb

@@ -194,14 +194,15 @@
     "\n",
     "    \n",
     "    def k_to_r(self,obj_k):\n",
-    "        \"\"\" Fourier transform from k-space to real space \"\"\"\n",
+    "        \"\"\" Fourier transform from k-space to real space, A(r) = sum_k exp(-ik.r) A(k), without a factor 1/nk \"\"\"\n",
     "        obj_k = obj_k.reshape(-1, self.nk1, self.nk2)\n",
     "        obj_r = np.fft.fftn(obj_k,axes=(1,2))\n",
     "        obj_r = obj_r.reshape(-1, self.nk)\n",
     "        return obj_r\n",
     "\n",
     "    def r_to_k(self,obj_r):\n",
-    "        \"\"\" Fourier transform from real space to k-space \"\"\"\n",
+    "        \"\"\" Fourier transform from real space to k-space, A(k) = (1/nk^2) sum_r exp(ik.r) A(r) (1/nk from ifftn and 1/nk here),\n",
+    "        so that the product of two k_to_r arrays A(r)B(r) gives the convolution (1/nk) sum_q A(q)B(k-q) \"\"\"\n",
     "        obj_r = obj_r.reshape(-1, self.nk1, self.nk2)\n",
     "        obj_k = np.fft.ifftn(obj_r,axes=(1,2))/self.nk\n",
     "        obj_k = obj_k.reshape(-1, self.nk)\n",

src/GW_py.ipynb

@@ -93,7 +93,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Above we have successfully derived $G^\\mathrm{F}(\\bar{\\tau}_k^\\mathrm{F})$ from $G^\\mathrm{F}(\\mathrm{i}\\bar{\\nu}_k^\\mathrm{F})$ by first obtaining the Green's function's basis representation $G_l^\\mathrm{F}$ and secondly evaluating $G_l^\\mathrm{F}$ on the $\\tau$-sampling points, which are given by the computations of sparse_ir. Throughout, $\\bar{\\tau}_k^\\mathrm{F}$ and $\\mathrm{i}\\bar{\\nu}_k^\\mathrm{F}$ ($\\bar{\\tau}_k^\\mathrm{B}$ and $\\mathrm{i}\\bar{\\nu}_k^\\mathrm{B}$) denote the fermionic (bosonic) sampling points in imaginary time and Matsubara frequency, where $\\nu^\\mathrm{F} = n\\pi/\\beta$ with odd $n$ and $\\nu^\\mathrm{B} = n\\pi/\\beta$ with even $n$. Here, the self-consistent second order Green's function theory (GF2) and the self-consistent GW method are are based this very mechanism (and the so called Hedin equations).[5][7]\n",
+    "Above, we have computed the non-interacting Green's function $G^\\mathrm{F}(\\mathrm{i}\\bar{\\nu}_k^\\mathrm{F})$ on the Matsubara sampling points. Below, we derive $G^\\mathrm{F}(\\bar{\\tau}_k^\\mathrm{F})$ from $G^\\mathrm{F}(\\mathrm{i}\\bar{\\nu}_k^\\mathrm{F})$ by first obtaining the Green's function's basis representation $G_l^\\mathrm{F}$ and secondly evaluating $G_l^\\mathrm{F}$ on the $\\tau$-sampling points, which are given by the computations of sparse_ir. Throughout, $\\bar{\\tau}_k^\\mathrm{F}$ and $\\mathrm{i}\\bar{\\nu}_k^\\mathrm{F}$ ($\\bar{\\tau}_k^\\mathrm{B}$ and $\\mathrm{i}\\bar{\\nu}_k^\\mathrm{B}$) denote the fermionic (bosonic) sampling points in imaginary time and Matsubara frequency, where $\\nu^\\mathrm{F} = n\\pi/\\beta$ with odd $n$ and $\\nu^\\mathrm{B} = n\\pi/\\beta$ with even $n$. Here, the self-consistent second order Green's function theory (GF2) and the self-consistent GW method are are based this very mechanism (and the so called Hedin equations).[5][7]\n",
     "\n",
     "The following image should give an overview on their iterative schemes[5]:\n",
     "\n",
@@ -394,13 +394,13 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "After changing back into fermionic statistics the next quantity that is being dealt with is the so-called self energy $\\Sigma$. When we first introduced the Dyson equation in the first chapter, we saw $V$ as some sort of pertubation. The self-energy is very much alike as it describes the correlation effects of a many-body system. Here we will calculate $\\tilde{\\Sigma}(\\bar{\\tau}_k^\\mathrm{F})$ using $\\Sigma^{GW}$ so that\n",
+    "After changing back into fermionic statistics the next quantity that is being dealt with is the so-called self energy $\\Sigma$. It enters the Dyson equation below and describes the correlation effects of a many-body system. Here we will calculate $\\tilde{\\Sigma}(\\bar{\\tau}_k^\\mathrm{F})$ using $\\Sigma^{GW}$ so that\n",
     "\n",
-    "$${\\Sigma}(\\bar{\\tau}_k^\\mathrm{F})=-G(\\bar{\\tau}_k^\\mathrm{F})*{W}(\\bar{\\tau}_k^\\mathrm{F}).$$\n",
+    "$${\\Sigma}(\\bar{\\tau}_k^\\mathrm{F})=G(\\bar{\\tau}_k^\\mathrm{F})*{W}(\\bar{\\tau}_k^\\mathrm{F}).$$\n",
     "\n",
     "This can again be rewritten into its equivalent form\n",
     "\n",
-    "$$\\Sigma=-G(U+UP(U+UP(...)))= -(GU+GUPU+GUPUPU+...)$$"
+    "$$\\Sigma=G(U+UP(U+UP(...)))= GU+GUPU+GUPUPU+...$$"
    ]
   },
   {

src/TPSC_py.ipynb

@@ -57,7 +57,7 @@
     "\n",
     "Afterwards we can calculate the double occupancy $\\langle n_{\\uparrow}n_{\\downarrow}\\rangle = \\frac{U_{\\mathrm{sp}}}{4U} n^2$ and then perform a similar root finding for $U_{\\mathrm{ch}}$ from the function \n",
     "\n",
-    "$$ g(U_{\\mathrm{ch}}) = 2 \\frac{T}{N_{\\boldsymbol{k}}} \\sum_{\\nu^\\mathrm{B},\\boldsymbol{q}} \\chi_{\\mathrm{ch}}(U_{\\mathrm{ch}}) - n - 2\\langle n_{\\uparrow}n_{\\downarrow}\\rangle^2 + n^2\\;. $$\n",
+    "$$ g(U_{\\mathrm{ch}}) = 2 \\frac{T}{N_{\\boldsymbol{k}}} \\sum_{\\nu^\\mathrm{B},\\boldsymbol{q}} \\chi_{\\mathrm{ch}}(U_{\\mathrm{ch}}) - n - 2\\langle n_{\\uparrow}n_{\\downarrow}\\rangle + n^2\\;. $$\n",
     "\n",
     "In TPSC, a self-energy $\\Sigma$ can be derived that is calculated from the interaction {cite:p}`Moukouri2000`\n",
     "\n",
@@ -218,14 +218,15 @@
     "\n",
     "\n",
     "    def k_to_r(self,obj_k):\n",
-    "        \"\"\" Fourier transform from k-space to real space \"\"\"\n",
+    "        \"\"\" Fourier transform from k-space to real space, A(r) = sum_k exp(-ik.r) A(k), without a factor 1/nk \"\"\"\n",
     "        obj_k = obj_k.reshape(-1, self.nk1, self.nk2)\n",
     "        obj_r = np.fft.fftn(obj_k,axes=(1,2))\n",
     "        obj_r = obj_r.reshape(-1, self.nk)\n",
     "        return obj_r\n",
     "\n",
     "    def r_to_k(self,obj_r):\n",
-    "        \"\"\" Fourier transform from real space to k-space \"\"\"\n",
+    "        \"\"\" Fourier transform from real space to k-space, A(k) = (1/nk^2) sum_r exp(ik.r) A(r) (1/nk from ifftn and 1/nk here),\n",
+    "        so that the product of two k_to_r arrays A(r)B(r) gives the convolution (1/nk) sum_q A(q)B(k-q) \"\"\"\n",
     "        obj_r = obj_r.reshape(-1, self.nk1, self.nk2)\n",
     "        obj_k = np.fft.ifftn(obj_r,axes=(1,2))/self.nk\n",
     "        obj_k = obj_k.reshape(-1, self.nk)\n",
@@ -380,7 +381,7 @@
     "        U_b = 100\n",
     "\n",
     "        chi_trace = self.chi_qtrace_calc\n",
-    "        sfc_eq = lambda U_ch : 2*chi_trace(-U_ch) - self.n + (1 - 2*self.docc)*self.n**2\n",
+    "        sfc_eq = lambda U_ch : 2*chi_trace(-U_ch) - self.n - 2*self.docc + self.n**2\n",
     "\n",
     "        self.U_ch = sc.optimize.brentq(sfc_eq, U_a, U_b, rtol = self.U_sfc_tol)\n",
     "\n",

src/analytic_continuation_py.ipynb

@@ -116,10 +116,10 @@
    "outputs": [],
    "source": [
     "def semicirc_dos(w):\n",
-    "    return 2/np.pi * np.sqrt((np.abs(w) < wmax) * (1 - np.square(w/wmax)))\n",
+    "    return 2/(np.pi*wmax) * np.sqrt((np.abs(w) < wmax) * (1 - np.square(w/wmax)))\n",
     "\n",
     "def insulator_dos(w):\n",
-    "    return semicirc_dos(8*w/wmax - 4) + semicirc_dos(8*w/wmax + 4)\n",
+    "    return 4/wmax * (semicirc_dos(8*w/wmax - 4) + semicirc_dos(8*w/wmax + 4))\n",
     "\n",
     "# For testing, compute exact coefficients g_l for two models\n",
     "rho1_l = basis.v.overlap(semicirc_dos)\n",
@@ -232,7 +232,7 @@
     "As a simple example, let us consider **Ridge regression** in the above problem.  We again expand the spectral function in the right singular functions with $L'=L$, but include a regularization term:\n",
     "\n",
     "$$\n",
-    "    f_\\mathrm{reg}[\\rho] = \\alpha \\sum_{l=0}^{L-1} | \\rho_l |^2,\n",
+    "    f_\\mathrm{reg}[\\rho] = \\alpha^2 \\sum_{l=0}^{L-1} | \\rho_l |^2,\n",
     "$$\n",
     "\n",
     "where $\\alpha$ is a hyperparameter (ideally tuned to the noise level).  This term prevents $\\rho_l$ becoming too large due to noise amplification. The regularized least squares problem (2) then amounts to:\n",

src/eliashberg_holstein_py.ipynb

@@ -399,12 +399,13 @@
    },
    "outputs": [],
    "source": [
-    "def add_noise(arr, noise):\n",
+    "def add_noise(arr, noise, seed=0):\n",
     "    \"\"\"\n",
-    "    Add Gaussian noise to an array\n",
+    "    Add Gaussian noise to an array (seeded, so that the results are reproducible)\n",
     "    \"\"\"\n",
-    "    arr += noise*np.random.randn(*arr.shape)\n",
-    "    arr += noise*1j*np.random.randn(*arr.shape)\n",
+    "    rng = np.random.default_rng(seed)\n",
+    "    arr += noise*rng.standard_normal(arr.shape)\n",
+    "    arr += noise*1j*rng.standard_normal(arr.shape)\n",
     "    return arr"
    ]
   },

src/liechtenstein_py.ipynb

@@ -201,7 +201,7 @@
    "outputs": [],
    "source": [
     "Lambda = 2*max(8*t, b_eff) * beta\n",
-    "basis = sparse_ir.FiniteTempBasis(statistics='F', beta=beta, wmax=Lambda, eps=1e-7)\n",
+    "basis = sparse_ir.FiniteTempBasis(statistics='F', beta=beta, wmax=Lambda/beta, eps=1e-7)\n",
     "l = liechtenstein(basis, t, beta, nk, b_eff)"
    ]
   },

src/second_order_perturbation_py.ipynb

@@ -41,7 +41,7 @@
     "The corresponding non-interacting Green's function is given by\n",
     "\n",
     "$$\n",
-    "    G(\\mathrm{i}\\nu, \\boldsymbol{k}) = \\frac{1}{\\mathrm{i}\\nu - \\epsilon(\\boldsymbol{k}) - \\mu},\n",
+    "    G(\\mathrm{i}\\nu, \\boldsymbol{k}) = \\frac{1}{\\mathrm{i}\\nu - \\epsilon(\\boldsymbol{k}) + \\mu},\n",
     "$$\n",
     "\n",
     "where $\\nu = n\\pi/\\beta$ with odd $n$ is a fermionic Matsubara frequency.\n",