Usage examples

3 orbital impurity with Hubbard-Kanamori interaction

  1from itertools import product
  2import numpy as np
  3from mpi4py import MPI
  4
  5# TRIQS many-body operator objects
  6from triqs.operators import c, c_dag, n
  7
  8# Compatibility layer between EDIpack and TRIQS
  9from edipack2triqs.solver import EDIpackSolver
 10
 11# Fundamental sets for impurity degrees of freedom
 12fops_imp_up = [('up', 0), ('up', 1), ('up', 2)]
 13fops_imp_dn = [('dn', 0), ('dn', 1), ('dn', 2)]
 14
 15# Fundamental sets for bath degrees of freedom
 16fops_bath_up = [('B_up', i) for i in range(3 * 2)]
 17fops_bath_dn = [('B_dn', i) for i in range(3 * 2)]
 18
 19# Define the Hamiltonian
 20orbs = range(3)
 21
 22## Non-interacting part of the impurity Hamiltonian
 23h_loc = np.diag([-0.7, -0.5, -0.7])
 24H = sum(h_loc[o1, o2] * c_dag(spin, o1) * c(spin, o2)
 25        for spin, o1, o2 in product(('up', 'dn'), orbs, orbs))
 26
 27## Interaction part
 28U = 3.0     # Local intra-orbital interactions U
 29Ust = 1.2   # Local inter-orbital interaction U'
 30Jh = 0.2    # Hund's coupling
 31Jx = 0.15   # Spin-exchange coupling constant
 32Jp = 0.1    # Pair-hopping coupling constant
 33
 34H += U * sum(n('up', o) * n('dn', o) for o in orbs)
 35H += Ust * sum(int(o1 != o2) * n('up', o1) * n('dn', o2)
 36               for o1, o2 in product(orbs, orbs))
 37H += (Ust - Jh) * sum(int(o1 < o2) * n(s, o1) * n(s, o2)
 38                      for s, o1, o2 in product(('up', 'dn'), orbs, orbs))
 39H -= Jx * sum(int(o1 != o2)
 40              * c_dag('up', o1) * c('dn', o1) * c_dag('dn', o2) * c('up', o2)
 41              for o1, o2 in product(orbs, orbs))
 42H += Jp * sum(int(o1 != o2)
 43              * c_dag('up', o1) * c_dag('dn', o1) * c('dn', o2) * c('up', o2)
 44              for o1, o2 in product(orbs, orbs))
 45
 46## Bath part
 47
 48# Matrix dimensions of eps and V: 3 orbitals x 2 bath states
 49eps = np.array([[-0.1, 0.1],
 50                [-0.2, 0.2],
 51                [-0.3, 0.3]])
 52V = np.array([[0.1, 0.2],
 53              [0.1, 0.2],
 54              [0.1, 0.2]])
 55
 56H += sum(eps[o, nu] * c_dag("B_" + s, nu * 3 + o) * c("B_" + s, nu * 3 + o)
 57         for s, o, nu in product(('up', 'dn'), orbs, range(2)))
 58H += sum(V[o, nu] * (c_dag(s, o) * c("B_" + s, nu * 3 + o)
 59                     + c_dag("B_" + s, nu * 3 + o) * c(s, o))
 60         for s, o, nu in product(('up', 'dn'), orbs, range(2)))
 61
 62# Create a solver object that wraps EDIpack functionality
 63# See help(EDIpackSolver) for a complete list of parameters
 64solver = EDIpackSolver(H,
 65                       fops_imp_up, fops_imp_dn, fops_bath_up, fops_bath_dn,
 66                       lanc_dim_threshold=16)
 67
 68# Solve the impurity model
 69beta = 100.0  # Inverse temperature
 70n_iw = 1024   # Number of Matsubara frequencies for impurity GF calculations
 71energy_window = (-2.0, 2.0)  # Energy window for real-frequency GF calculations
 72n_w = 4000    # Number of real-frequency points for impurity GF calculations
 73broadening = 0.005  # Broadening on the real axis for impurity GF calculations
 74
 75solver.solve(beta=beta,
 76             n_iw=n_iw,
 77             energy_window=energy_window,
 78             n_w=n_w,
 79             broadening=broadening)
 80
 81# On master node, output results of calculation
 82if MPI.COMM_WORLD.Get_rank() == 0:
 83    print("Using SciFortran", solver.scifor_version)
 84    print("Using EDIpack", solver.edipack_version)
 85
 86    print("Densities (per orbital):", solver.densities)
 87    print("Double occupancy (per orbital):", solver.double_occ)
 88    print("Magnetization M_x (per orbital):", solver.magnetization[:, 0])
 89    print("Magnetization M_y (per orbital):", solver.magnetization[:, 1])
 90    print("Magnetization M_z (per orbital):", solver.magnetization[:, 2])
 91
 92    from triqs.plot.mpl_interface import plt, oplot
 93
 94    # Extract the Matsubara GF from the solver and plot it
 95    g_iw = solver.g_iw
 96
 97    plt.figure()
 98    for spin in ('up', 'dn'):
 99        for orb in range(3):
