13764497520

Committed 10 Mar 2025 12:21PM UTC coverage: 86.835% (-0.09%) from 86.922%

Build # 13764497520

Build Type

push

github

Committed by

web-flow

Commit Message

Fix mistake in chemical potential (#100)

Run Details

87 of 100 branches covered (87.0%)

Branch coverage included in aggregate %.

3 of 4 new or added lines in 3 files covered. (75.0%)

1 existing line in 1 file now uncovered.

909 of 1047 relevant lines covered (86.82%)

0.87 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

16.46

/src/quant_met/dmft/dmft_loop.py

# SPDX-FileCopyrightText: 2024 Tjark Sievers
#
# SPDX-License-Identifier: MIT

"""Functions to run self-consistent calculation for the order parameter."""

import logging
from itertools import product

import numpy as np
import numpy.typing as npt
from edipack2triqs.fit import BathFittingParams
from edipack2triqs.solver import EDIpackSolver
from triqs.gf import BlockGf, Gf, MeshBrZone
from triqs.lattice.tight_binding import TBLattice
from triqs.operators import c, c_dag, dagger, n

from quant_met.mean_field.hamiltonians import BaseHamiltonian
from quant_met.parameters import GenericParameters

from .utils import _dmft_weiss_field, get_gloc

logger = logging.getLogger(__name__)


def dmft_loop(
    tbl: TBLattice,
    h: BaseHamiltonian[GenericParameters],
    h0_nambu_k: Gf,
    n_bath: float,
    n_iw: int,
    broadening: float,
    n_w: int,
    w_mixing: float,
    n_success: int,
    xmu: npt.NDArray[np.float64],
    kmesh: MeshBrZone,
    epsilon: float,
    max_iter: int,
) -> EDIpackSolver:
    """DMFT loop.

    Parameters
    ----------
    tbl
    h
    h0_nambu_k
    n_bath
    n_iw
    broadening
    n_w
    w_mixing
    n_success
    xmu
    kmesh
    epsilon
    max_iter

    Returns
    -------
    EDIpackSolver

    """
    energy_window = (-2.0 * h.hopping_gr, 2.0 * h.hopping_gr)

    spins = ("up", "dn")
    orbs = range(tbl.n_orbitals)

    # Fundamental sets for impurity degrees of freedom
    fops_imp_up = [("up", o) for o in orbs]
    fops_imp_dn = [("dn", o) for o in orbs]

    # Fundamental sets for bath degrees of freedom
    fops_bath_up = [("B_up", i) for i in range(tbl.n_orbitals * n_bath)]
    fops_bath_dn = [("B_dn", i) for i in range(tbl.n_orbitals * n_bath)]

    # Non-interacting part of the impurity Hamiltonian
    h_loc = -xmu * np.eye(tbl.n_orbitals)
    hamiltonian = sum(
        h_loc[o1, o2] * c_dag(spin, o1) * c(spin, o2) for spin, o1, o2 in product(spins, orbs, orbs)
    )

    # Interaction part
    hamiltonian += -h.hubbard_int_orbital_basis[0] * sum(n("up", o) * n("dn", o) for o in orbs)

    # Matrix dimensions of eps and V: 3 orbitals x 2 bath states
    eps = np.array([[-1.0, -0.5, 0.5, 1.0] for _ in range(tbl.n_orbitals)])
    v = 0.5 * np.ones((tbl.n_orbitals, n_bath))
    d = -0.2 * np.eye(tbl.n_orbitals * n_bath)

    # Bath
    hamiltonian += sum(
        eps[o, nu] * c_dag("B_" + s, o * n_bath + nu) * c("B_" + s, o * n_bath + nu)
        for s, o, nu in product(spins, orbs, range(n_bath))
    )

    hamiltonian += sum(
        v[o, nu]
        * (c_dag(s, o) * c("B_" + s, o * n_bath + nu) + c_dag("B_" + s, o * n_bath + nu) * c(s, o))
        for s, o, nu in product(spins, orbs, range(n_bath))
    )

    # Anomalous bath
    hamiltonian += sum(
        d[o, q] * (c("B_up", o) * c("B_dn", q)) + dagger(d[o, q] * (c("B_up", o) * c("B_dn", q)))
        for o, q in product(range(tbl.n_orbitals * n_bath), range(tbl.n_orbitals * n_bath))
    )

    # Create solver object
    fit_params = BathFittingParams(method="minimize", grad="numeric")
    solver = EDIpackSolver(
        hamiltonian,
        fops_imp_up,
        fops_imp_dn,
        fops_bath_up,
        fops_bath_dn,
        lanc_dim_threshold=1024,
        verbose=1,
        bath_fitting_params=fit_params,
    )

    gooditer = 0
    g0_prev = np.zeros((2, 2 * n_iw, tbl.n_orbitals, tbl.n_orbitals), dtype=complex)
    for iloop in range(max_iter):
        print(f"\nLoop {iloop + 1} of {max_iter}")

        # Solve the effective impurity problem
        solver.solve(
            beta=h.beta,
            n_iw=n_iw,
            energy_window=energy_window,
            n_w=n_w,
            broadening=broadening,
        )

        # Normal and anomalous components of computed self-energy
        s_iw = solver.Sigma_iw["up"]
        s_an_iw = solver.Sigma_an_iw["up_dn"]

        # Compute local Green's function
        g_iw, g_an_iw = get_gloc(s_iw, s_an_iw, h0_nambu_k, xmu, broadening, kmesh)
        # Compute Weiss field
        g0_iw, g0_an_iw = _dmft_weiss_field(g_iw, g_an_iw, s_iw, s_an_iw)

        # Bath fitting and mixing
        g0_iw_full = BlockGf(name_list=spins, block_list=[g0_iw, g0_iw])
        g0_an_iw_full = BlockGf(name_list=["up_dn"], block_list=[g0_an_iw])

        bath_new = solver.chi2_fit_bath(g0_iw_full, g0_an_iw_full)[0]
        solver.bath = w_mixing * bath_new + (1 - w_mixing) * solver.bath

        # Check convergence of the Weiss field
        g0 = np.asarray([g0_iw.data, g0_an_iw.data])
        errvec = np.real(np.sum(abs(g0 - g0_prev), axis=1) / np.sum(abs(g0), axis=1))
        # First iteration
        if iloop == 0:
            errvec = np.ones_like(errvec)
        errmin, err, errmax = np.min(errvec), np.average(errvec), np.max(errvec)

        g0_prev = np.copy(g0)

        if err < epsilon:
            gooditer += 1  # Increase good iterations count
        else:
            gooditer = 0  # Reset good iterations count

        conv_bool = ((err < epsilon) and (gooditer > n_success) and (iloop < max_iter)) or (
            iloop >= max_iter
        )

        # Print convergence message
        if iloop < max_iter:
            if errvec.size > 1:
                print(f"max error={errmax:.6e}")
            print("    " * (errvec.size > 1) + f"error={err:.6e}")
            if errvec.size > 1:
                print(f"min error={errmin:.6e}")
        else:
            if errvec.size > 1:
                print(f"max error={errmax:.6e}")
            print("    " * (errvec.size > 1) + f"error={err:.6e}")
            if errvec.size > 1:
                print(f"min error={errmin:.6e}")
            print(f"Not converged after {max_iter} iterations.")

        if conv_bool:
            break

    return solver

1	# SPDX-FileCopyrightText: 2024 Tjark Sievers
2	#
3	# SPDX-License-Identifier: MIT
4
5	"""Functions to run self-consistent calculation for the order parameter."""
6
7	import logging	1✔
8	from itertools import product	1✔
9
10	import numpy as np	1✔
11	import numpy.typing as npt	1✔
12	from edipack2triqs.fit import BathFittingParams	1✔
13	from edipack2triqs.solver import EDIpackSolver	1✔
14	from triqs.gf import BlockGf, Gf, MeshBrZone	1✔
15	from triqs.lattice.tight_binding import TBLattice
16	from triqs.operators import c, c_dag, dagger, n	1✔
UNCOV 17		×
18	from quant_met.mean_field.hamiltonians import BaseHamiltonian	1✔
19	from quant_met.parameters import GenericParameters	1✔
20
21	from .utils import _dmft_weiss_field, get_gloc	1✔
22
23	logger = logging.getLogger(__name__)	1✔
24
25
26	def dmft_loop(	1✔
27	tbl: TBLattice,
28	h: BaseHamiltonian[GenericParameters],
29	h0_nambu_k: Gf,
30	n_bath: float,
31	n_iw: int,
32	broadening: float,
33	n_w: int,
34	w_mixing: float,
35	n_success: int,
36	xmu: npt.NDArray[np.float64],
37	kmesh: MeshBrZone,
38	epsilon: float,
39	max_iter: int,
40	) -> EDIpackSolver:
41	"""DMFT loop.
42
43	Parameters
44	----------
45	tbl
46	h
47	h0_nambu_k
48	n_bath
49	n_iw
50	broadening
51	n_w
52	w_mixing
53	n_success
54	xmu
55	kmesh
56	epsilon
57	max_iter
58
59	Returns
60	-------
61	EDIpackSolver
62
63	"""
64	energy_window = (-2.0 * h.hopping_gr, 2.0 * h.hopping_gr)	×
65
66	spins = ("up", "dn")	×
67	orbs = range(tbl.n_orbitals)	×
68
69	# Fundamental sets for impurity degrees of freedom
70	fops_imp_up = [("up", o) for o in orbs]	×
71	fops_imp_dn = [("dn", o) for o in orbs]	×
72
73	# Fundamental sets for bath degrees of freedom
74	fops_bath_up = [("B_up", i) for i in range(tbl.n_orbitals * n_bath)]	×
75	fops_bath_dn = [("B_dn", i) for i in range(tbl.n_orbitals * n_bath)]	×
76
77	# Non-interacting part of the impurity Hamiltonian
78	h_loc = -xmu * np.eye(tbl.n_orbitals)	×
79	hamiltonian = sum(	×
80	h_loc[o1, o2] * c_dag(spin, o1) * c(spin, o2) for spin, o1, o2 in product(spins, orbs, orbs)
81	)
82
83	# Interaction part
84	hamiltonian += -h.hubbard_int_orbital_basis[0] * sum(n("up", o) * n("dn", o) for o in orbs)	×
85
86	# Matrix dimensions of eps and V: 3 orbitals x 2 bath states
87	eps = np.array([[-1.0, -0.5, 0.5, 1.0] for _ in range(tbl.n_orbitals)])	×
88	v = 0.5 * np.ones((tbl.n_orbitals, n_bath))	×
89	d = -0.2 * np.eye(tbl.n_orbitals * n_bath)	×
90
91	# Bath
92	hamiltonian += sum(	×
93	eps[o, nu] * c_dag("B_" + s, o * n_bath + nu) * c("B_" + s, o * n_bath + nu)
94	for s, o, nu in product(spins, orbs, range(n_bath))
95	)
96
97	hamiltonian += sum(	×
98	v[o, nu]
99	* (c_dag(s, o) * c("B_" + s, o * n_bath + nu) + c_dag("B_" + s, o * n_bath + nu) * c(s, o))
100	for s, o, nu in product(spins, orbs, range(n_bath))
101	)
102
103	# Anomalous bath
104	hamiltonian += sum(	×
105	d[o, q] * (c("B_up", o) * c("B_dn", q)) + dagger(d[o, q] * (c("B_up", o) * c("B_dn", q)))
106	for o, q in product(range(tbl.n_orbitals * n_bath), range(tbl.n_orbitals * n_bath))
107	)
108
109	# Create solver object
110	fit_params = BathFittingParams(method="minimize", grad="numeric")	×
111	solver = EDIpackSolver(	×
112	hamiltonian,
113	fops_imp_up,
114	fops_imp_dn,
115	fops_bath_up,
116	fops_bath_dn,
117	lanc_dim_threshold=1024,
118	verbose=1,
119	bath_fitting_params=fit_params,
120	)
121
122	gooditer = 0	×
123	g0_prev = np.zeros((2, 2 * n_iw, tbl.n_orbitals, tbl.n_orbitals), dtype=complex)	×
124	for iloop in range(max_iter):	×
125	print(f"\nLoop {iloop + 1} of {max_iter}")	×
126
127	# Solve the effective impurity problem
128	solver.solve(	×
129	beta=h.beta,
130	n_iw=n_iw,
131	energy_window=energy_window,
132	n_w=n_w,
133	broadening=broadening,
134	)
135
136	# Normal and anomalous components of computed self-energy
137	s_iw = solver.Sigma_iw["up"]	×
138	s_an_iw = solver.Sigma_an_iw["up_dn"]	×
139
140	# Compute local Green's function
141	g_iw, g_an_iw = get_gloc(s_iw, s_an_iw, h0_nambu_k, xmu, broadening, kmesh)	×
142	# Compute Weiss field
143	g0_iw, g0_an_iw = _dmft_weiss_field(g_iw, g_an_iw, s_iw, s_an_iw)	×
144
145	# Bath fitting and mixing
146	g0_iw_full = BlockGf(name_list=spins, block_list=[g0_iw, g0_iw])	×
147	g0_an_iw_full = BlockGf(name_list=["up_dn"], block_list=[g0_an_iw])	×
148
149	bath_new = solver.chi2_fit_bath(g0_iw_full, g0_an_iw_full)[0]	×
150	solver.bath = w_mixing * bath_new + (1 - w_mixing) * solver.bath	×
151
152	# Check convergence of the Weiss field
153	g0 = np.asarray([g0_iw.data, g0_an_iw.data])	×
154	errvec = np.real(np.sum(abs(g0 - g0_prev), axis=1) / np.sum(abs(g0), axis=1))	×
155	# First iteration
156	if iloop == 0:	×
157	errvec = np.ones_like(errvec)	×
158	errmin, err, errmax = np.min(errvec), np.average(errvec), np.max(errvec)	×
159
160	g0_prev = np.copy(g0)	×
161
162	if err < epsilon:	×
163	gooditer += 1 # Increase good iterations count	×
164	else:
165	gooditer = 0 # Reset good iterations count	×
166
167	conv_bool = ((err < epsilon) and (gooditer > n_success) and (iloop < max_iter)) or (	×
168	iloop >= max_iter
169	)
170
171	# Print convergence message
172	if iloop < max_iter:	×
173	if errvec.size > 1:	×
174	print(f"max error={errmax:.6e}")	×
175	print(" " * (errvec.size > 1) + f"error={err:.6e}")	×
176	if errvec.size > 1:	×
177	print(f"min error={errmin:.6e}")	×
178	else:
179	if errvec.size > 1:	×
180	print(f"max error={errmax:.6e}")	×
181	print(" " * (errvec.size > 1) + f"error={err:.6e}")	×
182	if errvec.size > 1:	×
183	print(f"min error={errmin:.6e}")	×
184	print(f"Not converged after {max_iter} iterations.")	×
185
186	if conv_bool:	×
187	break	×
188
189	return solver	×

Ruberhauptmann / quant-met / 13764497520

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous