@@ -197,7 +197,7 @@
     "        smpl_wn  = {'F': self.IR_basis_set.smpl_wn_f,  'B': self.IR_basis_set.smpl_wn_b }[statistics]\n",
     "        return smpl_tau, smpl_wn\n",
     "\n",
-    "    \n",
+    "\n",
     "    def tau_to_wn(self, statistics, obj_tau):\n",
     "        \"\"\" Fourier transform from tau to iwn via IR basis \"\"\"\n",
     "        smpl_tau, smpl_wn = self.smpl_obj(statistics)\n",
@@ -216,7 +216,7 @@
     "        obj_tau = smpl_tau.evaluate(obj_l, axis=0).reshape((smpl_tau.tau.size, self.nk))\n",
     "        return obj_tau\n",
     "\n",
-    "    \n",
+    "\n",
     "    def k_to_r(self,obj_k):\n",
     "        \"\"\" Fourier transform from k-space to real space \"\"\"\n",
     "        obj_k = obj_k.reshape(-1, self.nk1, self.nk2)\n",
@@ -252,7 +252,7 @@
     "    def __init__(self, mesh, U, n, U_sfc_tol=1e-12, verbose=True):\n",
     "        \"\"\"\n",
     "        Solver class to calculate the TPSC method.\n",
-    "        After initializing the Solver by `solver = TPSCSolver(mesh, U, n, **kwargs)` it \n",
+    "        After initializing the Solver by `solver = TPSCSolver(mesh, U, n, **kwargs)` it\n",
     "        can be run by `solver.solve()`.\n",
     "        \"\"\"\n",
     "        ## set internal parameters for the solver\n",
@@ -261,35 +261,35 @@
     "        self.mesh = mesh\n",
     "        self.U_sfc_tol = U_sfc_tol\n",
     "        self.verbose = verbose\n",
-    "        \n",
+    "\n",
     "        ## set initial Green function and irreducible susceptibility\n",
     "        self.sigma = 0\n",
-    "        \n",
+    "\n",
     "        self.mu = 0\n",
     "        self.mu_calc()\n",
-    "        \n",
+    "\n",
     "        self.gkio_calc(self.mu)\n",
     "        self.grit_calc()\n",
     "        self.ckio_calc()\n",
-    "        \n",
+    "\n",
     "        # determine critical U_crit = 1/max(chi0) as an upper bound to U_sp\n",
     "        self.U_crit = 1/np.amax(self.ckio.real)\n",
-    "    \n",
-    "    \n",
+    "\n",
+    "\n",
     "    #%%%%%%%%%%% Solving instance\n",
     "    def solve(self):\n",
     "        \"\"\"\n",
     "        Determine spin and charge vertex self-consistently from sum rules and calculate self-energy.\n",
     "        \"\"\"\n",
     "        # determine spin vertex U_sp\n",
     "        self.spin_vertex_calc()\n",
-    "        \n",
+    "\n",
     "        # set double occupancy from Kanamori-Bruckner screening\n",
     "        self.docc_calc()\n",
-    " \n",
+    "\n",
     "        # determine charge vertex U_ch\n",
     "        self.charge_vertex_calc()\n",
-    "        \n",
+    "\n",
     "        # set spin and charge susceptibility\n",
     "        self.chi_spin   = self.RPA_term_calc( self.U_sp)\n",
     "        self.chi_charge = self.RPA_term_calc(-self.U_ch)\n",
@@ -299,7 +299,7 @@
     "        self.sigma_calc()\n",
     "        self.mu_calc()\n",
     "        self.gkio_calc(self.mu)\n",
-    "    \n",
+    "\n",
     "    #%%%%%%%%%%% Calculation steps for self.energy\n",
     "    def gkio_calc(self, mu):\n",
     "        \"\"\" Calculate Green function G(iw,k) \"\"\"\n",
@@ -313,8 +313,7 @@
     "\n",
     "    def ckio_calc(self):\n",
     "        \"\"\" Calculate irreducible susciptibility chi0(iv,q) \"\"\"\n",
-    "        # fix: from v1 to v2\n",
-    "        ckio = - self.grit * self.grit[::-1, :]\n",
+    "        ckio = self.grit * self.grit[::-1, :]\n",
     "\n",
     "        # Fourier transform\n",
     "        ckio = self.mesh.r_to_k(ckio)\n",
@@ -333,7 +332,7 @@
     "    def sigma_calc(self):\n",
     "        \"\"\" Calculate self-energy Sigma(iwn,k) \"\"\"\n",
     "        sigma = self.V * self.grit\n",