100            label = r"g_{%s,%i%i}(i\omega_n)$" % (spin, orb, orb)
101            oplot(g_iw[spin][orb, orb].real, label=r"$\Re " + label)
102            oplot(g_iw[spin][orb, orb].imag, label=r"$\Im " + label)
103    plt.legend(loc='lower center', ncols=3)
104    plt.savefig("g_iw.pdf")
105
106    # Extract the Matsubara self-energy from the solver and plot it
107    Sigma_iw = solver.Sigma_iw
108
109    plt.figure()
110    for spin in ('up', 'dn'):
111        for orb in range(3):
112            label = r"\Sigma_{%s,%i%i}(i\omega_n)$" % (spin, orb, orb)
113            oplot(Sigma_iw[spin][orb, orb].real, label=r"$\Re " + label)
114            oplot(Sigma_iw[spin][orb, orb].imag, label=r"$\Im " + label)
115    plt.legend(loc='lower center', ncols=3)
116    plt.savefig("Sigma_iw.pdf")
117
118    # Extract the real-frequency GF from the solver and plot it
119    g_w = solver.g_w
120
121    plt.figure()
122    for spin in ('up', 'dn'):
123        for orb in range(3):
124            label = r"g_{%s,%i%i}(\omega)$" % (spin, orb, orb)
125            oplot(g_w[spin][orb, orb].real, label=r"$\Re " + label)
126            oplot(g_w[spin][orb, orb].imag, label=r"$\Im " + label)
127    plt.legend(loc='lower center', ncols=3)
128    plt.savefig("g_w.pdf")
129
130    # Extract the real-frequency self-energy from the solver and plot it
131    Sigma_w = solver.Sigma_w
132
133    plt.figure()
134    for spin in ('up', 'dn'):
135        for orb in range(3):
136            label = r"\Sigma_{%s,%i%i}(\omega)$" % (spin, orb, orb)
137            oplot(Sigma_w[spin][orb, orb].real, label=r"$\Re " + label)
138            oplot(Sigma_w[spin][orb, orb].imag, label=r"$\Im " + label)
139    plt.legend(loc='lower center', ncols=3)
140    plt.savefig("Sigma_w.pdf")
141
142# It is possible to change some parameters of the Hamiltonian and to diagonalize
143# it again.
144
145# Change a matrix element of h_loc
146# Since the non-interacting part of the original Hamiltonian is spin-symmetric,
147# solver.hloc() returns an array of shape (1, 1, 3, 3). In presence of a spin
148# energy splitting the shape would be (2, 2, 3, 3).
149solver.hloc[0, 0, 1, 1] = -0.6  # spin1 = spin2 = up and down, orb1 = orb2 = 1
150
151# Change elements of the interaction matrix corresponding to Jp
152# The order of indices of solver.U is
153# orb1, spin1, orb2, spin2, orb3, spin3, orb4, spin4.
154Jp = 0.15
155for s, o1, o2 in product(range(2), orbs, orbs):
156    if o1 == o2:
157        continue
158    solver.U[o1, s, o1, 1 - s, o2, s, o2, 1 - s] = 0.5 * Jp
159    solver.U[o1, s, o1, 1 - s, o2, 1 - s, o2, s] = -0.5 * Jp
160
161# Update some bath parameters
162solver.bath.eps[0, 2, 1] = -0.4  # spin = up and down, orb1 = 2, nu = 1
163solver.bath.V[0, 0, 1] = 0.3     # spin = up and down, orb1 = 0, nu = 1
164
165solver.solve(beta=beta)

