Skip to content
Snippets Groups Projects

Interface refactoring

Merged Interface refactoring
interface-refactoring into main
Merged Kostas Vilkelis requested to merge interface-refactoring into main
Compare
and
Show latest version
1 file
+ 1
− 1
Compare changes
  • Side-by-side
  • Inline
examples/testInterface.ipynb
+ 1
− 1
%% Cell type:code id: tags:
``` python
import numpy as np
from scipy.integrate import quad, quad_vec
```
%% Cell type:code id: tags:
``` python
def tb2kfunc(tb_model):
"""
Fourier transforms a real-space tight-binding model to a k-space function.
Parameters
----------
tb_model : dict
A dictionary with real-space vectors as keys and complex np.arrays as values.
Returns
-------
function
A function that takes a k-space vector and returns a complex np.array.
"""
def bloch_ham(k):
ham = 0
for vector in tb_model.keys():
ham += tb_model[vector] * np.exp(
-1j * np.dot(k, np.array(vector, dtype=float))
)
return ham
return bloch_ham
def kfunc2tbfunc(kfunc):
"""
Inverse Fourier transforms a k-space function to a real-space function.
Parameters
----------
kfunc : function
A function that takes a k-space vector and returns a np.array.
Returns
-------
function
A function that takes a real-space integer vector and returns a np.array.
"""
def tbfunc(vector):
def integrand(k):
return kfunc(k) * np.exp(1j * np.dot(k, np.array(vector, dtype=float)))
return quad_vec(integrand, -np.pi, np.pi)[0]
return tbfunc
```
%% Cell type:code id: tags:
``` python
def addTb(tb1, tb2):
"""
Add up two tight-binding models together.
Parameters:
-----------
tb1 : dict
Tight-binding model.
tb2 : dict
Tight-binding model.
Returns:
--------
dict
Sum of the two tight-binding models.
"""
return {k: tb1.get(k, 0) + tb2.get(k, 0) for k in frozenset(tb1) | frozenset(tb2)}
def generateRho(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
-------
rho : function
Density matrix at the same k-point as hk
"""
def rhofunc(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 rhofunc
def fermiOnGrid(hkfunc, filling, Nk=100):
"""
Compute the Fermi energy on a grid of k-points.
Parameters
----------
hkfunc : function
Function that returns the Hamiltonian at a given k-point.
Nk : int
Number of k-points in the grid.
Returns
-------
E_F : float
Fermi energy
"""
ks = np.linspace(-np.pi, np.pi, Nk, endpoint=False)
hkarray = np.array([hkfunc(k) for k in ks])
vals = np.linalg.eigvalsh(hkarray)
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 totalEnergy(rho, hk):
def integrand(k):
return np.real(np.trace(rho(k) @ hk(k)))
return quad(integrand, -np.pi, np.pi)[0]
def meanField(rhoTbfunc, int_model):
"""
Compute the mean-field in real space.
Parameters
----------
rhoTbfunc : function
Function that returns the density matrix at a given real-space vector.
int_model : dict
Interaction tb model.
Returns
-------
dict
Mean-field tb model.
"""
exchange = {vec: -1*(int_model.get(vec, 0) + rhoTbfunc(vec)) for vec in frozenset(int_model)}
exchange = {vec: -1*(int_model.get(vec, 0)*rhoTbfunc(vec)) for vec in frozenset(int_model)}
# direct interaction is missing
return exchange
```
%% Cell type:code id: tags:
``` python
U0 = 1
nk = 100
filling = 2
hopp = np.kron(np.array([[0, 1], [0, 0]]), np.eye(2))
tb_model = {(0,): hopp + hopp.T.conj(), (1,): hopp, (-1,): hopp.T.conj()}
mf_model = {} # placeholder for guess
int_model = {
(0,): U0 * np.kron(np.eye(2), np.ones((2, 2))),
}
# attempting to do a single iteration of mean-field
hkfunc = tb2kfunc(addTb(tb_model, mf_model))
EF = fermiOnGrid(hkfunc, filling, nk)
rhofunc = generateRho(hkfunc, EF)
rhoTbfunc = kfunc2tbfunc(rhofunc)
mf_model = meanField(rhoTbfunc, int_model)
```