From f32a4656b6f076d244867b1abead21f9188e5958 Mon Sep 17 00:00:00 2001 From: "H. L. Nourse" <32994022+thenoursehorse@users.noreply.github.com> Date: Tue, 14 May 2024 10:25:26 +0900 Subject: [PATCH] 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 | 4 +- .github/dependabot.yml | 11 +++++ .github/workflows/ci.yml | 3 +- .gitignore | 4 +- docs/tutorials/self-consistent.md | 9 ++-- examples/hubbard.py | 35 +++++++++++++ pyproject.toml | 3 +- src/risb/helpers.py | 19 ++++--- src/risb/helpers_triqs.py | 68 +++++++++++++++++++++++++- src/risb/solve_lattice.py | 3 +- tests/test_helpers.py | 3 +- 11 files changed, 138 insertions(+), 24 deletions(-) create mode 100644 .github/dependabot.yml diff --git a/.github/actions/setup-triqs/action.yml b/.github/actions/setup-triqs/action.yml index 1d6e236..93f7856 100644 --- a/.github/actions/setup-triqs/action.yml +++ b/.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 }} diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..4556ff4 --- /dev/null +++ b/.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: + - "*" \ No newline at end of file diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9ea48cb..c9924a6 100644 --- a/.github/workflows/ci.yml +++ b/.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"} diff --git a/.gitignore b/.gitignore index 8700ad5..ad88fd9 100644 --- a/.gitignore +++ b/.gitignore @@ -16,4 +16,6 @@ _build .github.auth # Visual Studio Code -.vscode \ No newline at end of file +.vscode + +src/risb/version.py \ No newline at end of file diff --git a/docs/tutorials/self-consistent.md b/docs/tutorials/self-consistent.md index ad69898..84d7edb 100644 --- a/docs/tutorials/self-consistent.md +++ b/docs/tutorials/self-consistent.md @@ -230,9 +230,9 @@ 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 @@ -240,7 +240,8 @@ 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 diff --git a/examples/hubbard.py b/examples/hubbard.py index 3bff35a..10b2e94 100644 --- a/examples/hubbard.py +++ b/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') diff --git a/pyproject.toml b/pyproject.toml index 3a422e8..cfc7534 100644 --- a/pyproject.toml +++ b/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 = [ diff --git a/src/risb/helpers.py b/src/risb/helpers.py index d535880..20836b9 100644 --- a/src/risb/helpers.py +++ b/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, diff --git a/src/risb/helpers_triqs.py b/src/risb/helpers_triqs.py index 812ca19..72d5752 100644 --- a/src/risb/helpers_triqs.py +++ b/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 \ No newline at end of file + 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 \ No newline at end of file diff --git a/src/risb/solve_lattice.py b/src/risb/solve_lattice.py index 76fb0f2..1600eca 100644 --- a/src/risb/solve_lattice.py +++ b/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]) diff --git a/tests/test_helpers.py b/tests/test_helpers.py index 35191f0..8ba8a42 100644 --- a/tests/test_helpers.py +++ b/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"):