DMFT study of superconductivity in the attractive Hubbard model

  1from itertools import product
  2import numpy as np
  3
  4# TRIQS modules
  5from triqs.gf import (Gf, BlockGf, MeshImFreq, MeshProduct,
  6                      iOmega_n, conjugate, inverse)
  7from triqs.gf.tools import dyson
  8from triqs.operators import c, c_dag, n, dagger
  9from triqs.lattice.tight_binding import TBLattice
 10from h5 import HDFArchive
 11
 12# edipack2triqs modules
 13from edipack2triqs.solver import EDIpackSolver
 14from edipack2triqs.fit import BathFittingParams
 15
 16
 17# Parameters
 18Nspin = 1
 19Nloop = 100
 20Nsuccess = 1
 21threshold = 1e-6
 22wmixing = 0.5
 23Norb = 1
 24Nbath = 4
 25n_k = 20
 26U = -3.0      # Local intra-orbital interactions U
 27Ust = 1.2     # Local inter-orbital interaction U'
 28Jh = 0.2      # Hund's coupling
 29Jx = 0.15     # Spin-exchange coupling constant
 30Jp = 0.1      # Pair-hopping coupling constant
 31beta = 100.0  # Inverse temperature
 32n_iw = 1024   # Number of Matsubara frequencies for impurity GF calculations
 33t = 0.5
 34# Energy window for real-frequency GF calculations
 35energy_window = (-2.0 * t, 2.0 * t)
 36n_w = 4000    # Number of real-frequency points for impurity GF calculations
 37broadening = 0.005  # Broadening on the real axis for impurity GF calculations
 38
 39spins = ('up', 'dn')
 40orbs = range(Norb)
 41
 42# Fundamental sets for impurity degrees of freedom
 43fops_imp_up = [('up', o) for o in orbs]
 44fops_imp_dn = [('dn', o) for o in orbs]
 45
 46# Fundamental sets for bath degrees of freedom
 47fops_bath_up = [('B_up', i) for i in range(Norb * Nbath)]
 48fops_bath_dn = [('B_dn', i) for i in range(Norb * Nbath)]
 49
 50# Non-interacting part of the impurity Hamiltonian
 51xmu = U / 2 + (Norb - 1) * Ust / 2 + (Norb - 1) * (Ust - Jh) / 2
 52
 53h_loc = -xmu * np.eye(Norb)
 54H = sum(h_loc[o1, o2] * c_dag(spin, o1) * c(spin, o2)
 55        for spin, o1, o2 in product(spins, orbs, orbs))
 56
 57# Interaction part
 58H += U * sum(n('up', o) * n('dn', o) for o in orbs)
 59H += Ust * sum(int(o1 != o2) * n('up', o1) * n('dn', o2)
 60               for o1, o2 in product(orbs, orbs))
 61H += (Ust - Jh) * sum(int(o1 < o2) * n(s, o1) * n(s, o2)
 62                      for s, o1, o2 in product(spins, orbs, orbs))
 63H -= Jx * sum(int(o1 != o2)
 64              * c_dag('up', o1) * c('dn', o1) * c_dag('dn', o2) * c('up', o2)
 65              for o1, o2 in product(orbs, orbs))
 66H += Jp * sum(int(o1 != o2)
 67              * c_dag('up', o1) * c_dag('dn', o1) * c('dn', o2) * c('up', o2)
 68              for o1, o2 in product(orbs, orbs))
 69
 70# Lattice Hamiltonian
 71units = [(1, 0, 0), (0, 1, 0)]
 72orbital_positions = [(0.0, 0.0, 0.0) for i in range(Norb)]
 73hoppings = {(1, 0): t * np.eye(Norb),
 74            (-1, 0): t * np.eye(Norb),
 75            (0, 1): t * np.eye(Norb),
 76            (0, -1): t * np.eye(Norb)}
 77
 78TBL = TBLattice(units=units,
 79                hoppings=hoppings,
 80                orbital_positions=orbital_positions)
 81kmesh = TBL.get_kmesh(n_k=(n_k, n_k, 1))
 82enk = TBL.fourier(kmesh)
 83
 84nambu_shape = (2 * Norb, 2 * Norb)
 85h0_nambu_k = Gf(mesh=kmesh, target_shape=nambu_shape)
 86for k in kmesh:
 87    h0_nambu_k[k][:Norb, :Norb] = enk(k)
 88    h0_nambu_k[k][Norb:, Norb:] = -enk(-k)
 89
 90# Matrix dimensions of eps and V: 3 orbitals x 2 bath states
 91eps = np.array([[-1.0, -0.5, 0.5, 1.0] for i in range(Norb)])
 92V = 0.5 * np.ones((Norb, Nbath))
 93D = -0.2 * np.eye(Norb * Nbath)
 94
 95# Bath
 96H += sum(eps[o, nu]
 97         * c_dag("B_" + s, o * Nbath + nu) * c("B_" + s, o * Nbath + nu)
 98         for s, o, nu in product(spins, orbs, range(Nbath)))
 99
