From 3a22cd6df30d4bf0ea53b95068c1878b9d21c244 Mon Sep 17 00:00:00 2001
From: Johanna <johanna@zijderveld.de>
Date: Thu, 11 Apr 2024 22:41:36 +0200
Subject: [PATCH] remove old functions and remove kvector ending on new
functions
---
codes/mf.py | 96 ++++--------------------------------------
codes/model.py | 37 ++++------------
codes/solvers.py | 60 +-------------------------
codes/test_graphene.py | 6 +--
4 files changed, 19 insertions(+), 180 deletions(-)
diff --git a/codes/mf.py b/codes/mf.py
index 4b81464..6ae70ed 100644
--- a/codes/mf.py
+++ b/codes/mf.py
@@ -2,36 +2,6 @@ import numpy as np
from codes.tb.tb import addTb
-def densityMatrixGenerator(hkfunc, E_F):
- """
- Generate a function that returns the density matrix at a given k-point.
-
- Parameters
- ----------
- hkfunc : function
- Function that return Hamiltonian at a given k-point.
- E_F : float
- Fermi level
-
- Returns
- -------
- densityMatrixFunc : function
- Returns a density matrix at a given k-point (kx, kx, ...)
- """
-
- def densityMatrixFunc(k):
- hk = hkfunc(k)
- vals, vecs = np.linalg.eigh(hk)
- unocc_vals = vals > E_F
- occ_vecs = vecs
- occ_vecs[:, unocc_vals] = 0
-
- # Outter products between eigenvectors
- return occ_vecs @ occ_vecs.T.conj()
-
- return densityMatrixFunc
-
-
def densityMatrix(kham, E_F):
"""
Parameters
@@ -56,24 +26,7 @@ def densityMatrix(kham, E_F):
return densityMatrixKgrid
-def fermiOnGridkvector(kham, filling):
-
- vals = np.linalg.eigvalsh(kham)
-
- norbs = vals.shape[-1]
- vals_flat = np.sort(vals.flatten())
- ne = len(vals_flat)
- ifermi = int(round(ne * filling / norbs))
- if ifermi >= ne:
- return vals_flat[-1]
- elif ifermi == 0:
- return vals_flat[0]
- else:
- fermi = (vals_flat[ifermi - 1] + vals_flat[ifermi]) / 2
- return fermi
-
-
-def fermiOnGrid(hkfunc, filling, nK=100, ndim=1): # need to extend to 2D
+def fermiOnGrid(kham, filling):
"""
Compute the Fermi energy on a grid of k-points.
@@ -88,17 +41,8 @@ def fermiOnGrid(hkfunc, filling, nK=100, ndim=1): # need to extend to 2D
E_F : float
Fermi energy
"""
- ks = np.linspace(-np.pi, np.pi, nK, endpoint=False)
- if ndim == 1:
- hkarray = np.array([hkfunc(k) for k in ks])
- if ndim == 2:
- hkarray = np.array([[hkfunc((kx, ky)) for kx in ks] for ky in ks])
- elif ndim > 2:
- raise NotImplementedError(
- "Fermi energy calculation is not implemented for ndim > 2"
- )
- vals = np.linalg.eigvalsh(hkarray)
+ vals = np.linalg.eigvalsh(kham)
norbs = vals.shape[-1]
vals_flat = np.sort(vals.flatten())
@@ -113,32 +57,7 @@ def fermiOnGrid(hkfunc, filling, nK=100, ndim=1): # need to extend to 2D
return fermi
-def meanFieldFFTkvector(densityMatrixTb, h_int, n=2):
-
- localKey = tuple(np.zeros((n,), dtype=int))
-
- direct = {
- localKey: np.sum(
- np.array(
- [
- np.diag(
- np.einsum("pp,pn->n", densityMatrixTb[localKey], h_int[vec])
- )
- for vec in frozenset(h_int)
- ]
- ),
- axis=0,
- )
- }
-
- exchange = {
- vec: -1 * h_int.get(vec, 0) * densityMatrixTb[vec] # / (2 * np.pi)#**2
- for vec in frozenset(h_int)
- }
- return addTb(direct, exchange)
-
-
-def meanField(densityMatrix, int_model, n=2, nK=100):
+def meanField(densityMatrixTb, h_int, n=2):
"""
Compute the mean-field in k-space.
@@ -154,6 +73,7 @@ def meanField(densityMatrix, int_model, n=2, nK=100):
dict
Mean-field tb model.
"""
+
localKey = tuple(np.zeros((n,), dtype=int))
direct = {
@@ -161,9 +81,9 @@ def meanField(densityMatrix, int_model, n=2, nK=100):
np.array(
[
np.diag(
- np.einsum("pp,pn->n", densityMatrix[localKey], int_model[vec])
+ np.einsum("pp,pn->n", densityMatrixTb[localKey], h_int[vec])
)
- for vec in frozenset(int_model)
+ for vec in frozenset(h_int)
]
),
axis=0,
@@ -171,7 +91,7 @@ def meanField(densityMatrix, int_model, n=2, nK=100):
}
exchange = {
- vec: -1 * int_model.get(vec, 0) * densityMatrix[vec] # / (2 * np.pi)#**2
- for vec in frozenset(int_model)
+ vec: -1 * h_int.get(vec, 0) * densityMatrixTb[vec] # / (2 * np.pi)#**2
+ for vec in frozenset(h_int)
}
return addTb(direct, exchange)
diff --git a/codes/model.py b/codes/model.py
index de0b5c1..cad2a5b 100644
--- a/codes/model.py
+++ b/codes/model.py
@@ -1,13 +1,10 @@
# %%
from codes.tb.tb import addTb
-from codes.tb.transforms import tb2kfunc, tb2kham, kdens2tbFFT, kfunc2tb, ifftn2tb
+from codes.tb.transforms import tb2kham, kdens2tbFFT, ifftn2tb
from codes.mf import (
- densityMatrixGenerator,
densityMatrix,
- fermiOnGridkvector,
- meanFieldFFTkvector,
- meanField,
fermiOnGrid,
+ meanField,
)
import numpy as np
@@ -34,37 +31,19 @@ class Model:
_check_hermiticity(h_0)
_check_hermiticity(h_int)
- def calculateEF(self, nK=200):
- self.EF = fermiOnGrid(self.hkfunc, self.filling, nK=nK, ndim=self._ndim)
+ def calculateEF(self):
+ self.EF = fermiOnGrid(self.kham, self.filling)
def makeDensityMatrix(self, mf_model, nK=200):
- self.hkfunc = tb2kfunc(addTb(self.h_0, mf_model))
- self.calculateEF(nK=nK)
- return kfunc2tb(
- densityMatrixGenerator(self.hkfunc, self.EF), nSamples=nK, ndim=self._ndim
- )
-
- def mfield(self, mf_model, nK=200):
- self.densityMatrix = self.makeDensityMatrix(mf_model, nK=nK)
- return addTb(
- meanField(self.densityMatrix, self.h_int, n=self._ndim, nK=nK),
- {self._localKey: -self.EF * np.eye(self._size)},
- )
-
- #######################
- def calculateEFkvector(self):
- self.EF = fermiOnGridkvector(self.kham, self.filling)
-
- def makeDensityMatrixkvector(self, mf_model, nK=200):
self.kham = tb2kham(addTb(self.h_0, mf_model), nK=nK, ndim=self._ndim)
- self.calculateEFkvector()
+ self.calculateEF()
return densityMatrix(self.kham, self.EF)
- def mfieldFFTkvector(self, mf_model, nK=200):
- densityMatrix = self.makeDensityMatrixkvector(mf_model, nK=nK)
+ def mfield(self, mf_model, nK=200):
+ densityMatrix = self.makeDensityMatrix(mf_model, nK=nK)
densityMatrixTb = ifftn2tb(kdens2tbFFT(densityMatrix, self._ndim))
return addTb(
- meanFieldFFTkvector(densityMatrixTb, self.h_int, n=self._ndim),
+ meanField(densityMatrixTb, self.h_int, n=self._ndim),
{self._localKey: -self.EF * np.eye(self._size)},
)
diff --git a/codes/solvers.py b/codes/solvers.py
index 9de9afb..c0acf56 100644
--- a/codes/solvers.py
+++ b/codes/solvers.py
@@ -27,28 +27,6 @@ def cost(mf_param, Model, nK=100):
return mf_params_new - mf_param
-def costkvector(mf_param, Model, nK=100):
- """
- Define the cost function for fixed point iteration.
- The cost function is the difference between the input mean-field real space parametrisation
- and a new mean-field.
-
- Parameters
- ----------
- mf_param : numpy.array
- The mean-field real space parametrisation.
- Model : Model
- The model object.
- nK : int, optional
- The number of k-points to use in the grid. The default is 100.
- """
- shape = Model._size
- mf_tb = rParams2mf(mf_param, list(Model.h_int), shape)
- mf_tb_new = Model.mfieldFFTkvector(mf_tb, nK=nK)
- mf_params_new = mf2rParams(mf_tb_new)
- return mf_params_new - mf_param
-
-
def solver(
Model, mf_guess, nK=100, optimizer=scipy.optimize.anderson, optimizer_kwargs={}
):
@@ -80,42 +58,6 @@ def solver(
result = rParams2mf(
optimizer(f, mf_params, **optimizer_kwargs), list(Model.h_int), shape
)
- Model.calculateEF(nK=nK)
- localKey = tuple(np.zeros((Model._ndim,), dtype=int))
- return addTb(result, {localKey: -Model.EF * np.eye(shape)})
-
-
-def solverkvector(
- Model, mf_guess, nK=100, optimizer=scipy.optimize.anderson, optimizer_kwargs={}
-):
- """
- Solve the mean-field self-consistent equation.
-
- Parameters
- ----------
- Model : Model
- The model object.
- mf_guess : numpy.array
- The initial guess for the mean-field tight-binding model.
- nK : int, optional
- The number of k-points to use in the grid. The default is 100.
- optimizer : scipy.optimize, optional
- The optimizer to use to solve for fixed-points. The default is scipy.optimize.anderson.
- optimizer_kwargs : dict, optional
- The keyword arguments to pass to the optimizer. The default is {}.
-
- Returns
- -------
- result : numpy.array
- The mean-field tight-binding model.
- """
-
- shape = Model._size
- mf_params = mf2rParams(mf_guess)
- f = partial(costkvector, Model=Model, nK=nK)
- result = rParams2mf(
- optimizer(f, mf_params, **optimizer_kwargs), list(Model.h_int), shape
- )
- Model.calculateEFkvector()
+ Model.calculateEF()
localKey = tuple(np.zeros((Model._ndim,), dtype=int))
return addTb(result, {localKey: -Model.EF * np.eye(shape)})
diff --git a/codes/test_graphene.py b/codes/test_graphene.py
index dc6da60..9dbf414 100644
--- a/codes/test_graphene.py
+++ b/codes/test_graphene.py
@@ -1,7 +1,7 @@
# %%
import numpy as np
from codes.model import Model
-from codes.solvers import solverkvector
+from codes.solvers import solver
from codes import kwant_examples
from codes.kwant_helper import utils
from codes.tb.utils import compute_gap
@@ -49,9 +49,7 @@ def gap_prediction(U, V):
guess = utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
model = Model(h_0, h_int, filling)
- mf_sol = solverkvector(
- model, guess, nK=nK, optimizer_kwargs={"verbose": True, "M": 0}
- )
+ mf_sol = solver(model, guess, nK=nK, optimizer_kwargs={"verbose": True, "M": 0})
gap = compute_gap(addTb(h_0, mf_sol), n=100)
# Check if the gap is predicted correctly
--
GitLab