-    "    \n",
+    "\n",
     "        # Fourier transform\n",
     "        sigma = self.mesh.r_to_k(sigma)\n",
     "        self.sigma = self.mesh.tau_to_wn('F', sigma)\n",
@@ -343,8 +342,8 @@
     "    def RPA_term_calc(self, U):\n",
     "        \"\"\" Set RPA-like susceptibility \"\"\"\n",
     "        chi_RPA = self.ckio / (1 - U*self.ckio)\n",
-    "        return chi_RPA     \n",
-    "    \n",
+    "        return chi_RPA\n",
+    "\n",
     "    def chi_qtrace_calc(self, U):\n",
     "        \"\"\" Calculate (iv_m, q) trace of chi_RPA term \"\"\"\n",
     "        # chi_qtrace = sum_(m,q) chi(iv_m,q)\n",
@@ -353,47 +352,47 @@
     "        chi_trace_l  = self.mesh.IR_basis_set.smpl_wn_b.fit(chi_trace)\n",
     "        chi_trace = self.mesh.IR_basis_set.basis_b.u(0)@chi_trace_l\n",
     "        return chi_trace.real\n",
-    "    \n",
+    "\n",
     "    def docc_calc(self):\n",
     "        \"\"\" Calculate double occupancy from Kanamori-Bruckner type screening \"\"\"\n",
     "        self.docc = 0.25 * self.U_sp/self.U * self.n**2\n",
-    "    \n",
+    "\n",
     "    def spin_vertex_calc(self):\n",
     "        \"\"\" Determine U_sp self-consistently from local sum rule \"\"\"\n",
     "        # interval [U_a, U_b] for root finding\n",
     "        U_a = 0\n",
     "        U_b = np.floor(self.U_crit*100)/100\n",
-    "        \n",
+    "\n",
     "        chi_trace = self.chi_qtrace_calc\n",
     "        sfc_eq = lambda U_sp : 2*chi_trace(U_sp) - self.n + 0.5*(U_sp/self.U)*self.n**2\n",
     "\n",
-    "        if sfc_eq(U_b) > 0:        \n",
+    "        if sfc_eq(U_b) > 0:\n",
     "            self.U_sp = sc.optimize.brentq(sfc_eq, U_a, U_b, rtol = self.U_sfc_tol)\n",
     "            print(\"Setting U_sp = {}\".format(self.U_sp))\n",
     "        else:\n",
     "            print(\"Error\")\n",
     "            warn(\"System underwent phase transition, U^sp > U_crit = {}! U is too large or T too low for given doping.\".format(self.U_crit))\n",
-    "    \n",
+    "\n",
     "    def charge_vertex_calc(self):\n",
     "        \"\"\" Determine U_ch self-consistently from local sum rule \"\"\"\n",
     "        # interval [U_a, U_b] for root finding\n",
     "        U_a = 0\n",
     "        U_b = 100\n",
-    "        \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",
     "\n",
     "        self.U_ch = sc.optimize.brentq(sfc_eq, U_a, U_b, rtol = self.U_sfc_tol)\n",
     "\n",
-    "        \n",
+    "\n",
     "    #%%%%%%%%%%% Setting chemical potential mu\n",
     "    def calc_electron_density(self, mu):\n",
     "        \"\"\" Calculate chemical potential mu from Green function \"\"\"\n",
     "        self.gkio_calc(mu)\n",
     "        gio  = np.sum(self.gkio,axis=1)/self.mesh.nk\n",
     "        g_l  = self.mesh.IR_basis_set.smpl_wn_f.fit(gio)\n",
     "        g_tau0 = self.mesh.IR_basis_set.basis_f.u(0)@g_l\n",
-    "    \n",
+    "\n",
     "        n  = 1 + np.real(g_tau0)\n",
     "        n  = 2*n #for spin\n",
     "        return n\n",
@@ -455,12 +454,12 @@
     "# plot 2D k-dependence of lowest Matsubara frequency of e.g. green function\n",
     "plt.pcolormesh(2*mesh.k1.reshape(nk1,nk2), 2*mesh.k2.reshape(nk1,nk2), np.real(solver.gkio[mesh.iw0_f].reshape(mesh.nk1,mesh.nk2)), shading='auto')\n",
     "ax = plt.gca()\n",
-    "ax.set_xlabel('$k_x/\\pi$')\n",
+    "ax.set_xlabel(r'$k_x/\\pi$')\n",
     "ax.set_xlim([0,2])\n",
