13840136129

Committed 13 Mar 2025 04:50PM UTC coverage: 82.713% (-4.1%) from 86.824%

Build # 13840136129

Build Type

Pull #104

github

Committed by

web-flow

Commit Message

Merge 170f405cf into 2a2a1a05f

Pull Request Pull Request #104: Fix TBLattice and DMFT convergence script

Run Details

86 of 106 branches covered (81.13%)

Branch coverage included in aggregate %.

7 of 100 new or added lines in 4 files covered. (7.0%)

1 existing line in 1 file now uncovered.

914 of 1103 relevant lines covered (82.86%)

0.83 hits per line

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

22.03

/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 _check_convergence, _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)
    )

    ust = 0
    jh = 0
    jx = 0
    jp = 0

    # Interaction part
    hamiltonian += h.hubbard_int_orbital_basis[0] * sum(n("up", o) * n("dn", o) for o in orbs)
    hamiltonian += ust * sum(
        int(o1 != o2) * n("up", o1) * n("dn", o2) for o1, o2 in product(orbs, orbs)
    )
    hamiltonian += (ust - jh) * sum(
        int(o1 < o2) * n(s, o1) * n(s, o2) for s, o1, o2 in product(spins, orbs, orbs)
    )
    hamiltonian -= jx * sum(
        int(o1 != o2) * c_dag("up", o1) * c("dn", o1) * c_dag("dn", o2) * c("up", o2)
        for o1, o2 in product(orbs, orbs)
    )
    hamiltonian += jp * sum(
        int(o1 != o2) * c_dag("up", o1) * c_dag("dn", o1) * c("dn", o2) * c("up", o2)
        for o1, o2 in product(orbs, 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,
    )

    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])
        # Check convergence of the Weiss field
        g0 = np.asarray([g0_iw.data, g0_an_iw.data])
        err, converged = _check_convergence(g0, epsilon, n_success, max_iter)

        if converged:
            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✔
17		×
18	from quant_met.mean_field.hamiltonians import BaseHamiltonian	1✔
19	from quant_met.parameters import GenericParameters	1✔
20
21	from .utils import _check_convergence, _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
NEW 83	ust = 0	×
NEW 84	jh = 0	×
NEW 85	jx = 0	×
NEW 86	jp = 0	×
87
88	# Interaction part
NEW 89	hamiltonian += h.hubbard_int_orbital_basis[0] * sum(n("up", o) * n("dn", o) for o in orbs)	×
NEW 90	hamiltonian += ust * sum(	×
91	int(o1 != o2) * n("up", o1) * n("dn", o2) for o1, o2 in product(orbs, orbs)
92	)
NEW 93	hamiltonian += (ust - jh) * sum(	×
94	int(o1 < o2) * n(s, o1) * n(s, o2) for s, o1, o2 in product(spins, orbs, orbs)
95	)
NEW 96	hamiltonian -= jx * sum(	×
97	int(o1 != o2) * c_dag("up", o1) * c("dn", o1) * c_dag("dn", o2) * c("up", o2)
98	for o1, o2 in product(orbs, orbs)
99	)
NEW 100	hamiltonian += jp * sum(	×
101	int(o1 != o2) * c_dag("up", o1) * c_dag("dn", o1) * c("dn", o2) * c("up", o2)
102	for o1, o2 in product(orbs, orbs)
103	)
104
105	# Matrix dimensions of eps and V: 3 orbitals x 2 bath states
106	eps = np.array([[-1.0, -0.5, 0.5, 1.0] for _ in range(tbl.n_orbitals)])	×
107	v = 0.5 * np.ones((tbl.n_orbitals, n_bath))	×
108	d = -0.2 * np.eye(tbl.n_orbitals * n_bath)	×
109
110	# Bath
111	hamiltonian += sum(	×
112	eps[o, nu] * c_dag("B_" + s, o * n_bath + nu) * c("B_" + s, o * n_bath + nu)
113	for s, o, nu in product(spins, orbs, range(n_bath))
114	)
115
116	hamiltonian += sum(	×
117	v[o, nu]
118	* (c_dag(s, o) * c("B_" + s, o * n_bath + nu) + c_dag("B_" + s, o * n_bath + nu) * c(s, o))
119	for s, o, nu in product(spins, orbs, range(n_bath))
120	)
121
122	# Anomalous bath
123	hamiltonian += sum(	×
124	d[o, q] * (c("B_up", o) * c("B_dn", q)) + dagger(d[o, q] * (c("B_up", o) * c("B_dn", q)))
125	for o, q in product(range(tbl.n_orbitals * n_bath), range(tbl.n_orbitals * n_bath))
126	)
127
128	# Create solver object
129	fit_params = BathFittingParams(method="minimize", grad="numeric")	×
130	solver = EDIpackSolver(	×
131	hamiltonian,
132	fops_imp_up,
133	fops_imp_dn,
134	fops_bath_up,
135	fops_bath_dn,
136	lanc_dim_threshold=1024,
137	verbose=1,
138	bath_fitting_params=fit_params,
139	)
140
141	for iloop in range(max_iter):	×
142	print(f"\nLoop {iloop + 1} of {max_iter}")	×
143
144	# Solve the effective impurity problem
145	solver.solve(	×
146	beta=h.beta,
147	n_iw=n_iw,
148	energy_window=energy_window,
149	n_w=n_w,
150	broadening=broadening,
151	)
152
153	# Normal and anomalous components of computed self-energy
154	s_iw = solver.Sigma_iw["up"]	×
155	s_an_iw = solver.Sigma_an_iw["up_dn"]	×
156
157	# Compute local Green's function
158	g_iw, g_an_iw = get_gloc(s_iw, s_an_iw, h0_nambu_k, xmu, broadening, kmesh)	×
159	# Compute Weiss field
160	g0_iw, g0_an_iw = _dmft_weiss_field(g_iw, g_an_iw, s_iw, s_an_iw)	×
161
162	# Bath fitting and mixing
163	g0_iw_full = BlockGf(name_list=spins, block_list=[g0_iw, g0_iw])	×
164	g0_an_iw_full = BlockGf(name_list=["up_dn"], block_list=[g0_an_iw])	×
165
166	bath_new = solver.chi2_fit_bath(g0_iw_full, g0_an_iw_full)[0]	×
167	solver.bath = w_mixing * bath_new + (1 - w_mixing) * solver.bath	×
168
169	# Check convergence of the Weiss field
170	g0 = np.asarray([g0_iw.data, g0_an_iw.data])	×
171	# Check convergence of the Weiss field
NEW 172	g0 = np.asarray([g0_iw.data, g0_an_iw.data])	×
NEW 173	err, converged = _check_convergence(g0, epsilon, n_success, max_iter)	×
174
NEW 175	if converged:	×
176	break	×
177
178	return solver	×

Ruberhauptmann / quant-met / 13840136129

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