import numpy as np
import tbmodels as tbm
import pythtb as ptb
import wannierberri as wberri
from .classes import Model
from . import _tbmodels
from . import _pythtb
from . import _wberri
import scipy.linalg as la
from opt_einsum import contract
from .utils import *
[docs]
class Supercell(Model):
r"""
A class creating a supercell given a model built from either TBmodels or PythTB instances.
It also can receive as input a model which is initialized as WannierBerri instance :python:`wannierberri.System_w90` from the Wannier functions generated by Wannier90 code starting from a :math:`\Gamma`-only *ab initio* calculation.
It contains methods to add disorder and vacancies to the supercell and to calculate single-point topological invariants and local topological markers.
Parameters
----------
tbmodel :
Tight-binding model constructed from TBmodels or PythTB, or *ab initio* tight-binding model inizialized as :python:`wannierberri.System_w90` from WannierBerri.
Lx :
Number of unit cells repeated along the :math:`\mathbf{a}_1` direction in the supercell. In case of :python:`wannierberri.System_w90` instance, please set :python:`Lx = 1`.
Ly :
Number of unit cells repeated along the :math:`\mathbf{a}_2` direction in the supercell. In case of :python:`wannierberri.System_w90` instance, please set :python:`Ly = 1`.
spinful :
Whether the model should be interpreted as spinful or not. Default is :python:`False`.
file_spn :
Whether the file seedname.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.
"""
def __init__(self, tbmodel, Lx : int, Ly : int, spinful : bool = False,
file_spn : bool = False, n_occ : int = None):
# Store local variables
self.Lx = Lx
self.Ly = Ly
self.Lz = 1
# Create supercell
self.model = self._make_supercell(tbmodel)
super().__init__(tbmodel = self.model,
spinful = spinful,
states_uc = super()._calc_states_uc(tbmodel, spinful),
Lx = self.Lx,
Ly = self.Ly,
file_spn = file_spn,
n_occ = n_occ
)
def _make_supercell(self, model):
r"""
Create :math:`Lx \times Ly` supercell for a 2D tight-binding model. If a 3D tight-binding model is provided, the supercell is created in the *xy*-plane.
For *ab initio* tight-binding model, the model is expected to be provided already in a supercell (:math:`Gamma`-only *ab initio* calculation).
Parameters
----------
model :
Model instance.
"""
if (self.Lx != 1 or self.Ly != 1):
if isinstance(model, tbm.Model):
return model.supercell([self.Lx,self.Ly]) if model.dim == 2 else model.supercell([self.Lx,self.Ly,self.Lz])
elif isinstance(model, ptb.tb_model):
return model.make_supercell([[self.Lx,0],[0,self.Ly]]) if model._dim_r == 2 else model.make_supercell([[self.Lx,0,0],[0,self.Ly,0],[0,0,self.Lz]])
else:
raise NotImplementedError("Invalid model instance.")
else:
return model
def _reciprocal_vec(self):
"""
Returns reciprocal lattice vectors in cartesian coordinates.
"""
if isinstance(self.model, tbm.Model):
return _tbmodels._reciprocal_vec(self.model)
elif isinstance(self.model, ptb.tb_model):
return _pythtb._reciprocal_vec(self.model)
elif isinstance(self.model, wberri.System_w90):
return _wberri._reciprocal_vec(self.model)
else:
raise NotImplementedError("Invalid model instance.")
def _pszp_matrix (self, u_n0 = None):
r"""
Returns the matrix of spin operator projected on the Hamiltonian occupied states at :math:`\Gamma`-point. The matrix elements are :math:`\langle u_{n \Gamma} | \hat{S}_{z} | u_{m \Gamma} \rangle`, for :math:`| u_{n \Gamma} \rangle` occupied eigenstates.
If the model is obtained from the Wannier functions generated by Wannier90 code starting from a :math:`\Gamma`-only *ab initio* calculation and ``seedname.spn`` file is provided, the spin matrix in the *z*-direction in the Wannier basis is constructed and the projection onto occupied eigenstates is considered.
Otherwise, the assumption of a tight-binding basis which is diagonal in the spin operator in the *z*-direction is considered.
Parameters
----------
u_n0 :
Matrix of Hamiltonian eigenstates at :math:`\Gamma`-point.
"""
pszp = np.ndarray([self.n_occ,self.n_occ],dtype=complex)
if self.wannier and self.file_spn:
pszp = self.sz[:self.n_occ,:self.n_occ]
else:
sz = self.sz
pszp = u_n0[:self.n_occ,:].conjugate() @ (sz @ u_n0[:self.n_occ,:].T)
return pszp
def _periodic_gauge(self, u_n0, b, n_occ = None):
r"""
Returns the matrix of occupied eigenvectors at the edge of the Brillouin zone imposing periodic gauge upon the eigenvectors at :math:`\Gamma`.
Parameters
----------
u_n0 :
Matrix of Hamiltonian eigenstates at :math:`\Gamma`-point.
:math:`\mathbf{b}` :
Reciprocal lattice vector indicating the edge of the Brillouin zone.
n_occ :
Number of occupied states.
"""
if n_occ is None: n_occ = self.n_occ
orb_c = self.cart_positions
vec_scal_b = orb_c @ b
vec_exp_b = np.exp(-1.j*vec_scal_b)
u_nb = vec_exp_b.T * u_n0[:n_occ,:]
return u_nb
def _dual_state(self, u_n0, u_nb, spin = None, n_sub = None):
r"""
Returns the "dual states" reported in Eq. (8) in Ref. `Ceresoli-Resta (2007) <https://journals.aps.org/prb/abstract/10.1103/PhysRevB.76.012405>`_ for the calculation of the single-point Chern number and those in Eq. (10) in Ref. `Favata-Marrazzo (2023) <https://iopscience.iop.org/article/10.1088/2516-1075/acba6f/meta>`_ for the calculation of the single-point spin Chern number.
Parameters
----------
u_n0 :
Matrix of Hamiltonian eigenstates at :math:`\Gamma`-point.
u_nb :
Matrix of Hamiltonian eigenstates at :math:`\mathbf{b}`-point at the edge of the Brillouin zone.
spin :
In case of :python:`spinful` system, indicates which sector of spin projected operator spectra is considered in the calculation of the single-point spin Chern number.
If :python:`spin` is :python:`'up'` (:python:`'down'`), eigenstates of spin projected operator with positive (negative) eigevalues are considered.
n_sub :
Number of states used in the single-point invariant calculation.
If :python:`n_sub == None`, the number of occupied eigenstates is used for the calculation of the single-point Chern number and the number of eigenstates of spin projected operator with eigenvalues of a given sign (:math:`\pm`) is used for the calculation of the single-point spin Chern number.
"""
if n_sub is None:
if self.spinful:
n_sub = self.n_occ//2
else:
n_sub = self.n_occ
s_matrix_b = np.zeros([n_sub,n_sub], dtype=np.complex128)
udual_nb = np.zeros((n_sub, self.n_orb), dtype=np.complex128)
if (spin == 'down' or spin == None):
s_matrix_b = np.conjugate(u_n0[:n_sub,:]) @ (u_nb[:n_sub,:]).T
s_inv_b = np.linalg.pinv(s_matrix_b)
udual_nb = (s_inv_b.T) @ u_nb[:n_sub,:]
elif spin == 'up':
s_matrix_b = np.conjugate(u_n0[n_sub:,:]) @ (u_nb[n_sub:,:]).T
s_inv_b = np.linalg.pinv(s_matrix_b)
udual_nb = (s_inv_b.T) @ u_nb[n_sub:,:]
return udual_nb
#################################################
# Single-point invariants
#################################################
[docs]
def single_point_chern(self, formula : str = 'both', return_ham_gap : bool = False):
r"""
Evaluate the single-point Chern number provided in Ref. `Ceresoli-Resta (2007) <https://journals.aps.org/prb/abstract/10.1103/PhysRevB.76.012405>`_.
Parameters
----------
formula :
Two possible formulas are used for the calculation of the single-point Chern number: the :python:`'asymmetric'` one is the one originally derived in Ref. `Ceresoli-Resta (2007)` and uses right-hand discrete derivative; the :python:`'symmetric'` one uses symmetric derivative centered in :math:`\Gamma`-point and converges faster to the exact value with the supercell size. Default value is :python:`'both'`, which returns the invariants computed with both symmetric and asymmetric formulas.
return_ham_gap :
If :python:`True`, returns the value of the gap between the highest occupied and the lowest unoccupied eigenstates of the Hamiltonian at :math:`\Gamma`-point. Default value is :python:`False`.
Returns
-------
chern :
Dictionary with results from :python:`'asymmetric'` and/or :python:`'symmetric'` formula for the single-point Chern number.
hamiltonian_gap :
Gap between the lowest unoccupied and the highest occupied eigenstates of the Hamiltonian at :math:`\Gamma`-point if :python:`return_ham_gap == True`.
"""
eig, u_n0 = la.eigh(self.hamiltonian)
u_n0 = u_n0.T
hamiltonian_gap = eig[self.n_occ] - eig[self.n_occ-1]
b1, b2 = self._reciprocal_vec()
u_nb1 = self._periodic_gauge(u_n0, b1)
u_nb2 = self._periodic_gauge(u_n0, b2)
udual_nb1 = self._dual_state(u_n0, u_nb1)
udual_nb2 = self._dual_state(u_n0, u_nb2)
chern = {}
if (formula=='asymmetric' or formula =='both'):
sum_occ = 0.
for i in range(self.n_occ):
sum_occ += np.vdot(udual_nb1[i],udual_nb2[i])
chern['asymmetric'] = -np.imag(sum_occ)/np.pi
if (formula=='symmetric' or formula =='both'):
u_nmb1 = self._periodic_gauge(u_n0, -b1)
u_nmb2 = self._periodic_gauge(u_n0, -b2)
udual_nmb1 = self._dual_state(u_n0, u_nmb1)
udual_nmb2 = self._dual_state(u_n0, u_nmb2)
sum_occ = 0.
for i in range(self.n_occ):
sum_occ += np.vdot((udual_nmb1[i]-udual_nb1[i]),(udual_nmb2[i]-udual_nb2[i]))
chern['symmetric'] = -np.imag(sum_occ)/(4*np.pi)
if return_ham_gap:
return chern, hamiltonian_gap
else:
return chern
[docs]
def single_point_spin_chern(self, spin : str = 'down', formula : str = 'both', return_pszp_gap : bool = False, return_ham_gap : bool = False):
r"""
Evaluate the single-point spin Chern number derived in Ref. `Favata-Marrazzo (2023) <https://iopscience.iop.org/article/10.1088/2516-1075/acba6f/meta>`_.
Parameters
----------
spin :
Indicates which sector of the spin projected operator spectra is considered in the calculation of the single-point spin Chern number. If :python:`spin` is :python:`'up'` (:python:`'down'`), eigenstates of spin projected operator with positive (negative) eigevalues are considered. Default value is :python:`'down'`.
formula :
Two possible formulas are used for the calculation of the single-point spin Chern number: the :python:`'asymmetric'` one uses right-hand discrete derivative; the :python:`'symmetric'` one uses symmetric derivative centered in :math:`\Gamma`-point and converges faster to the exact value with the supercell size. Default value is :python:`'both'`, which returns the invariants computed with both symmetric and asymmetric formulas.
return_pszp_gap :
If :python:`True`, returns the value of the gap between positive and negative eigenvalues of the spin operator projected onto occupied states at :math:`\Gamma`-point. Default value is :python:`False`.
return_ham_gap :
If :python:`True`, returns the value of the gap between the lowest unoccupied and the highest occupied eigenstates of the Hamiltonian at :math:`\Gamma`-point. Default value is :python:`False`.
Returns
-------
spin_chern :
Dictionary with results from :python:`'asymmetric'` and/or :python:`'symmetric'` formula for the single-point spin Chern number.
pszp_gap :
Gap between positive and negative eigenvalues of the spin operator projected onto occupied states at :math:`\Gamma`-point if :python:`return_pszp_gap == True`.
hamiltonian_gap :
Gap between the lowest unoccupied and the highest occupied eigenstates of the Hamiltonian at :math:`\Gamma`-point if :python:`return_ham_gap == True`.
"""
n_sub = self.n_occ//2
if self.wannier and self.file_spn:
eig = self.eig
u_n0 = self.un_0
else:
eig, u_n0 = la.eigh(self.hamiltonian)
u_n0 = u_n0.T
hamiltonian_gap = eig[self.n_occ] - eig[self.n_occ-1]
pszp = self._pszp_matrix(u_n0)
eval_pszp, eig_pszp = la.eigh(pszp)
eig_pszp = eig_pszp.T
pszp_gap = eval_pszp[n_sub] - eval_pszp[n_sub-1]
spin_chern = {}
if (pszp_gap < 10**(-14)):
raise Exception('Closing P Sz P gap!!')
elif (eval_pszp[n_sub]*eval_pszp[n_sub-1]>0) :
#check symmetry of P Sz P spectrum
raise Exception('P Sz P spectrum is NOT symmetric!!!')
else :
q_0 = np.zeros([self.n_occ, self.n_orb], dtype=complex)
q_0 = eig_pszp @ u_n0[:self.n_occ,:]
b1, b2 = self._reciprocal_vec()
q_b1 = self._periodic_gauge(q_0, b1)
q_b2 = self._periodic_gauge(q_0, b2)
qdual_b1 = self._dual_state(q_0, q_b1, spin)
qdual_b2 = self._dual_state(q_0, q_b2, spin)
if (formula=='asymmetric' or formula=='both'):
sum_occ = 0.
for i in range(n_sub):
sum_occ += np.vdot(qdual_b1[i],qdual_b2[i])
spin_chern['asymmetric'] = -np.imag(sum_occ)/np.pi
if (formula=='symmetric' or formula=='both'):
q_mb1 = self._periodic_gauge(q_0,-b1)
q_mb2 = self._periodic_gauge(q_0,-b2)
qdual_mb1 = self._dual_state(q_0, q_mb1, spin)
qdual_mb2 = self._dual_state(q_0, q_mb2, spin)
sum_occ = 0.
for i in range(n_sub):
sum_occ += np.vdot((qdual_mb1[i]-qdual_b1[i]),(qdual_mb2[i]-qdual_b2[i]))
spin_chern['symmetric'] = -np.imag(sum_occ)/(4*np.pi)
if return_pszp_gap and return_ham_gap:
return spin_chern, pszp_gap, hamiltonian_gap
elif return_pszp_gap:
return spin_chern, pszp_gap
elif return_ham_gap:
return spin_chern, hamiltonian_gap
else:
return spin_chern
#################################################
# Local PBC marker
#################################################
[docs]
def pbc_local_chern_marker(self, direction : int = None, start : int = 0, return_projector : bool = False, input_projector = None, formula : str = "symmetric", macroscopic_average : bool = False, cutoff : float = 0.8, smearing_temperature : float = 0, fermidirac_cutoff : float = 0.1):
r"""
Evaluate the local Chern marker provided in Ref. `Baù-Marrazzo (2023) <https://arxiv.org/abs/2310.15783>`_ on the whole supercell if :python:`direction == None`. If ``direction`` is not :python:`None` evaluates the PBC local Chern marker along :python:`direction` starting from :python:`start`. Allowed directions are ``0`` (meaning along :math:`\mathbf{a}_1`), and ``1`` (meaning along :math:`\mathbf{a}_2`).
Parameters
----------
direction :
Direction along which to compute the PBC local Chern marker. Default is :python:`None` (returns the marker on the whole supercell). Allowed directions are ``0`` (meaning along :math:`\mathbf{a}_1`), and ``1`` (meaning along :math:`\mathbf{a}_2`).
start :
If :python:`direction` is not :python:`None`, is the coordinate of the unit cell from which the evaluation of the PBC local Chern marker starts. For instance, if interested on the value of the PBC local marker along the :math:`\mathbf{a}_1` direction at half height, it should be set :python:`direction = 0` and :python:`start = Ly // 2`.
return_projector :
If :python:`True`, returns the ground state projector at the end of the calculation. Default is :python:`False`.
input_projector :
Input the list of projectors :math:`[\mathcal{P}_{\Gamma},\mathcal{P}_{\mathbf{b}_1},\mathcal{P}_{\mathbf{b}_2}]` (or :math:`[\mathcal{P}_{\Gamma},\mathcal{P}_{\mathbf{b}_1},\mathcal{P}_{\mathbf{b}_2},\mathcal{P}_{-\mathbf{b}_1},\mathcal{P}_{-\mathbf{b}_2}]` if :python:`formula == 'symmetric'`) to be used in the calculation. Default is :python:`None`, which means that it is computed from the model stored in the class.
formula :
Formula to be used. Default is :python:`'symmetric'`, which is computationally more demanding but converges faster. Any other input will result in the :python:`'asymmetric'` formulation.
macroscopic_average :
If :python:`True`, returns the PBC local Chern marker averaged in real space over a radius equal to :python:`cutoff`. Default is :python:`False`.
cutoff :
Cutoff set for the calculation of the macroscopic average in real space of the PBC local Chern marker.
smearing_temperature :
Set a fictitious temperature :math:`T_s` to be used when weighting the eigenstates of the Hamiltonian comprising the ground state projector. In particular, the ground state projector is computed as :math:`\mathcal P=\sum_{n}f(\epsilon_n, T_s, \mu)|u_n\rangle\langle u_n|` where :math:`f(\epsilon_n, T_s, \mu)` is the Fermi-Dirac distribution, :math:`\mu` is the chemical potential and :math:`\mathcal{H}_{\mathbf{k}}|u_n\rangle=\epsilon_n|u_n\rangle`. Introducing some smearing is particularly useful when dealing with heterojunctions o inhomogeneous models whose insulating gap is small in order to improve the convergence of the local marker. Default is ``0``, so no smearing is introduced and a model half-filled is implied.
fermidirac_cutoff :
Cutoff imposed on the Fermi-Dirac distribution to further improve the convergence, mostly when :math:`T_s\neq0`. Default is ``0.1``, which is appropriate in most cases.
Returns
-------
lattice_chern :
PBC local Chern marker evaluated on the whole lattice if :python:`direction == None`.
lcm_direction :
PBC local Chern marker evaluated along :python:`direction` starting from :python:`start`.
return_proj :
List of projectors :math:`[\mathcal{P}_{\Gamma},\mathcal{P}_{\mathbf{b}_1},\mathcal{P}_{\mathbf{b}_2}]` (or :math:`[\mathcal{P}_{\Gamma},\mathcal{P}_{\mathbf{b}_1},\mathcal{P}_{\mathbf{b}_2},\mathcal{P}_{-\mathbf{b}_1},\mathcal{P}_{-\mathbf{b}_2}]` if :python:`formula == 'symmetric'`) used in the calculation if :python:`return_projector == True`.
"""
return_proj = []
if input_projector is None:
eigenvals, eigenvecs = la.eigh(self.hamiltonian)
# Find the chemical potential
mu = chemical_potential(eigenvals, smearing_temperature, self.n_occ)
# Evaluate the effective number of occupied states whose occupation is greater than the Fermi-Dirac cutoff
rank = np.sum(fermidirac(eigenvals, smearing_temperature, mu) > fermidirac_cutoff)
eigenvecs_use = eigenvecs.T
# Reciprocal lattice vectors
b1, b2 = self._reciprocal_vec()
# Periodic gauge along b_1 and b_2
u_nb1 = self._periodic_gauge(eigenvecs_use, b1, n_occ = rank)
u_nb2 = self._periodic_gauge(eigenvecs_use, b2, n_occ = rank)
# Dual states at b_1 and b_2
udual_b1 = self._dual_state(eigenvecs_use, u_nb1, n_sub = rank)
udual_b2 = self._dual_state(eigenvecs_use, u_nb2, n_sub = rank)
# Evaulate the ground state projector, and projectors P_b1 and P_b2
gsp = contract("ij,ik->jk", eigenvecs_use[:rank, :], (fermidirac(eigenvals[:rank], smearing_temperature, mu) * eigenvecs_use[:rank, :].conjugate().T).T)
pb1 = contract("ij,ik->jk", udual_b1, (fermidirac(eigenvals[:rank], smearing_temperature, mu) * (udual_b1.conjugate().T)).T)
pb2 = contract("ij,ik->jk", udual_b2, (fermidirac(eigenvals[:rank], smearing_temperature, mu) * (udual_b2.conjugate().T)).T)
return_proj.append(gsp); return_proj.append(pb1); return_proj.append(pb2)
p = pb1 @ pb2 - pb2 @ pb1
# If I want the symmetric formula I need to do the same also for -b_1 and -b_2
if formula == "symmetric":
u_nmb1 = self._periodic_gauge(eigenvecs_use, -1 * b1, n_occ = rank)
u_nmb2 = self._periodic_gauge(eigenvecs_use, -1 * b2, n_occ = rank)
udual_mb1 = self._dual_state(eigenvecs_use, u_nmb1, n_sub = rank)
udual_mb2 = self._dual_state(eigenvecs_use, u_nmb2, n_sub = rank)
pmb1 = contract("ij,ik->jk", udual_mb1, (fermidirac(eigenvals[:rank], smearing_temperature, mu) * (udual_mb1.conjugate().T)).T)
pmb2 = contract("ij,ik->jk", udual_mb2, (fermidirac(eigenvals[:rank], smearing_temperature, mu) * (udual_mb2.conjugate().T)).T)
return_proj.append(pmb1)
return_proj.append(pmb2)
p += ( (pmb1 @ pmb2 - pmb2 @ pmb1) - (pb1 @ pmb2 - pmb2 @ pb1) - (pmb1 @ pb2 - pb2 @ pmb1) )
else:
gsp = input_projector[0]
pb1 = input_projector[1]
pb2 = input_projector[2]
p = pb1 @ pb2 - pb2 @ pb1
if formula == 'symmetric':
pmb1 = input_projector[3]
pmb2 = input_projector[4]
p += ( (pmb1 @ pmb2 - pmb2 @ pmb1) - (pb1 @ pmb2 - pmb2 @ pb1) - (pmb1 @ pb2 - pb2 @ pmb1) )
# PBC Chern marker operator
chern_operator = -np.imag(p @ gsp) * float((self.Lx * self.Ly) / (np.pi * (8 if formula == "symmetric" else 2)))
# If macroscopic_average I have to compute the lattice values with the averages first (explicit PBC contraction passed)
if macroscopic_average or self.disordered:
contraction = self._PBC_lattice_contraction(cutoff)
pbclcm = self._average_over_radius(np.diag(chern_operator), cutoff, contraction = contraction)
if direction is not None:
# Evaluate index of the selected direction
indices = self._xy_to_line('x' if direction == 1 else 'y', start)
# If macroscopic_average consider the averaged lattice, else the Chern operators
if macroscopic_average or self.disordered:
pbclcm_line = [pbclcm[int(indices[self.states_uc * i] / self.states_uc)] for i in range(int(len(indices) / self.states_uc))]
else:
pbclcm_line = [np.sum([chern_operator[indices[self.states_uc * i + j], indices[self.states_uc * i + j]] for j in range(self.states_uc)]) for i in range(int(len(indices) / self.states_uc))]
if not return_projector:
return np.array(pbclcm_line)
else:
return np.array(pbclcm_line), np.array(return_proj)
if not macroscopic_average and not self.disordered:
pbclcm = [np.sum([chern_operator[self.states_uc * i + j, self.states_uc * i + j] for j in range(self.
states_uc)]) for i in range(int(len(chern_operator) / self.states_uc))]
# Repeat to ensure that the dimension of the marker matches the dimension of the position matrices, since if not macroscopic_average the value of the marker is defined per unit cell
pbclcm = np.repeat(pbclcm, self.states_uc)
if not return_projector:
return np.array(pbclcm)
else:
return np.array(pbclcm), np.array(return_proj)
