import numpy as np
import tbmodels as tbm
import pythtb as ptb
import wannierberri as wberri
from . import _tbmodels
from . import _pythtb
from . import _wberri
import scipy.linalg as la
[docs]
class Model():
r"""
Generic :python:`Model` class from which :python:`Supercell` and :python:`FiniteModel` classes are constructed. Contains the functions which are needed to load spin, positions, Hamiltonian matrices, and other useful information.
Parameters
----------
tbmodel :
Model instance.
spinful :
Whether the model is spinful or not. Default is :python:`False`.
states_uc :
Number of states per unit cell (if the model is spinful, this should be twice the dimension of the basis). Default is :python:`None`, which means it is computed from :python:`tbmodel`.
Lx :
Number of unit cells repeated along the :math:`\mathbf{a}_1` direction. Default is ``1``.
Ly :
Number of unit cells repeated along the :math:`\mathbf{a}_2` direction. Default is ``1``.
file_spn :
Whether the file ``.spn`` is provided if the model is obtained from Wannier functions generated by Wannier90 code starting from a :math:`\Gamma`-only *ab initio* calculation. Default is :python:`False`.
n_occ :
Number of occupied states. This information is required *only* if the model is obtained from Wannier functions generated by Wannier90 code starting from a :math:`\Gamma`-only *ab initio* calculation. Otherwise, half-filling is automatically assumed.
.. warning::
This class should **not** be called explicitly by the user since it is already called from both :python:`Supercell` and :python:`FiniteModel` upon creation.
"""
def __init__(self, tbmodel = None, spinful : bool = False, states_uc : int = None, Lx : int = 1, Ly : int = 1,
file_spn : bool = False, n_occ : int = None):
# Store local variables
self.model = tbmodel
self.Lx = Lx
self.Ly = Ly
self.spinful = spinful
# Wannier tight-binding model variables
if isinstance(self.model, wberri.System_w90):
self.wannier = True
self.file_spn = file_spn
if n_occ == None:
raise NotImplementedError("Provide number of occupied states as n_occ.")
else:
self.wannier = False
self.file_spn = False
# Generate the positions of the atoms in the lattice
self.cart_positions = self._get_positions()
self.n_orb = self.cart_positions[:,0].size
self.dim_r = self._get_dim()
if self.wannier:
self.hamiltonian, self.data = self._get_hamiltonian()
self.n_occ = n_occ
else:
self.hamiltonian, self.n_occ = self._get_hamiltonian()
# Generate the spin matrices if the model is spinful
if self.spinful:
if self.wannier and self.file_spn:
self.eig, self.un_0 = self._diagonalize()
self.sz = self._read_spn()
else:
self.sz = self._get_spin()
else:
self.sz = None
# Generate the position matrices
self.r = []
for d in range (self.dim_r):
self.r.append( np.diag(self.cart_positions[:,d]) )
# Number of states per unit cell
if states_uc == None:
self.states_uc = self._calc_states_uc(self.model, self.spinful)
else:
self.states_uc = states_uc
# Disordered system
self.disordered = False
self._mask = self._initialize_mask()
#################################################
# Load functions
#################################################
def _get_positions(self):
r"""
Returns the cartesian coordinates of the states in the sample.
"""
if isinstance(self.model, tbm.Model):
return _tbmodels.get_positions(self.model)
elif isinstance(self.model, ptb.tb_model):
return _pythtb.get_positions(self.model, self.spinful)
elif isinstance(self.model, wberri.System_w90):
return _wberri.get_positions(self.model)
else:
raise NotImplementedError("Invalid model instance.")
def _get_hamiltonian(self):
r"""
Returns the hamiltonian matrix at the :math:`\Gamma`-point.
"""
if isinstance(self.model, tbm.Model):
gamma_point = np.zeros(self.dim_r)
return _tbmodels.get_hamiltonian(self.model, gamma_point)
elif isinstance(self.model, ptb.tb_model):
gamma_point = np.zeros(self.dim_r)
return _pythtb.get_hamiltonian(self.model, self.spinful, gamma_point, self.dim_r)
elif isinstance(self.model, wberri.System_w90):
return _wberri.get_hamiltonian(self.model)
else:
raise NotImplementedError("Invalid model instance.")
def _diagonalize(self):
"""
Returns the eigenvalues and the eigenstates of the Hamiltonian matrix at the :math:`Gamma`-point along rows.
"""
eig, u_n0 = la.eigh(self.hamiltonian)
u_n0 = u_n0.T
return eig, u_n0
def _get_spin(self):
r"""
Returns the spin matrix elements in the basis of tight-binding orbitals which are diagonal in the spin operator :math:`S_z`.
"""
if isinstance(self.model, tbm.Model) or isinstance(self.model, ptb.tb_model) or isinstance(self.model, wberri.System_w90):
return np.diag([1, -1] * (self.n_orb//2))
else:
raise NotImplementedError("Invalid model instance.")
def _read_spn(self):
"""
Returns the spin matrix elements in the basis of Wannier functions
"""
if isinstance(self.model, wberri.System_w90):
return _wberri.read_spn(self.model,self.data,self.un_0)
else:
raise NotImplementedError("Invalid model instance.")
def _get_dim(self):
r"""
Returns the real-space dimensionality of the model.
Parameters
----------
model :
Model instance.
"""
if isinstance(self.model, tbm.Model):
return self.model.dim
elif isinstance(self.model, ptb.tb_model):
return self.model._dim_r
elif isinstance(self.model, wberri.System_w90):
return(self.cart_positions[0,:].size)
else:
raise NotImplementedError("Invalid model instance.")
def _calc_states_uc(self, model, spinful):
r"""
Returns number of states per unit cell.
Parameters
----------
model :
Model instance.
spinful :
Whether the model is spinful or not.
"""
if isinstance(model, tbm.Model):
return _tbmodels.calc_states_uc(model)
elif isinstance(model, ptb.tb_model):
return _pythtb.calc_states_uc(model, spinful)
elif isinstance(model, wberri.System_w90):
return _wberri.calc_states_uc(model)
else:
raise NotImplementedError("Invalid model instance.")
def _initialize_mask(self):
r"""
Returns a list of :python:`True` values with length the dimension of the lattice.
"""
if isinstance(self.model, tbm.Model):
return _tbmodels.initialize_mask(self.model)
elif isinstance(self.model, ptb.tb_model):
return _pythtb.initialize_mask(self.model, self.spinful)
elif isinstance(self.model, wberri.System_w90):
return _wberri.calc_states_uc(self.model)
else:
raise NotImplementedError("Invalid model instance.")
#################################################
# Lattice functions
#################################################
[docs]
def add_onsite_disorder(self, w : float = 0, seed : int = None):
r"""
Add on-site Anderson disorder to the model. Anderson disorder consists in a diagonal term in the Hamiltonian :math:`\mathcal H_{Anderson}=\sum_i w^{}_i c_i^{\dagger}c^{}_i` where :math:`w^{}_i` are random variables uniformly distributed in the interval :math:`\bigl[-\frac{W}{2},\frac{W}{2} \bigr]`.
Parameters
----------
w :
Disorder amplitude, defining the extrema of the interval :math:`\bigl[-\frac{W}{2},\frac{W}{2} \bigr]`. Default is ``0``.
seed :
Seed to be used in the random number generation. Default is :python:`None`, which result in the default Numpy seed.
"""
if w != 0: self.disordered = True
# Set the seed for the random number generator
if seed is not None:
np.random.seed(seed)
if self.spinful:
spinstates = 2
else:
spinstates = 1
# Create an array of random numbers uniformly distributed in (-w/2, w/2)
d = 0.5 * w * (2 * np.random.rand(self.n_orb // spinstates) - 1.0)
dis = np.repeat(d,spinstates)
# This function returns None, modify self.hamiltonian
np.fill_diagonal(self.hamiltonian,self.hamiltonian.diagonal() + dis)
[docs]
def add_vacancies(self, vacancies_list):
r"""
Add vacancies in the system by removing a site in the lattice.
Parameters
----------
vacancies_list:
A list of [cell_x, cell_y, basis] pointing to the cell and atom to be removed. Multiple values at once are allowed. For instance, consider a :math:`L_x\times L_y` lattice with a basis of :math:`N_b` atoms, and :math:`M` sites to be removed. Then, the list :python:`[[np.random.randint(L_x), np.random.randint(L_y), np.random.randint(N_b)] for _ in range(M)]` can be given as argument.
.. note::
In order to enforce half-filling, if the model is not spinful, an even number of vacancies is required.
"""
ham = self.hamiltonian
pos = self.r
cart_pos = self.cart_positions.T
sz = self.sz
nocc = self.n_occ
norbs = self.n_orb
# Convert vacancies_list to internal indexing
if np.array(vacancies_list).shape == (3,):
vacancies = self._xy_to_index(vacancies_list[0], vacancies_list[1], vacancies_list[2])
else:
vacancies = []
for hl in vacancies_list:
vacancies.append(self._xy_to_index(hl[0], hl[1], hl[2]))
# Remove the selected sites
if np.array(vacancies).shape == ():
self._mask[vacancies] = False
if not self.spinful:
raise RuntimeError("Adding an odd number of vacancies causes the model to no longer be half-filled.")
else:
self._mask[vacancies + 1] = False
for _ in range(2):
ham = np.delete(ham, vacancies, 0)
ham = np.delete(ham, vacancies, 1)
for dim in range(len(self.r)):
pos[dim] = np.delete(pos[dim], vacancies, 0)
pos[dim] = np.delete(pos[dim], vacancies, 1)
cart_pos = np.delete(cart_pos, vacancies, 1)
sz = np.delete(sz, vacancies, 0)
sz = np.delete(sz, vacancies, 1)
nocc -= 1
norbs -= 2
else:
# Avoid relabelling
vacancies = np.sort(vacancies)[::-1]
# Check if the number of vacancies is odd if the model is not spinful
if (not self.spinful) and (len(vacancies) % 2 == 1):
raise RuntimeError("Adding an odd number of vacancies causes the model to no longer be half-filled.")
for h in vacancies:
self._mask[h] = False
if not self.spinful:
ham = np.delete(ham, h, 0)
ham = np.delete(ham, h, 1)
for dim in range(len(self.r)):
pos[dim] = np.delete(pos[dim], h, 0)
pos[dim] = np.delete(pos[dim], h, 1)
cart_pos = np.delete(cart_pos, h, 1)
nocc -= 0.5
norbs -= 1
else:
self._mask[h] = False
self._mask[h + 1] = False
for _ in range(2):
ham = np.delete(ham, h, 0)
ham = np.delete(ham, h, 1)
for dim in range(len(self.r)):
pos[dim] = np.delete(pos[dim], h, 0)
pos[dim] = np.delete(pos[dim], h, 1)
cart_pos = np.delete(cart_pos, h, 1)
sz = np.delete(sz, h, 0)
sz = np.delete(sz, h, 1)
nocc -= 1
norbs -= 2
self.hamiltonian = ham
self.sz = sz
self.n_occ = int(nocc)
self.disordered = True
self.cart_positions = cart_pos.T
[docs]
def make_heterostructure(self, model2, direction : int = 0, start : int = 0, stop : int = 0):
r"""
Create a heterostructure starting from :python:`model2`. The system will be split in the direction :python:`direction` starting from the unit cell whose index is :python:`start` and stopping at the cell indexed as :python:`end`.
Parameters
----------
model2 :
The model whose Hamiltonian has to be used when constructing the heterostructure.
direction :
Direction in which the splitting is realized, allowed values are ``0`` for the :math:`\mathbf{a}_1` direction or ``1`` for the :math:`\mathbf{a}_2` direction. Default is ``0``.
start :
Starting point for the splitting in the :python:`direction` direction. Default is ``0``.
end :
End point of the splitting in the :python:`direction` direction. Default is ``0``.
.. warning::
The current instance of the model will be modified.
"""
# Ensure start <= stop
if start > stop:
tmp = start
start = stop
stop = tmp
# Check input data are ok
if not (start >= 0 and start < (self.Lx if direction == 0 else self.Ly)):
raise RuntimeError("Start point value not allowed.")
if not (stop > 0 and stop < (self.Lx if direction == 0 else self.Ly)):
raise RuntimeError("End point value not allowed.")
if direction not in [0, 1]:
raise RuntimeError("Direction not allowed: insert 0 for 'x' and 1 for 'y'.")
if not issubclass(type(model2), Model):
raise RuntimeError("The two models must be instances of Model, Supercell or FiniteModel.")
if not (self.Lx == model2.Lx and self.Ly == model2.Ly and np.allclose(self.r, model2.r)):
raise RuntimeError("You can only build heterostructures of the same model.")
# Generate the Hamiltonian of the second model
hamilt_model2 = np.copy(model2.hamiltonian)
# Check if only the on-site terms are changed (if so, avoid checking all the matrix elements)
onsite_only = False
if np.allclose( self.hamiltonian - np.diag(np.diag(self.hamiltonian)), hamilt_model2 - np.diag(np.diag(hamilt_model2)) ):
onsite_only = True
# Remove on-site terms from the model
onsite1 = np.diag(self.hamiltonian).copy()
onsite2 = np.diag(hamilt_model2).copy()
self.hamiltonian -= np.diag(onsite1)
if direction == 0:
# Splitting along the x direction
ind = np.array([[(start + i) * self.Ly * self.states_uc + j * self.states_uc for j in range(self.Ly)] for i in range(stop - start + 1)]).flatten()
else:
# Splitting along the y direction
ind = np.array([[i * self.Ly * self.states_uc + start * self.states_uc + j * self.states_uc for j in range(stop - start + 1)] for i in range(self.Lx)]).flatten()
for i in ind:
for j in range(self.states_uc):
onsite1[i + j] = onsite2[i + j]
# Add the new on-site terms
self.hamiltonian += np.diag(onsite1)
# If other matrix element are changed
if not onsite_only:
# Indices of every atom in the selected cells, not only of the initial atom of the cell
indices = np.array([[i + j for j in range(self.states_uc)] for i in ind]).flatten()
# Cycle over the rows of the hopping matrix
for k in range(self.hamiltonian.shape[0]):
# Cycle over the columns of the hopping matrix
for l in range(self.hamiltonian.shape[1]):
if k == l: continue
# Hopping amplitudes
amplitude1 = self.hamiltonian[k][l]
amplitude2 = hamilt_model2[k][l]
if k in indices:
if np.absolute(amplitude2) < 1e-10: continue
self.hamiltonian[k][l] = amplitude2
else:
if np.absolute(amplitude1) < 1e-10: continue
self.hamiltonian[k][l] = amplitude1
#################################################
# Utility functions
#################################################
def _xy_to_index(self, cellx : int, celly : int, basis : int):
r"""
Convert :python:`[cell_x, cell_y, basis]` to the internal indexing of the lattice sites.
Parameters
----------
cellx :
Index of the unit cell along the :math:`\mathbf{a}_1` direction.
celly :
Index of the unit cell along the :math:`\mathbf{a}_2` direction.
basis :
Index of the atom in the basis of the lattice. This index is defined when creating the model.
Returns
-------
index :
The index of an atom in the lattice according to the internal indexing.
.. note::
The conventional internal indexing is the one also used in PythTB and TBmodels when dealing with supercells: the first index is the bottom-left corner of the lattice and the last one is the top-right. The indices increase first along the :math:`\mathbf{a}_2` direction, then along :math:`\mathbf{a}_1`.
"""
return self.Ly * self.states_uc * cellx + self.states_uc * celly + basis * (2 if self.spinful else 1)
def _xy_to_line(self, fixed_coordinate : str, xy : int):
r"""
Returns the indices of the sites of the lattice along the fixed direction :python:`fixed_coordinate` starting from :python:`xy`.
Parameters
----------
fixed_coordinate :
Coordinate that has to be kept fixed. If one is interested in the indices along the :math:`\mathbf{a}_1` direction, this has to be set :python:`'y'` else :python:`'x'`.
xy :
Index of the unit cell along the direction :python:`fixed_coordiante` from which the indices are returned.
Returns
-------
index_list :
The list of indices of the atoms in the lattice keeping fixed :python:`fixed_direction` from the cell :python:`xy`.
"""
if fixed_coordinate == 'x':
return np.array([self.states_uc * xy * self.Ly + i for i in range(self.states_uc * self.Ly) if self._mask[self.states_uc * xy * self.Ly + i]]).flatten().tolist()
elif fixed_coordinate == 'y':
indices = []
for i in range(self.Lx):
for j in range(self.states_uc):
if self._mask[self.states_uc * xy + self.states_uc * self.Ly * i + j]:
indices.append(self.states_uc * xy + self.states_uc * self.Ly * i + j)
return np.array(indices)
else:
raise RuntimeError("Direction not allowed, only 'x' or 'y'.")
#################################################
# Macroscopic average functions
#################################################
def _lattice_contraction(self, cutoff : float):
r"""
Defines which atomic sites must be contracted on one site, for each site of the lattice.
Parameters
----------
cutoff :
Real-space cutoff for the contraction window.
Returns
-------
contraction :
List of the indices of the atoms that are within cutoff from each other, per atom.
"""
contraction = []
# Diagonal position matrix elements
rx = self.cart_positions[:, 0]; ry = self.cart_positions[:, 1]
# Return True if current is close to trial within cutoff
def within_range(current, trial):
return True if (rx[current] - rx[trial]) ** 2 + (ry[current] - ry[trial]) ** 2 - cutoff ** 2 < 1e-6 else False
# Store the lattice sites closer than cutoff for each lattice site
for current in range(len(rx)):
contraction.append([trial for trial in range(len(rx)) if within_range(current, trial)])
return contraction
def _PBC_lattice_contraction(self, cutoff):
r"""
Defines which atomic sites must be contracted on one site, for each site of the lattice within PBCs (minimum convention image).
Parameters
----------
cutoff :
Real-space cutoff for the contraction window.
Returns
-------
contraction :
List of the indices of the atoms that are within cutoff from each other, per atom.
"""
contraction = []
# Compute the primitive lattice vectors
if isinstance(self.model, tbm.Model):
lvecs = np.copy(self.model.uc)
elif isinstance(self.model, ptb.tb_model):
lvecs = np.copy(self.model._lat)
else:
raise NotImplementedError("Invalid model instance.")
lvecs[0] /= self.Lx
lvecs[1] /= self.Ly
# Compute reduced coordinates
reds = []
linv = np.linalg.inv(lvecs)
for r in range(len(self.cart_positions)):
reds.append(np.dot(np.array([self.cart_positions[r, 0], self.cart_positions[r, 1]]), linv))
reds = np.array(reds)
rx = reds[:, 0]; ry = reds[:, 1]
# Return True if current is close to trial within cutoff (trial comprise minimum image convention)
def within_range(current, trial):
trialorbitals = np.array([
[rx[trial], ry[trial]], [rx[trial], ry[trial] + self.Ly], [rx[trial], ry[trial] - self.Ly], [rx[trial] + self.Lx, ry[trial]],
[rx[trial] + self.Lx, ry[trial] + self.Ly], [rx[trial] + self.Lx, ry[trial] - self.Ly], [rx[trial] - self.Lx, ry[trial]],
[rx[trial] - self.Lx, ry[trial] + self.Ly], [rx[trial] - self.Lx, ry[trial] - self.Ly]
])
for to in trialorbitals:
# Go back to real space coordinates to compute the distance
dr = np.dot(np.array([rx[current] - to[0], ry[current] - to[1]]), lvecs)
if np.dot(dr, dr) - cutoff ** 2 < 1e-6:
return True
return False
# Store the lattice sites closer than cutoff for each lattice site
for current in range(len(rx)):
contraction.append([trial for trial in range(len(rx)) if within_range(current, trial)])
return contraction
def _average_over_radius(self, vals, cutoff : float, contraction = None):
r"""
Evaluate the average of :python:`vals` according to a certain :python:`cutoff` in real space.
Parameters
----------
vals :
Values that have to be averaged.
cutoff :
Real space cutoff for the real space average.
contraction :
List of atoms that have to be considered in the macroscopic average, per atom. Default is :python:`None`, which means it is computed from the positions of the atoms in the lattice.
Returns
-------
return_vals :
List of averaged values per lattice site within a cutoff.
"""
return_vals = []
# If no contraction list is passed, it is computed from the lattice position (default OBC contraction)
if contraction is None:
contraction = self._lattice_contraction(cutoff)
# Macroscopic average within a certain radius = cutoff
for current in range(len(self.cart_positions[:, 0])):
tmp = [vals[ind] for ind in contraction[current]]
if not len(tmp) == 0:
return_vals.append(np.sum(tmp) / (len(tmp) / self.states_uc))
else:
raise RuntimeError("Unexpected error occourred in counting the neighbors of a lattice site, there are none.")
return np.array(return_vals)
