fix: actions and pytest strictness (#37)

* fix: make get_h_qp helper only return the matrix

* docs: reflect change to h_qp helper

* feat: helper functions to construct triqs green's functions

* chore: added green's function to example

* chore: ignore version file

* fix: update actions, use dependabot

* fix: typo in gcc-version

* fix: ignore triqs errors in pytest

* check if ignore warnings works on github

* only filter triqs warnings

* fix: only make warnings errors in risb, ignore third party

.github/actions/setup-triqs/action.yml

@@ -25,7 +25,7 @@ runs:
 
   steps:
   - name: Setup Python ${{ inputs.python-version }}
-    uses: actions/setup-python@v4
+    uses: actions/setup-python@v5
     with:
       python-version: ${{ inputs.python-version }}
       allow-prereleases: true
@@ -106,7 +106,7 @@ runs:
       
   - name: Cache build TRIQS
     id: build-triqs
-    uses: actions/cache@v3
+    uses: actions/cache@v4
     with:
       path: triqs
       key: triqs-${{ inputs.os }}-${{ inputs.python-version }} ${{ inputs.cc }}-${{ inputs.llvm }}-${{ inputs.gcc-version }}

.github/dependabot.yml

@@ -0,0 +1,11 @@
+version: 2
+updates:
+  # Maintain dependencies for GitHub Actions
+  - package-ecosystem: "github-actions"
+    directory: "/"
+    schedule:
+      interval: "weekly"
+    groups:
+      actions:
+        patterns:
+          - "*"

.github/workflows/ci.yml

@@ -19,7 +19,8 @@ jobs:
       matrix:
         include:
           - {os: "ubuntu-22.04", python-version: "3.10", cc: "gcc", cxx: "g++", llvm: "15", gcc-version: "12"}
-          #- {os: "ubuntu-22.04", python-version: "3.11", cc: "gcc", cxx: "g++", llvm: "15", gccversion: "12"}
+          - {os: "ubuntu-22.04", python-version: "3.11", cc: "gcc", cxx: "g++", llvm: "15", gcc-version: "12"}
+          - {os: "ubuntu-22.04", python-version: "3.12", cc: "gcc", cxx: "g++", llvm: "15", gcc-version: "12"}
           #- {os: "ubuntu-22.04", python-version: "3.10", cc: "clang", cxx: "clang++", llvm: "15", gcc-version: "12"}
           #- {os: "ubuntu-22.04", python-version: "3.11", cc: "clang", cxx: "clang++", llvm: "15", gcc-version: "12"}
 

.gitignore

@@ -16,4 +16,6 @@ _build
 .github.auth
 
 # Visual Studio Code
-.vscode
+.vscode
+
+src/risb/version.py

docs/tutorials/self-consistent.md

@@ -230,17 +230,18 @@ input initial guesses to the solution to the self-consistent loop. Often a
 reasonable initial guess is to choose $\mathcal{R}$ as the identity and 
 $\lambda$ as the zero matrix.
 
-There is a helper function that constructs $\hat{H}^{\mathrm{qp}}$ and 
-returns its eigenvalues and eigenvectors at each $k$-point on the finite 
-grid.
+There is a helper function that constructs $\hat{H}^{\mathrm{qp}}$, and 
+it is simple to get its eigenvalues and eigenvectors at each $k$-point on the 
+finite grid.
 
 ```python
 from risb import helpers
 
 energies_qp = dict()
 bloch_vector_qp = dict()
 for bl, bl_size in gf_struct
-    energies_qp[bl], bloch_vector_qp[bl] = helpers.get_h_qp(R[bl], Lambda[bl], h0_kin_k[bl])
+    h_qp = helpers.get_h_qp(R[bl], Lambda[bl], h0_kin_k[bl])
+    energies_qp[bl], bloch_vector_qp[bl] = np.linalg.eigh(h_qp)
 ```
 
 ## D: Setup k-integrator function

examples/hubbard.py

@@ -4,10 +4,12 @@
 from triqs.operators import Operator, n
 from triqs.operators.util.op_struct import set_operator_structure
 from triqs.operators.util.observables import S2_op, N_op
+from triqs.gf import BlockGf, MeshImFreq, MeshProduct
 
 from risb import LatticeSolver
 from risb.kweight import SmearingKWeight
 from risb.embedding import EmbeddingAtomDiag
+from risb.helpers_triqs import get_sigma_w, get_g_qp_k_w, get_g_k_w, get_g_w_loc
 
 import matplotlib.pyplot as plt
 
@@ -82,7 +84,40 @@ def force_paramagnetic(A):
     total_number.append( embedding.overlap(total_number_Op) )
     total_spin.append( embedding.overlap(total_spin_Op) )
     Z.append(S.Z[0]['up'][0,0])
+
+    if np.abs(U - 3.5) < 1e-10:
+        # Some different ways to construct some Green's functions
+        mu = kweight.mu
+
+        # Non-interacting lattice Green's function
+        iw_mesh = MeshImFreq(beta=beta, S='Fermion', n_max=64)
+        k_iw_mesh = MeshProduct(mk, iw_mesh)
+        G0_k_iw = BlockGf(mesh=k_iw_mesh, gf_struct=gf_struct)
+        for bl, gf in G0_k_iw:
+            e_k = tbl.fourier(mk)
+            for k, iw in gf.mesh:
+                gf[k,iw] = 1 / (iw - e_k[k] + mu)
+
+        # Quasiparticle lattice Green's function, local self-energy, lattice Green's function
+        G_qp_k_iw = get_g_qp_k_w(gf_struct=gf_struct, mesh=k_iw_mesh, h0_kin_k=S.h0_kin_k, Lambda=S.Lambda[0], R=S.R[0], mu=mu)
+        Sigma_iw = get_sigma_w(mesh=iw_mesh, gf_struct=gf_struct, Lambda=S.Lambda[0], R=S.R[0], mu=mu)
+        G_k_iw = get_g_k_w(g0_k_w=G0_k_iw, sigma_w=Sigma_iw)
+        G_k_iw2 = get_g_k_w(g_qp_k_w=G_qp_k_iw, R=S.R[0])
 
+        # Local Green's functions integrated over k
+        G0_iw_loc = get_g_w_loc(G0_k_iw)
+        G_qp_iw_loc = get_g_w_loc(G_qp_k_iw)
+        G_iw_loc = get_g_w_loc(G_k_iw)
+        G_iw_loc2 = get_g_w_loc(G_k_iw2)
+
+        # Filling of physical electron scales with Z
+        print("G0:", G0_iw_loc.total_density().real)
+        print("G_qp:", G_qp_iw_loc.total_density().real)
+        print("G:", G_iw_loc.total_density().real)
+        print("G2:", G_iw_loc2.total_density().real)
+        print("Z:", S.Z[0]['up'][0,0])
+        print()
+
 fig, axis = plt.subplots(1,2)
 axis[0].plot(U_arr, Z, '-ok')
 axis[0].plot(U_arr, -0.5 + 0.5 * np.sqrt( 1 + 4*np.array(total_spin) ), '-ok')

pyproject.toml

@@ -55,8 +55,7 @@ minversion = "7.0"
 addopts = ["-ra", "--showlocals", "--strict-markers", "--strict-config"]
 xfail_strict = true
 filterwarnings = [
-    "ignore:invalid escape sequence:DeprecationWarning",
-    "error",
+  "error:::risb",
 ]
 log_cli_level = "info"
 testpaths = [

src/risb/helpers.py

@@ -248,26 +248,24 @@ def get_h_qp(R : np.ndarray,
              h0_kin_k : np.ndarray, 
              mu : float = 0) -> tuple[ np.ndarray, np.ndarray ]:
     """
-    Construct eigenvalues and eigenvectors of the quasiparticle Hamiltonian.
+    Construct the quasiparticle Hamiltonian.
     
     Parameters
     ----------
     R : numpy.ndarray
-        Renormalization matrix) from electrons to quasiparticles.
+        Renormalization matrix) from electrons to quasiparticles
     Lambda : numpy.ndarray
         Correlation potential of quasiparticles.
     h0_kin_k : numpy.ndarray
         Single-particle dispersion between local clusters. Indexed as 
-        k, orb_i, orb_j.
+        k, orb_i, orb_j
     mu : float, optional
-        Chemical potential.
+        Chemical potential
 
     Return
     ------
-    eigenvalues : numpy.ndarray
-        Indexed as k, band.
-    eigenvectors : numpy.ndarray
-        Indexed as k, each column an eigenvector.
+    h_qp : numpy.ndarray
+        Indexed as k, orb_i, orb_j
     
     Notes
     -----
@@ -280,8 +278,9 @@ def get_h_qp(R : np.ndarray,
         (Lambda - mu*np.eye(Lambda.shape[0]))
     if not np.allclose(h_qp, np.swapaxes(h_qp, 1, 2).conj()):
         warnings.warn("H_qp is not Hermitian !", RuntimeWarning)
-    eig, vec = np.linalg.eigh(h_qp)
-    return (eig, vec)
+    #eig, vec = np.linalg.eigh(h_qp)
+    #return (eig, vec)
+    return h_qp
 
 def get_h0_kin_k_R(R : np.ndarray, 
                    h0_kin_k : np.ndarray, 

src/risb/helpers_triqs.py

@@ -18,7 +18,8 @@
 
 import numpy as np
 from triqs.operators import Operator, c, c_dag
-from risb.helpers import get_h0_loc_matrix
+from triqs.gf import BlockGf, inverse
+from risb.helpers import get_h0_loc_matrix, get_h_qp
 
 def get_C_Op(gf_struct : list[tuple[str,int]], dagger : bool = False) -> dict[list[Operator]]:
     """
@@ -150,4 +151,67 @@ def get_h0_loc(h0_k : dict[np.ndarray],
     h0_loc = Operator()
     for Op in h0_loc_blocks.values():
         h0_loc += Op
-    return h0_loc
+    return h0_loc
+
+def get_gf_struct_from_g(block_gf):
+    gf_struct = []
+    for bl, gf in block_gf:
+        gf_struct.append( [bl, gf.data.shape[-1]] )
+    return gf_struct
+
+def get_sigma_w(gf_struct, mesh, Lambda, R, mu=0, h_loc=None):
+    sigma_w = BlockGf(mesh=mesh, gf_struct=gf_struct)
+    for bl, gf in sigma_w:
+        id = np.eye(*gf.data.shape[-2::]) # Last two indices are the orbital structure
+        Z = R[bl] @ R[bl].conj().T
+        id_Z = (id - np.linalg.inv(Z))
+        id_Z_mu = id_Z * mu
+        hf = np.linalg.inv(R[bl]) @ Lambda[bl] @ np.linalg.inv(R[bl].conj().T)
+        if h_loc is not None:
+            for w in gf.mesh:
+                gf[w] = id_Z * w + hf + id_Z_mu - h_loc[bl]
+        else:
+            for w in gf.mesh:
+                gf[w] = id_Z * w + hf + id_Z_mu
+    return sigma_w
+
+# FIXME have to check h0_kin_k shares the same mesh 
+def get_g_qp_k_w(gf_struct, mesh, h0_kin_k, Lambda, R, mu=0):
+    g_qp_k_w = BlockGf(mesh=mesh, gf_struct=gf_struct)
+    for bl, gf in g_qp_k_w:
+        h_qp = get_h_qp(R=R[bl], Lambda=Lambda[bl], h0_kin_k=h0_kin_k[bl], mu=mu)
+        for k, w in g_qp_k_w.mesh:
+            gf[k,w] = inverse(w - h_qp[k.data_index])
+    return g_qp_k_w
+
+def get_g_k_w(**kwargs):
+    if ('g0_k_w' in kwargs) and ('sigma_w' in kwargs):
+        g0_k_w = kwargs['g0_k_w']
+        sigma_w = kwargs['sigma_w']
+        g_k_w = g0_k_w.copy()
+        g_k_w.zero()
+        for bl, gf in g_k_w:
+            for k, w in gf.mesh:            
+                gf[k,w] = inverse(inverse(g0_k_w[bl][k,w]) - sigma_w[bl][w])
+    elif ('g_qp_k_w' in kwargs) and ('R' in kwargs):
+        g_qp_k_w = kwargs['g_qp_k_w']
+        R = kwargs['R']
+        g_k_w = g_qp_k_w.copy()
+        g_k_w.zero()
+        for bl, gf in g_qp_k_w:
+            g_k_w[bl] = R[bl].conj().T @ gf @ R[bl]
+    else:
+        msg = 'Required kwargs are g0_k_w and sigma_w, or g_qp_k_w and R !'
+        raise ValueError(msg)
+    return g_k_w
+
+def get_g_w_loc(g_k_w):
+    k_mesh = g_k_w.mesh[0]
+    w_mesh = g_k_w.mesh[1]
+    gf_struct = get_gf_struct_from_g(g_k_w)
+    g_w_loc = BlockGf(mesh=w_mesh, gf_struct=gf_struct)
+    for bl, gf in g_w_loc:
+        for k in k_mesh:
+            gf += g_k_w[bl][k,:]
+        gf /= np.prod(k_mesh.dims)
+    return g_w_loc

src/risb/solve_lattice.py

@@ -436,7 +436,8 @@ def one_cycle(self,
         h0_k_R = dict()
         #R_h0_k_R = dict()
         for bl in self.h0_kin_k.keys():
-            self.energies_qp[bl], self.bloch_vector_qp[bl] = helpers.get_h_qp(self.R_full[bl], self.Lambda_full[bl], self.h0_kin_k[bl])
+            h_qp = helpers.get_h_qp(self.R_full[bl], self.Lambda_full[bl], self.h0_kin_k[bl])
+            self.energies_qp[bl], self.bloch_vector_qp[bl] = np.linalg.eigh(h_qp)
             h0_k_R[bl] = helpers.get_h0_kin_k_R(self.R_full[bl], self.h0_kin_k[bl], self.bloch_vector_qp[bl])
             #R_h0_k_R[bl] = helpers.get_R_h0_kin_kR(R_full[bl], self.h0_kin_k[bl], self.bloch_vector_qp[bl])
 

tests/test_helpers.py

@@ -20,7 +20,8 @@ def test_get_h_qp(subtests):
     np.fill_diagonal(R, 1)
     Lambda = np.zeros(shape=(2,2))
     np.fill_diagonal(Lambda, 0.5)
-    eig, vec = helpers.get_h_qp(R, Lambda, h0_k)
+    h_qp = helpers.get_h_qp(R, Lambda, h0_k)
+    eig, vec = np.linalg.eigh(h_qp)
     with subtests.test(msg="eigenvalues"):
         assert eig == approx(eig_expected, abs=abs)
     with subtests.test(msg="eigenvectors"):