-    "ax.set_ylabel('$k_y/\\pi$')\n",
+    "ax.set_ylabel(r'$k_y/\\pi$')\n",
     "ax.set_ylim([0,2])\n",
     "ax.set_aspect('equal')\n",
-    "ax.set_title('Re $G(k,i\\omega_0)$')\n",
+    "ax.set_title(r'Re $G(k,i\\omega_0)$')\n",
     "plt.colorbar()\n",
     "plt.show()"
    ]
@@ -475,12 +474,12 @@
     "# plot 2D k-dependence of lowest Matsubara frequency of e.g. self-energy\n",
     "plt.pcolormesh(2*mesh.k1.reshape(nk1,nk2), 2*mesh.k2.reshape(nk1,nk2), np.imag(solver.sigma[mesh.iw0_f].reshape(mesh.nk1,mesh.nk2)), shading='auto')\n",
     "ax = plt.gca()\n",
-    "ax.set_xlabel('$k_x/\\pi$')\n",
+    "ax.set_xlabel(r'$k_x/\\pi$')\n",
     "ax.set_xlim([0,2])\n",
-    "ax.set_ylabel('$k_y/\\pi$')\n",
+    "ax.set_ylabel(r'$k_y/\\pi$')\n",
     "ax.set_ylim([0,2])\n",
     "ax.set_aspect('equal')\n",
-    "ax.set_title('Im $\\Sigma(k,i\\omega_0)$')\n",
+    "ax.set_title(r'Im $\\Sigma(k,i\\omega_0)$')\n",
     "plt.colorbar()\n",
     "plt.show()"
    ]
@@ -495,12 +494,12 @@
     "# plot 2D k-dependence of lowest Matsubara frequency of e.g. chi_spin\n",
     "plt.pcolormesh(2*mesh.k1.reshape(nk1,nk2), 2*mesh.k2.reshape(nk1,nk2), np.real(solver.chi_spin[mesh.iw0_b].reshape(mesh.nk1,mesh.nk2)), shading='auto')\n",
     "ax = plt.gca()\n",
-    "ax.set_xlabel('$k_x/\\pi$')\n",
+    "ax.set_xlabel(r'$k_x/\\pi$')\n",
     "ax.set_xlim([0,2])\n",
-    "ax.set_ylabel('$k_y/\\pi$')\n",
+    "ax.set_ylabel(r'$k_y/\\pi$')\n",
     "ax.set_ylim([0,2])\n",
     "ax.set_aspect('equal')\n",
-    "ax.set_title('$\\chi_{\\mathrm{sp}}(k,i\\nu_0)$')\n",
+    "ax.set_title(r'$\\chi_{\\mathrm{sp}}(k,i\\nu_0)$')\n",
     "plt.colorbar()\n",
     "plt.show()"
    ]
@@ -559,24 +558,24 @@
     "print(\"Start TPSC loop...\")\n",
     "for U_it, U in enumerate(U_array):\n",
     "    #print(\"Now: U = {:.1f}\".format(U))\n",
-    "    \n",
+    "\n",
     "    # TPSC solver\n",
     "    solver = TPSCSolver(mesh, U, n, verbose=False)\n",
     "    solver.solve()\n",
-    "    \n",
+    "\n",
     "    # save data for plotting\n",
     "    U_sp_array[U_it] = solver.U_sp\n",
     "    U_ch_array[U_it] = solver.U_ch\n",
     "print(\"Finished. Plotting now.\")\n",
     "\n",
     "\n",
     "#%%%%%%%%%%%%%%%% Plot results\n",
-    "plt.plot(U_array, U_ch_array, '-', label='$U_{\\mathrm{ch}}$')\n",
-    "plt.plot(U_array, U_sp_array, '-', label='$U_{\\mathrm{sp}}$')\n",
+    "plt.plot(U_array, U_ch_array, '-', label=r'$U_{\\mathrm{ch}}$')\n",
+    "plt.plot(U_array, U_sp_array, '-', label=r'$U_{\\mathrm{sp}}$')\n",
     "ax = plt.gca()\n",
-    "ax.set_xlabel('$U$', fontsize=12)\n",
+    "ax.set_xlabel(r'$U$', fontsize=12)\n",
     "ax.set_xlim([0,5])\n",
-    "ax.set_ylabel('Interaction', fontsize=12)\n",
+    "ax.set_ylabel(r'Interaction', fontsize=12)\n",
     "ax.set_ylim([0,20])\n",
     "ax.legend(frameon=False, fontsize=12)"
    ]