From b1edb1c37799b06e33298832e3f59d9feee3fee3 Mon Sep 17 00:00:00 2001
From: antoniolrm <am@antoniomanesco.org>
Date: Fri, 12 May 2023 15:05:24 +0200
Subject: [PATCH] more packaging and some clean up
---
codes/hf.py | 53 ++++++++++++++++++++++++++++++----
codes/utils.py | 77 +++++++++++++++++++++++++++++++++++++++++---------
2 files changed, 110 insertions(+), 20 deletions(-)
diff --git a/codes/hf.py b/codes/hf.py
index 173f68c..be6d6a9 100644
--- a/codes/hf.py
+++ b/codes/hf.py
@@ -1,9 +1,11 @@
from scipy.signal import convolve2d
import numpy as np
import utils
+from functools import partial
+from scipy.optimize import anderson, minimize
-def mean_field_F(vals, vecs, E_F):
- N_ks = vecs.shape[0]
+
+def mean_field_F(vals, vecs, E_F, nk):
unocc_vals = vals > E_F
def mf_generator(i, j):
@@ -12,9 +14,10 @@ def mean_field_F(vals, vecs, E_F):
F_ij = occ_vecs @ occ_vecs.conj().T
return F_ij
- F = np.array([[mf_generator(i, j) for i in range(N_ks)] for j in range(N_ks)])
+ F = np.array([[mf_generator(i, j) for i in range(nk)] for j in range(nk)])
return F
+
def convolution(M1, M2):
cell_size = M2.shape[-1]
V_output = np.array(
@@ -29,12 +32,50 @@ def convolution(M1, M2):
V_output = np.transpose(V_output, axes=(2, 3, 0, 1))
return V_output
+
def compute_mf(vals, vecs, filling, nk, H_int):
- H0_int = H_int[0,0]
+ H0_int = H_int[0, 0]
E_F = utils.get_fermi_energy(vals, filling)
- F = mean_field_F(vals, vecs, E_F=E_F)
+ F = mean_field_F(vals, vecs, E_F=E_F, nk=nk)
rho = np.diag(np.average(F, axis=(0, 1)))
exchange_mf = convolution(F, H_int) * nk ** (-2)
direct_mf = np.diag(np.einsum("i,ij->j", rho, H0_int))
- return direct_mf - exchange_mf
\ No newline at end of file
+ return direct_mf - exchange_mf
+
+
+def dm(mf0, mf):
+ return np.mean(np.abs(mf - mf0))
+
+
+def scf_loop(mf, H_int, nk, filling, hamiltonians_0, tol):
+ if np.linalg.norm(mf) < tol:
+ return 0
+ # Generate the Hamiltonian
+ hamiltonians = hamiltonians_0 + mf
+ vals, vecs = np.linalg.eigh(hamiltonians)
+ vecs = np.linalg.qr(vecs)[0]
+
+ mf_new = compute_mf(vals=vals, vecs=vecs, filling=filling, nk=nk, H_int=H_int)
+
+ diff = mf_new - mf
+
+ if np.linalg.norm(mf_new) < tol:
+ return 0
+ else:
+ return diff
+
+
+def find_groundstate_ham(
+ H_int, nk, filling, hamiltonians_0, tol, guess, mixing=0.5, order=1
+):
+ fun = partial(
+ scf_loop,
+ H_int=H_int,
+ nk=nk,
+ filling=filling,
+ hamiltonians_0=hamiltonians_0,
+ tol=tol,
+ )
+ mf = anderson(fun, guess, f_tol=tol, w0=mixing, M=order)
+ return hamiltonians_0 + mf
diff --git a/codes/utils.py b/codes/utils.py
index a2fbf84..6d629c8 100644
--- a/codes/utils.py
+++ b/codes/utils.py
@@ -7,12 +7,12 @@ def get_fermi_energy(vals, filling):
ne = len(vals_flat)
ifermi = int(round(ne * filling / norbs))
if ifermi >= ne:
- ifermi = ne - 1
- sorte = np.sort(vals_flat) # sorted eigenvalues
- if ifermi == 0:
- return sorte[0]
- fermi = (sorte[ifermi - 1] + sorte[ifermi]) / 2.0 # fermi energy
- return fermi
+ return vals_flat[-1]
+ elif ifermi == 0:
+ return vals_flat[0]
+ else:
+ fermi = (vals_flat[ifermi - 1] + vals_flat[ifermi]) / 2
+ return fermi
def syst2hamiltonian(kxs, kys, syst, params={}):
def h_k(kx, ky):
@@ -31,22 +31,25 @@ def potential2hamiltonian(
wrapped_V = kwant.wraparound.wraparound(V).finalized()
return syst2hamiltonian(kxs=ks, kys=ks, syst=wrapped_V, params=params)
-def generate_guess(max_neighbor, norbs, lattice):
+def generate_guess(max_neighbor, norbs, lattice, kxs, kys, dummy_syst):
n_sub = len(lattice.sublattices)
guess = np.zeros((n_sub + max_neighbor, 2, norbs, norbs))
for i in range(n_sub):
# Real part
- guess[i, 0] = np.random.rand(norbs, norbs)
+ guess[i, 0] = np.random.rand(norbs, norbs) - 0.5
guess[i, 0] += guess[i, 0].T
# Imag part
- guess[i, 1] = np.random.rand(norbs, norbs)
+ guess[i, 1] = np.random.rand(norbs, norbs) - 0.5
guess[i, 1] -= guess[i, 1].T
for neighbor in range(max_neighbor):
# Real part
- guess[n_sub + neighbor, 0] = np.random.rand(norbs, norbs)
+ guess[n_sub + neighbor, 0] = np.random.rand(norbs, norbs) - 0.5
# Imag part
- guess[n_sub + neighbor, 1] = np.random.rand(norbs, norbs)
- return guess
+ guess[n_sub + neighbor, 1] = np.random.rand(norbs, norbs) - 0.5
+ guess_k = syst2hamiltonian(
+ kxs=kxs, kys=kys, syst=dummy_syst, params=dict(mat=guess)
+ )
+ return guess_k
def extract_hopping_vectors(builder):
keep=None
@@ -69,7 +72,6 @@ def generate_scf_syst(max_neighbor, syst, lattice):
scf = kwant.Builder(kwant.TranslationalSymmetry(*lattice.prim_vecs))
scf[syst.sites()] = scf_onsite
for neighbor in range(max_neighbor):
-
def scf_hopping(site1, site2, mat):
return (
mat[len(lattice.sublattices) + neighbor, 0]
@@ -79,4 +81,51 @@ def generate_scf_syst(max_neighbor, syst, lattice):
scf[lattice.neighbors(neighbor + 1)] = scf_hopping
deltas = extract_hopping_vectors(scf)
wrapped_scf = kwant.wraparound.wraparound(scf).finalized()
- return wrapped_scf, deltas
\ No newline at end of file
+ return wrapped_scf, deltas
+
+def hk2hop(hk, deltas, ks, dk):
+ kxx, kyy = np.meshgrid(ks, ks)
+ kxy = np.array([kxx, kyy])
+ hopps = (
+ np.sum(
+ np.einsum(
+ "ijk,jklm->ijklm",
+ np.exp(1j * np.einsum("ij,jkl->ikl", deltas, kxy)),
+ hk,
+ ),
+ axis=(1, 2),
+ )
+ * (dk) ** 2
+ / (2 * np.pi) ** 2
+ )
+ return hopps
+
+
+def hk2syst(deltas, hk, ks, dk, max_neighbor, norbs, lattice):
+ hopps = hk2hop(hk, deltas, ks, dk)
+ bulk_scf = kwant.Builder(kwant.TranslationalSymmetry(*lattice.prim_vecs))
+ for i, delta in enumerate(deltas):
+ for j, sublattice1 in enumerate(lattice.sublattices):
+ for k, sublattice2 in enumerate(lattice.sublattices):
+ if np.allclose(delta, [0, 0]):
+ bulk_scf[sublattice1.shape((lambda pos: True), (0, 0))] = hopps[
+ i, j * norbs : (j + 1) * norbs, j * norbs : (j + 1) * norbs
+ ]
+ if k != j:
+ hopping = (delta, sublattice1, sublattice2)
+ bulk_scf[kwant.builder.HoppingKind(*hopping)] = hopps[
+ i, j * norbs : (j + 1) * norbs, k * norbs : (k + 1) * norbs
+ ]
+ else:
+ for k, sublattice2 in enumerate(lattice.sublattices):
+ hopping = (delta, sublattice1, sublattice2)
+ bulk_scf[kwant.builder.HoppingKind(*hopping)] = hopps[
+ i, j * norbs : (j + 1) * norbs, k * norbs : (k + 1) * norbs
+ ]
+ wrapped_scf_syst = kwant.wraparound.wraparound(bulk_scf).finalized()
+ return wrapped_scf_syst
+
+def calc_gap(vals, E_F):
+ emax = np.max(vals[vals <= E_F])
+ emin = np.min(vals[vals > E_F])
+ return np.abs(emin - emax)
\ No newline at end of file
--
GitLab