Merge pull request #16 from SpM-lab/fix/author-issues-jl

Fix the Julia notebooks reported in the author issues

Author: shinaoka

src/BgG_jl.ipynb

@@ -1094,18 +1094,22 @@
     "end\n",
     "\n",
     "\n",
+    "# RSCG stops when |r|^2 max_j|rho_j| < eps. The default eps = 1e-12 is too loose at eta = 0.01:\n",
+    "# the local DOS then deviates from exact diagonalization by up to 0.06. With eps = 1e-60 it agrees to 1e-11.\n",
+    "eps_rscg = 1e-60\n",
+    "\n",
     "ix = Nx ÷ 2\n",
     "iy = Ny ÷ 2\n",
     "ispin  =1\n",
     "\n",
     "i = (( iy-1)*Nx + ix-1)*2 + ispin\n",
     "j = i\n",
-    "Gij1 = greensfunctions(i,j,σ,Htsc) \n",
+    "Gij1 = greensfunctions(i,j,σ,Htsc; eps=eps_rscg)\n",
     "\n",
     "ispin  =2\n",
     "i = (( iy-1)*Nx + ix-1)*2 + ispin\n",
     "j = i\n",
-    "Gij2 = greensfunctions(i,j,σ,Htsc) \n",
+    "Gij2 = greensfunctions(i,j,σ,Htsc; eps=eps_rscg)\n",
     "plot(real.(σ),(-1/π)*imag.(Gij1 .+ Gij2),\n",
     "    xlabel = \"Energy [t]\",\n",
     "    ylabel = \"Local DOS\",)\n"
@@ -1132,12 +1136,12 @@
     "\n",
     "i = (( iy-1)*Nx + ix-1)*2 + ispin\n",
     "j = i\n",
-    "Gij1 = greensfunctions(i,j,σ,Htsc) \n",
+    "Gij1 = greensfunctions(i,j,σ,Htsc; eps=eps_rscg)\n",
     "\n",
     "ispin  =2\n",
     "i = (( iy-1)*Nx + ix-1)*2 + ispin\n",
     "j = i\n",
-    "Gij2 = greensfunctions(i,j,σ,Htsc) \n",
+    "Gij2 = greensfunctions(i,j,σ,Htsc; eps=eps_rscg)\n",
     "plot(real.(σ),(-1/π)*imag.(Gij1 .+ Gij2),\n",
     "    xlabel = \"Energy [t]\",\n",
     "    ylabel = \"Local DOS\",)"

src/DMFT_IPT_jl.ipynb

@@ -42,13 +42,16 @@
     "$\\mathcal{G}$ is the cavity Green function defined as\n",
     "\n",
     "$$\n",
-    "\\mathcal{G}(\\mathrm{i}\\nu) = \\frac{1}{\\mathrm{i}\\nu - \\mu - \\Delta(\\mathrm{i}\\nu)},\n",
+    "\\mathcal{G}(\\mathrm{i}\\nu) = \\frac{1}{\\mathrm{i}\\nu + \\mu - \\Delta(\\mathrm{i}\\nu)},\n",
     "$$\n",
     "\n",
     "$\\Delta$ is the hybrization function,\n",
     "$\\tau$ is the imaginary time, $\\beta$ is the inverse temperature.\n",
     "Here, we compute the self-energy within the interative perturbation theory (IPT) [1].\n",
     "\n",
+    "At half filling, $\\mu = U/2$ cancels the Hartree term $\\Sigma_\\mathrm{H} = U n_\\sigma = U/2$ of the self-energy. The code therefore contains neither $\\mu$ nor $\\Sigma_\\mathrm{H}$: its $\\Sigma$ is the second-order term, the denominator of $G_{\\mathrm{loc}}$ reads $\\mathrm{i}\\nu-\\omega-\\Sigma(\\mathrm{i}\\nu)$, and $\\mathcal{G}^{-1} = G_{\\mathrm{loc}}^{-1} + \\Sigma$ gives $\\mathrm{i}\\nu - \\Delta(\\mathrm{i}\\nu)$, the cavity function shifted by $\\Sigma_\\mathrm{H}$ that enters IPT.\n",
+    "The semicircular density of states above is that of the Bethe lattice with hopping $t_\\star = 1$, for which these equations are equivalent to $\\mathcal{G}^{-1}(\\mathrm{i}\\nu) = \\mathrm{i}\\nu - t_\\star^2 G_{\\mathrm{loc}}(\\mathrm{i}\\nu)$. Particle-hole symmetry makes $\\mathcal{G}(\\mathrm{i}\\nu)$ purely imaginary (the code keeps only its imaginary part), so that $\\mathcal{G}(\\beta-\\tau) = \\mathcal{G}(\\tau)$ and $\\Sigma(\\tau) = U^2\\mathcal{G}^3(\\tau)$. The equations are therefore the same as in the [Python version](DMFT_IPT_py.ipynb) of this tutorial.\n",
+    "\n",
     "### Renormalization factors\n",
     "\n",
     "The Mott transition can be detected by monitoring the renormalization factor\n",
@@ -161,7 +164,7 @@
     "        smpl_matsu = MatsubaraSampling(basis)\n",
     "        smpl_tau = TauSampling(basis)\n",
     "\n",
-    "        new(U, basis, SparseIR.beta(basis), rho_omega, omega_range, quad_rule, smpl_matsu, smpl_tau)\n",
+    "        new(U, basis, SparseIR.beta(basis), rho, omega_range, quad_rule, smpl_matsu, smpl_tau)\n",
     "    end\n",
     "end\n",
     ";"

src/FLEX_jl.ipynb

@@ -188,13 +188,14 @@
     "end\n",
     "\n",
     "function k_to_r(mesh::Mesh, 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_r = fft(obj_k, [2, 3])\n",
     "    return obj_r\n",
     "end\n",
     "\n",
     "function r_to_k(mesh::Mesh, 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 ifft 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_k = ifft(obj_r, [2, 3]) / mesh.nk\n",
     "    return obj_k\n",
     "end\n",