100H += sum(V[o, nu] * (c_dag(s, o) * c("B_" + s, o * Nbath + nu)
101                     + c_dag("B_" + s, o * Nbath + nu) * c(s, o))
102         for s, o, nu in product(spins, orbs, range(Nbath)))
103
104# Anomalous bath
105H += sum(D[o, q] * (c('B_up', o) * c('B_dn', q))
106         + dagger(D[o, q] * (c('B_up', o) * c('B_dn', q)))
107         for o, q in product(range(Norb * Nbath), range(Norb * Nbath)))
108
109# Create solver object
110fit_params = BathFittingParams(method="minimize", grad="numeric")
111solver = EDIpackSolver(H,
112                       fops_imp_up, fops_imp_dn, fops_bath_up, fops_bath_dn,
113                       lanc_dim_threshold=1024,
114                       verbose=1,
115                       bath_fitting_params=fit_params)
116
117
118# Compute local GF from bare lattice Hamiltonian and self-energy
119def get_gloc(s, s_an):
120    z = Gf(mesh=s.mesh, target_shape=nambu_shape)
121    if isinstance(s.mesh, MeshImFreq):
122        z[:Norb, :Norb] << iOmega_n + xmu - s
123        z[:Norb, Norb:] << -s_an
124        z[Norb:, :Norb] << -s_an
125        z[Norb:, Norb:] << iOmega_n - xmu + conjugate(s)
126    else:
127        z[:Norb, Norb:] << -s_an
128        z[Norb:, :Norb] << -s_an
129        for w in z.mesh:
130            z[w][:Norb, :Norb] = \
131                (w + 1j * broadening + xmu) * np.eye(Norb) - s[w]
132            z[w][Norb:, Norb:] = \
133                (w + 1j * broadening - xmu) * np.eye(Norb) + conjugate(s(-w))
134
135    g_k = Gf(mesh=MeshProduct(kmesh, z.mesh), target_shape=nambu_shape)
136    for k in kmesh:
137        g_k[k, :] << inverse(z - h0_nambu_k[k])
138
139    g_loc_nambu = sum(g_k[k, :] for k in kmesh) / len(kmesh)
140
141    g_loc = s.copy()
142    g_loc_an = s_an.copy()
143    g_loc[:] = g_loc_nambu[:Norb, :Norb]
144    g_loc_an[:] = g_loc_nambu[:Norb, Norb:]
145    return g_loc, g_loc_an
146
147
148# Compute Weiss field from local GF and self-energy
149def dmft_weiss_field(g_iw, g_an_iw, s_iw, s_an_iw):
150    g_nambu_iw = Gf(mesh=g_iw.mesh, target_shape=nambu_shape)
151    s_nambu_iw = Gf(mesh=s_iw.mesh, target_shape=nambu_shape)
152
153    g_nambu_iw[:Norb, :Norb] = g_iw
154    g_nambu_iw[:Norb, Norb:] = g_an_iw
155    g_nambu_iw[Norb:, :Norb] = g_an_iw
156    g_nambu_iw[Norb:, Norb:] = -conjugate(g_iw)
157
158    s_nambu_iw[:Norb, :Norb] = s_iw
159    s_nambu_iw[:Norb, Norb:] = s_an_iw
160    s_nambu_iw[Norb:, :Norb] = s_an_iw
161    s_nambu_iw[Norb:, Norb:] = -conjugate(s_iw)
162
163    g0_nambu_iw = dyson(G_iw=g_nambu_iw, Sigma_iw=s_nambu_iw)
164
165    g0_iw = g_iw.copy()
166    g0_an_iw = g_an_iw.copy()
167    g0_iw[:] = g0_nambu_iw[:Norb, :Norb]
168    g0_an_iw[:] = g0_nambu_iw[:Norb, Norb:]
169    return g0_iw, g0_an_iw
170
171
172#
173# DMFT loop
174#
175
176gooditer = 0
177g0_prev = np.zeros((2, 2 * n_iw, Norb, Norb), dtype=complex)
178for iloop in range(Nloop):
179    print(f"\nLoop {iloop + 1} of {Nloop}")
180
181    # Solve the effective impurity problem
182    solver.solve(beta=beta,
183                 n_iw=n_iw,
184                 energy_window=energy_window,
185                 n_w=n_w,
186                 broadening=broadening)
187
188    # Normal and anomalous components of computed self-energy
189    s_iw = solver.Sigma_iw["up"]
190    s_an_iw = solver.Sigma_an_iw["up_dn"]
191
192    # Compute local Green's function
193    g_iw, g_an_iw = get_gloc(s_iw, s_an_iw)
194    # Compute Weiss field
195    g0_iw, g0_an_iw = dmft_weiss_field(g_iw, g_an_iw, s_iw, s_an_iw)
196
197    # Bath fitting and mixing
198    G0_iw = BlockGf(name_list=spins, block_list=[g0_iw, g0_iw])
199    G0_an_iw = BlockGf(name_list=["up_dn"], block_list=[g0_an_iw])
200
201    bath_new = solver.chi2_fit_bath(G0_iw, G0_an_iw)[0]
202    solver.bath = wmixing * bath_new + (1 - wmixing) * solver.bath
203
204    # Check convergence of the Weiss field
205
206    g0 = np.asarray([g0_iw.data, g0_an_iw.data])
207    errvec = np.real(np.sum(abs(g0 - g0_prev), axis=1)
208                     / np.sum(abs(g0), axis=1))
209    # First iteration
210    if iloop == 0:
211        errvec = np.ones_like(errvec)
212    errmin, err, errmax = np.min(errvec), np.average(errvec), np.max(errvec)
213
214    g0_prev = np.copy(g0)
215
216    if err < threshold:
217        gooditer += 1  # Increase good iterations count
218    else:
219        gooditer = 0  # Reset good iterations count
220    iloop += 1
221    conv_bool = ((err < threshold) and (gooditer > Nsuccess)
222                 and (iloop < Nloop)) or (iloop >= Nloop)
223
224    # Print convergence message
225    if iloop < Nloop:
226        if errvec.size > 1:
227            print(f"max error={errmax:.6e}")
228        print("    " * (errvec.size > 1) + f"error={err:.6e}")
229        if errvec.size > 1:
230            print(f"min error={errmin:.6e}")
231    else:
232        if errvec.size > 1:
233            print(f"max error={errmax:.6e}")
234        print("    " * (errvec.size > 1) + f"error={err:.6e}")
235        if errvec.size > 1:
236            print(f"min error={errmin:.6e}")
237        print(f"Not converged after {Nloop} iterations.")
238
239    if conv_bool:
240        break
241
242
243# Calculate local Green's function on the real axis
244s_w = solver.Sigma_w["up"]
245s_an_w = solver.Sigma_an_w["up_dn"]
246g_w, g_an_w = get_gloc(s_w, s_an_w)
247
248# Save calculation results
249with HDFArchive('ahm.h5', 'w') as ar:
250    ar["s_iw"] = s_iw
251    ar["s_an_iw"] = s_an_iw
252    ar["g_iw"] = g_iw
253    ar["g_an_iw"] = g_an_iw
254    ar["g_w"] = g_w
255    ar["g_an_w"] = g_an_w
256
257print("Done...")