Usage examples

Hubbard model on the Bethe lattice

Input file and converged bath parameters are available in the EDIpack2py repository

 1# INIT MPI
 2comm = MPI.COMM_WORLD
 3rank = comm.Get_rank()
 4print("I am process", rank, "of", comm.Get_size())
 5master = rank == 0
 6
 7
 8def dens_bethe(x, d):
 9    root = (1 - (x / d) ** 2) + 0.0j
10    root = np.sqrt(root)
11    dens_bethe = (2 / (np.pi * d)) * root
12    return dens_bethe.real
13
14
15# READ ED INPUT:
16ed.read_input("inputED.conf")
17if ed.Nspin != 1 or ed.Norb != 1:
18    print("This test code is for Nspin=1 + Norb=1.")
19    ed.Nspin = 1
20    ed.Norb = 1
21Le = 1000
22wmixing = 0.3
23wband = 1.0
24
25# BUILD Density of States:
26Eband, de = np.linspace(-wband, wband, Le, retstep=True)
27Dband = dens_bethe(Eband, wband)
28
29# BUILD frequency array:
30wm = np.pi / ed.beta * (2 * np.arange(ed.Lmats) + 1)
31wr = np.linspace(ed.wini, ed.wfin, ed.Lreal)
32
33
34# Allocate minimal required arrays:
35# ALL functions must have shape [Nspin,Nspin,Norb,Norb(,L)]:
36Delta = np.zeros((ed.Nspin, ed.Nspin, ed.Norb, ed.Norb, ed.Lmats), dtype="complex")
37Hloc = np.zeros((ed.Nspin, ed.Nspin, ed.Norb, ed.Norb), dtype="complex")
38
39
40# SETUP SOLVER
41ed.set_hloc(Hloc)
42Nb = ed.get_bath_dimension()
43bath = ed.init_solver()
44bath_prev = np.copy(bath)
45
46# DMFT CYCLE
47converged = False
48iloop = 0
49while not converged and iloop < ed.Nloop:
50    iloop = iloop + 1
51    print("DMFT-loop:", iloop, "/", ed.Nloop)
52
53    # Solve the EFFECTIVE IMPURITY PROBLEM (first w/ a guess for the bath)
54    ed.solve(bath)
55
56    # Get Self-energies
57    Smats = ed.get_sigma(axis="m")
58    Sreal = ed.get_sigma(axis="r")
59    Gimp_mats = ed.get_gimp(axis="m")
60    Gimp_real = ed.get_gimp(axis="r")
61    Gmats = np.zeros_like(Smats)
62    Greal = np.zeros_like(Sreal)
63
64    zeta = wm * 1j - Smats[0, 0, 0, 0, :]
65    Gmats[0, 0, 0, 0, :] = np.sum(Dband / (zeta[:, np.newaxis] - Eband), axis=1) * de
66
67    zeta = wr + ed.eps * 1j - Sreal[0, 0, 0, 0, :]
68    Greal[0, 0, 0, 0, :] = np.sum(Dband / (zeta[:, np.newaxis] - Eband), axis=1) * de
69
70    if rank == 0:
71        np.savetxt(
72            "Gloc_iw.dat",
73            np.transpose([wm, Gmats[0, 0, 0, 0, :].imag, Gmats[0, 0, 0, 0, :].real]),
74        )
75        np.savetxt(
76            "Gloc_realw.dat",
77            np.transpose([wr, Greal[0, 0, 0, 0, :].imag, Greal[0, 0, 0, 0, :].real]),
78        )
79
80    # Get the Delta function and FIT:
81    Delta[0, 0, 0, 0, :] = 0.25 * wband * Gmats[0, 0, 0, 0, :]
82    if rank == 0:
83        np.savetxt(
84            "Delta_iw.dat",
85            np.transpose([wm, Delta[0, 0, 0, 0, :].imag, Delta[0, 0, 0, 0, :].real]),
86        )
87    bath = ed.chi2_fitgf(Delta, bath, ispin=0, iorb=0)
88
89    if iloop > 1:
90        bath = wmixing * bath + (1.0 - wmixing) * bath_prev
91    bath_prev = np.copy(bath)
92
93    err, converged = ed.check_convergence(Delta, ed.dmft_error)
94
95ed.finalize_solver()
96print("Done...")

2-orbital attractive Hubbard model on the square lattice

Input file is available in the EDIpack2py repository

  1# INIT MPI
  2comm = MPI.COMM_WORLD
  3rank = comm.Get_rank()
  4print("I am process", rank, "of", comm.Get_size())
  5master = rank == 0
  6
  7
  8# Auxiliary functions
  9def superconductive_zeta(warray, xmu, Sigma_all, axis):
 10    Ntot = np.shape(Sigma_all)[2]
 11    Lfreq = np.shape(warray)[0]
 12    zi = np.zeros(
 13        (2, 2, Ntot, Ntot, Lfreq), dtype=complex
 14    )
 15
 16    if axis == "m":
 17        zi[0, 0, :, :, :] = (warray + xmu) * np.eye(Ntot)[..., None] - Sigma_all[0]
 18        zi[0, 1, :, :, :] = -Sigma_all[1]
 19        zi[1, 0, :, :, :] = -Sigma_all[1]
 20        zi[1, 1, :, :, :] = (warray - xmu) * np.eye(Ntot)[..., None] + np.conj(
 21            Sigma_all[0]
 22        )
 23    elif axis == "r":
 24        warray_bar = warray[::-1, ...]
 25        Sigma_bar = Sigma_all[0, ..., ::-1]
 26        zi[0, 0, :, :, :] = (warray + xmu) * np.eye(Ntot)[..., None] - Sigma_all[0]
 27        zi[0, 1, :, :, :] = -Sigma_all[1]
 28        zi[1, 0, :, :, :] = -Sigma_all[1]
 29        zi[1, 1, :, :, :] = -np.conj(warray_bar + xmu) * np.eye(Ntot)[
 30            ..., None
 31        ] + np.conj(Sigma_bar)
 32    return zi
 33
 34
 35def get_gloc(warray, xmu, Hk, Sigma_all, axis):
 36    """
 37    Z has dimension  [Nambu,Nambu,Nso,Nso,Nfreq]
 38    Hk has dimension [Nnambu,Nk,Nso,Nso]
 39    Gk has dimension [Nk,Nnambu,Nso,Nso]
 40    Gmatrix has dimension [Nnambu*Ntot,Nnambu*Ntot,Nfreq]
 41    Returns an object of dimension [2,Ntot,Nfreq]
 42    """
 43
 44    try:
 45        import mpi4py
 46        from mpi4py import MPI
 47
 48        comm = MPI.COMM_WORLD
 49        rank = comm.Get_rank()
 50        size = comm.Get_size()
 51        mpiflag = True
 52    except:
 53        mpiflag = False
 54        rank = 0
 55        size = 1
 56
 57    master = rank == 0
 58
 59    if master:
 60        print("Calculating local G axis " + axis + ":")
 61
 62    Z = superconductive_zeta(warray, xmu, Sigma_all, axis)
 63    Ntot = np.shape(Sigma_all)[2]
 64    Nfreq = np.shape(Z)[-1]
 65    Nk = np.shape(Hk)[1]
 66
 67    if Nk >= Nfreq:
 68        base = int(Nk // size)
 69        leftover = int(Nk % size)
 70        chunks = np.ones(size, dtype=int) * base
 71        chunks[:leftover] += 1
 72        offsets = np.zeros(size, dtype=int)
 73        offsets[1:] = np.cumsum(chunks)[:-1]
 74        ilow = offsets[rank]
 75        ihigh = ilow + chunks[rank]
 76        # print(rank,ilow,ihigh,chunks[rank])
 77
 78        Gtmp = np.zeros((chunks[rank], 2 * Ntot, 2 * Ntot, Nfreq), dtype=complex)
 79
 80        Gtmp[:, 0:Ntot, 0:Ntot, :] = Z[0, 0, :, :, :] - Hk[0, ilow:ihigh, :, :, None]
 81        Gtmp[:, 0:Ntot, Ntot : 2 * Ntot, :] = Z[0, 1, :, :, :]
 82        Gtmp[:, Ntot : 2 * Ntot, 0:Ntot, :] = Z[1, 0, :, :, :]
 83        Gtmp[:, Ntot : 2 * Ntot, Ntot : 2 * Ntot, :] = (
 84            Z[1, 1, :, :, :] - Hk[1, ilow:ihigh, :, :, None]
 85        )
 86
 87        Gtmp = Gtmp.transpose(0, 3, 1, 2)
 88        Gtmp = np.linalg.inv(Gtmp)
 89        Gtmp = Gtmp.transpose(0, 2, 3, 1)
 90        if mpiflag:
 91            Gtmp = np.ascontiguousarray(np.sum(Gtmp, axis=0) / Nk)
 92            Gloc = np.zeros_like(Gtmp)
 93            comm.Allreduce(Gtmp, Gloc, op=MPI.SUM)
 94        else:
 95            Gloc = np.sum(Gtmp, axis=0) / Nk
 96    else:
 97        base = int(Nfreq // size)
 98        leftover = int(Nfreq % size)
 99        chunks = np.ones(size, dtype=int) * base
100        chunks[:leftover] += 1
101        offsets = np.zeros(size, dtype=int)
102        offsets[1:] = np.cumsum(chunks)[:-1]
103        ilow = offsets[rank]
104        ihigh = ilow + chunks[rank]
105        # print(rank,ilow,ihigh,chunks[rank])
106
107        Gtmp = np.zeros((Nk, 2 * Ntot, 2 * Ntot, Nfreq), dtype=complex)
108        Gloc = np.zeros((2 * Ntot, 2 * Ntot, Nfreq), dtype=complex)
109
110        Gtmp[:, 0:Ntot, 0:Ntot, ilow:ihigh] = (
111            Z[0, 0, :, :, ilow:ihigh] - Hk[0, :, :, :, None]
112        )
113        Gtmp[:, 0:Ntot, Ntot : 2 * Ntot, ilow:ihigh] = Z[0, 1, :, :, ilow:ihigh]
114        Gtmp[:, Ntot : 2 * Ntot, 0:Ntot, ilow:ihigh] = Z[1, 0, :, :, ilow:ihigh]
115        Gtmp[:, Ntot : 2 * Ntot, Ntot : 2 * Ntot, ilow:ihigh] = (
116            Z[1, 1, :, :, ilow:ihigh] - Hk[1, :, :, :, None]
117        )
118
119        Gtmp = Gtmp.transpose(0, 3, 1, 2)
120        Gtmp[:, ilow:ihigh, :, :] = np.linalg.inv(Gtmp[:, ilow:ihigh, :, :])
121        Gtmp = Gtmp.transpose(0, 2, 3, 1)
122
123        if mpiflag:
124            Gtmp = np.ascontiguousarray(np.sum(Gtmp, axis=0) / Nk)
125            comm.Allreduce(Gtmp, Gloc, op=MPI.SUM)
126        else:
127            Gloc = np.sum(Gtmp, axis=0) / Nk
128
129    return np.stack((Gloc[:Ntot, :Ntot, :], Gloc[:Ntot, Ntot:, :]), axis=0)
130
131
132def get_weiss_field(G, Sigma):
133    """
134    complex(8),dimension(2,Ntot,Ntot,Nfreq)              :: G
135    complex(8),dimension(2,Ntot,Ntot,Nfreq)              :: Sigma
136    """
137    try:
138        import mpi4py
139        from mpi4py import MPI
140
141        comm = MPI.COMM_WORLD
142        rank = comm.Get_rank()
143        size = comm.Get_size()
144        mpiflag = True
145    except:
146        mpiflag = False
147        rank = 0
148        size = 1
149
150    master = rank == 0
151
152    if master:
153        print("Calculating Weiss Field")
154
155    Ntot = np.shape(G)[1]
156    Nfreq = np.shape(G)[-1]
157
158    Gnambu = np.zeros((2 * Ntot, 2 * Ntot, Nfreq), dtype=complex)
159    Snambu = np.zeros((2 * Ntot, 2 * Ntot, Nfreq), dtype=complex)
160
161    Snambu[:Ntot, :Ntot, :] = Sigma[0, :, :, :]
162    Snambu[:Ntot, Ntot:, :] = Sigma[1, :, :, :]
163    Snambu[Ntot:, :Ntot, :] = Sigma[1, :, :, :]
164    Snambu[Ntot:, Ntot:, :] = -np.conj(Sigma[0, :, :, :])
165
166    Gnambu[:Ntot, :Ntot, :] = G[0, :, :, :]
167    Gnambu[:Ntot, Ntot:, :] = G[1, :, :, :]
168    Gnambu[Ntot:, :Ntot, :] = G[1, :, :, :]
169    Gnambu[Ntot:, Ntot:, :] = -np.conj(G[0, :, :, :])
170
171    Weiss = np.linalg.inv(
172        np.linalg.inv(Gnambu.transpose(2, 0, 1)) + Snambu.transpose(2, 0, 1)
173    ).transpose(1, 2, 0)
174
175    return np.stack((Weiss[:Ntot, :Ntot, :], Weiss[:Ntot, Ntot:, :]), axis=0)
176
177
178def generate_kgrid(Nk):
179    b1 = 2 * np.pi * np.array([1.0, 0.0])
180    b2 = 2 * np.pi * np.array([0.0, 1.0])
181    n1, n2 = np.meshgrid(np.arange(Nk), np.arange(Nk))
182    n1 = n1 / Nk
183    n2 = n2 / Nk
184    gridout = np.stack([n1.ravel(), n2.ravel()], axis=-1)
185    return np.dot(gridout, [b1, b2])
186
187
188def h_square2d(k, t):
189    return (
190        -2
191        * t
192        * (
193            np.cos(k[..., 0, np.newaxis, np.newaxis])
194            + np.cos(k[..., 1, np.newaxis, np.newaxis])
195        )
196        * np.eye(ed.Norb)
197    )
198
199
200def test_dos(hk, plot=False):
201    result = np.histogram(
202        hk[:, 0, 0], bins=100, density=True
203    )  # normalized to one, I will have to multiply
204    dx = result[1][2] - result[1][1]
205    np.savetxt(
206        "dos.dat",
207        np.transpose(
208            [np.delete(np.add(result[1], dx / 2), np.size(result[1]) - 1), result[0]]
209        ),
210    )
211
212    if plot:
213        plt.xlabel("E")
214        plt.ylabel("D(E)")
215        plt.xlim(-3, 3)
216        plt.ylim(0, 1)
217        plt.plot(
218            np.delete(np.add(result[1], dx / 2), np.size(result[1]) - 1), result[0]
219        )
220        plt.show()
221
222
223# READ ED INPUT:
224ed.read_input("inputAHM.conf")
225
226
227# Parameters
228wmixing = 0.5
229try:
230    Nk = int(sys.argv[1])
231except:
232    Nk = 30
233
234print(ed.Uloc)
235t_hop = 0.25
236
237# BUILD frequency arrays and k grid:
238wm = np.pi / ed.beta * (2 * np.arange(ed.Lmats) + 1)
239wr = np.linspace(ed.wini, ed.wfin, ed.Lreal, dtype=complex)
240
241kgrid = generate_kgrid(Nk)
242
243
244# Generate hk and hloc
245Hk = h_square2d(kgrid, t_hop)
246HkNambu = np.array([h_square2d(kgrid, t_hop), -np.conj(h_square2d(-kgrid, t_hop))])
247Hloc = np.sum(Hk, axis=0) / Nk**2
248Hloc = Hloc.astype(complex)
249
250
251# Generate dos and plot it
252test_dos(Hk)
253
254# SETUP SOLVER
255ed.set_hloc(Hloc)
256Nb = ed.get_bath_dimension()
257bath = ed.init_solver()
258bath_prev = np.copy(bath)
259
260
261# DMFT CYCLE
262converged = False
263iloop = 0
264while not converged and iloop < ed.Nloop:
265    iloop = iloop + 1
266    print("DMFT-loop:", iloop, "/", ed.Nloop)
267
268    ed.solve(bath)
269
270    Smats = np.array([ed.get_sigma(axis="m", typ="n"), ed.get_sigma(axis="m", typ="a")])
271    Sreal = np.array(
272        [
273            ed.build_sigma(wr + 1j * ed.eps, typ="n"),
274            ed.build_sigma(wr + 1j * ed.eps, typ="a"),
275        ]
276    )
277
278    Gmats = get_gloc(wm * 1j, ed.xmu, HkNambu, Smats, axis="m")
279    Greal = get_gloc(wr + 1j * ed.eps, ed.xmu, HkNambu, Sreal, axis="r")
280    Weiss = get_weiss_field(Gmats, Smats)
281
282    # Print
283
284    if rank == 0:
285        for ispin in range(ed.Nspin):
286            for iorb in range(ed.Norb):
287                io = iorb + ed.Norb * ispin
288                for typ in ["G", "F"]:
289                    for axis in ["i", "real"]:
290                        name = (
291                            typ
292                            + "_l"
293                            + str(iorb + 1)
294                            + "_s"
295                            + str(ispin + 1)
296                            + "_"
297                            + axis
298                            + "w.dat"
299                        )
300                        comp = 0 if typ == "G" else 1
301                        Grfx = Gmats if axis == "i" else Greal
302                        freq = wm if axis == "i" else np.real(wr)
303                        np.savetxt(
304                            name,
305                            np.transpose(
306                                [
307                                    freq,
308                                    Grfx[comp, io, io, :].imag,
309                                    Grfx[comp, io, io, :].real,
310                                ]
311                            ),
312                        )
313
314    # Fit
315    bath = ed.chi2_fitgf(Weiss[0], Weiss[1], bath)
316
317    if iloop > 1:
318        bath = wmixing * bath + (1.0 - wmixing) * bath_prev
319    bath_prev = np.copy(bath)
320
321    err, converged = ed.check_convergence(Weiss, ed.dmft_error)
322
323ed.finalize_solver()
324print("Done...")