"""Hamiltonian classes implementation for large systems.
Copyright © 2019-2023 HQS Quantum Simulations GmbH. All Rights Reserved.
"""
from typing import Optional
from lattice_builder import builderTools as BT
from .abstractHamiltonian_LB import HamiltonianLB
from .qolossal_tools import blocks2full
import numpy as np
from scipy.sparse import csr_array
from scipy.sparse import eye as speye
[docs]
class HamiltonianSparse(HamiltonianLB):
r"""Hamiltonian implementation in full sparse matrix representation.
This implementation is based on the representation of the Hamiltonian in sparse matrix form.
The parent class takes care of the interface with the Lattice Builder, saving the lattice
builder instance in the object, for easy access and extracting all the important properties of
the system to be used within the package's methods.
"""
[docs]
def __init__(self, input_dict: dict, disDict: Optional[dict] = None) -> None:
r"""Initialization of the HamiltonianSparse class.
Args:
input_dict (dict): Input dictionary for the Lattice Builder to construct the
Hamiltonian
disDict (Optional[dict]): dictionary specifying the disorder.
"""
# Call the init of the parent class. This runs the whole interface with the LB. All that
# is left is to specify the Hamiltonian and/or all quantities needed to implement the dot
# product.
super().__init__(input_dict, disDict)
self.set_Hamiltonian()
[docs]
def set_Hamiltonian(
self, kpt: Optional[np.ndarray] = None, mu: float = 0, Bz: float = 0
) -> None:
"""Hamiltonian setter method.
Args:
kpt (Optional[np.ndarray]): Bloch vector determining phase-shift at periodic
boundaries.
mu (float): additional chemical potential.
Bz (float): additional Bz field.
"""
if kpt is not None:
# if kpt is set, the hamiltonian will be complex (except at some special k points)
self.is_complex = True
else:
# NOTE: is this enough?
kpt = np.zeros((3,))
self.is_complex = self.lb.is_complex
# NOTE: this type of code structure is the same as for the velocity in the parent class,
# maybe this process could/should be made a function
self.is_spinflip = False
# extract hamiltonian (and possibly adjust velocity if needed)
tmp_H = self.lb.getH0(k=kpt)
if tmp_H is None:
raise RuntimeError("Hamiltonian not available. The system was not properly set up.")
tmp_D = self.lb.getD(k=kpt)
if tmp_D is None:
raise RuntimeError(
"BCS part of the Hamiltonian not available. The system was not properly set up."
)
if self.is_spin_matrix:
# if the Hamiltonian is spinful, the Hamiltonian will be made of the spin resolved
# blocks as
# ((H_uu, Hud)
# (H_du, H_dd))
# The helper function builds the full Hamiltonian starting from the blocks, which are
# the output of the LB.getH0 method. These are shifted up or down depending on the
# additional Bz contribution.
# spinflip flag!
self.is_spinflip = False if len(tmp_H[2].data) == 0 else True
# subtract/add B field induced splitting
splitting = Bz / 2 * speye(m=tmp_H[0].shape[0], dtype=self.dtype, format="csr")
tmp_H = blocks2full((tmp_H[0] - splitting, tmp_H[1] + splitting, tmp_H[2]))
tmp_D = blocks2full(tmp_D)
self.H = tmp_H
self.H -= mu * speye(m=self.H.shape[0], dtype=self.dtype, format="csr")
# if the system has BCS coupling terms, build the Hamiltonian in Nambu space.
if self.is_BCS:
# NOTE: the factor of two in front of D is needed because the pairing matrix returned
# by the LB is non-antisymmetrized, and I use this matrix in my code instead of the
# (correct) antisymmetrized one. To make up for this it is enough to add this factor 2.
self.H = blocks2full((self.H, -self.H, 2 * tmp_D))
# Only keep the real part if the Hamiltonian is real.
# NOTE: within the building of the matrix, we will always need enough RAM to build the
# matrix as complex. Only after building it we can cast it to real if possible. This
# limitation comes from the Lattice_Builder since the function getH0 always returns a
# complex.
if not self.is_complex:
self.H = self.H.real
[docs]
def dot(self, v: np.ndarray, out: np.ndarray) -> np.ndarray:
r"""Calculate the scalar product of H.
Args:
v (np.ndarray): vector to dot the Hamiltonian into.
out (np.ndarray): preallocated output array.
Returns:
np.ndarray: Dot product of H with v
"""
out[:] = self.H @ v + self.disLocPot.reshape(v.shape) * v
return out
[docs]
class HamiltonianDense(HamiltonianLB):
r"""Hamiltonian implementation in full dense matrix representation.
This implementation is based on the representation of the Hamiltonian in dense matrix form. It
is the exact same as the HamiltonianSparse implementation, but with dense matrices and
corresponding arithmetics methods.
The parent class takes care of the interface with the Lattice Builder, saving the lattice
builder instance in the object, for easy access and extracting all the important properties of
the system to be used within the package's methods.
"""
[docs]
def __init__(self, input_dict: dict, disDict: Optional[dict] = None) -> None:
r"""Initialization of the HamiltonianDense class.
Args:
input_dict (dict): Input dictionary for the Lattice Builder to construct the
Hamiltonian
disDict (Optional[dict]): dictionary specifying the disorder.
"""
# Call the init of the parent class. This runs the whole interface with the LB. All that
# is left is to specify the Hamiltonian and/or all quantities needed to implement the dot
# product.
super().__init__(input_dict, disDict)
self.set_Hamiltonian()
[docs]
def set_Hamiltonian(
self, kpt: Optional[np.ndarray] = None, mu: float = 0, Bz: float = 0
) -> None:
"""Hamiltonian setter method.
Args:
kpt (Optional[np.ndarray]): Bloch vector determining phase-shift at periodic
boundaries.
mu (float): additional chemical potential.
Bz (float): additional Bz field.
"""
if kpt is not None:
# if kpt is set, the hamiltonian will be complex (except at some special k points)
self.is_complex = True
else:
# NOTE: is this enough?
kpt = np.zeros((3,))
self.is_complex = self.lb.is_complex
self.is_spinflip = False
# extract hamiltonian
tmp_H = self.lb.getH0(k=kpt)
if tmp_H is None:
raise RuntimeError("Hamiltonian not available. The system was not properly set up.")
tmp_D = self.lb.getD(k=kpt)
if tmp_D is None:
raise RuntimeError(
"BCS part of the Hamiltonian not available. The system was not properly set up."
)
if self.lb.is_spin_matrix:
# create full Hamiltonian in the spinful space from the spin-resolved blocks
self.is_spinflip = False if len(tmp_H[2].data) == 0 else True
# subtract/add B field induced splitting
splitting = Bz / 2 * speye(m=tmp_H[0].shape[0], dtype=self.dtype, format="csr")
tmp_H = blocks2full((tmp_H[0] - splitting, tmp_H[1] + splitting, tmp_H[2]))
tmp_D = blocks2full(tmp_D)
self.H = tmp_H.toarray()
tmp_D = tmp_D.toarray()
self.H -= mu * np.eye(self.H.shape[0], dtype=self.dtype)
# create Hamiltonian in Nambu space if the system has BCS pairing.
if self.is_BCS:
# see HamiltonianSparse for explaination on factor of 2 in front of D.
self.H = blocks2full((self.H, -self.H, 2 * tmp_D))
# take only real part if H is real
if not self.is_complex:
self.H = self.H.real
[docs]
def dot(self, v: np.ndarray, out: np.ndarray) -> np.ndarray:
r"""Calculate the scalar product of H.
Args:
v (np.ndarray): vector to dot the Hamiltonian into.
out (np.ndarray): preallocated output array.
Returns:
np.ndarray: Dot product of H with v
"""
out[:] = self.H @ v + self.disLocPot.reshape(v.shape) * v
return out
[docs]
class HamiltonianTensor(HamiltonianLB):
r"""Hamiltonian in tensor representation.
This implementation expresses the Hamiltonian exactly like in the Lattice Builder, namely by
storing, for each bond within the lattice definition, the intra(inter) unit-cell bond matrix
and the expander matrix needed to map this bond matrix onto the whole space.
This means that the contribution to the full Hamiltonian of each bond is:
.. math::
{\rm E}_x \otimes {\rm E}_y \otimes {\rm E}_z \otimes B
where :math:`{\rm E_i}` is the expander in the i-th direction and :math:`B` is the bond matrix.
Because of this, calculating the contribution to the dot product :math:`H\cdot v` from each
bond is
.. math::
({\rm E}_x \otimes {\rm E}_y \otimes {\rm E}_z \otimes B) \cdot v
which is simple to calculate once we express (reshape) v as a multidimensional array with the
correct matching dimension to x, y, z and unit-cell.
"""
[docs]
def __init__(self, input_dict: dict, disDict: Optional[dict] = None) -> None:
r"""Initialization of the HamiltonianTensor class.
Args:
input_dict (dict): Input dictionary for the Lattice Builder to construct the
Hamiltonian
disDict (Optional[dict]): dictionary specifying the disorder.
Raises:
NotImplementedError: if the system is spinful or has BCS terms.
"""
# Call the init of the parent class. This runs the whole interface with the LB. All that
# is left is to specify the Hamiltonian and/or all quantities needed to implement the dot
# product.
super().__init__(input_dict, disDict)
if self.is_BCS:
raise NotImplementedError("Hamiltonian class not implemented yet for BCS systems.")
# in this case, the situation is a bit tricky, if the system is defined as spinful but the
# spin channels are completely equivalent, I still need to access self.lb.bond_h0_uu
# instead of self.lb.bond_h0, which does not exist.
if self.is_spin_matrix:
raise NotImplementedError("Hamiltonian class not implemented yet for spinful case.")
self.set_Hamiltonian()
[docs]
def set_Hamiltonian(
self, kpt: Optional[np.ndarray] = None, mu: float = 0, Bz: float = 0
) -> None:
"""Hamiltonian setter method.
Raises:
NotImplementedError: if Bz != 0.
Args:
kpt (Optional[np.ndarray]): Bloch vector determining phase-shift at periodic
boundaries.
mu (float): additional chemical potential.
Bz (float): additional Bz field.
"""
# NOTE: Bz not used. This Hamiltonian has not been implemented yet for spinful systems
if Bz != 0:
raise NotImplementedError("Nonzero Bz not implemented yet for Tensor Hamiltonian.")
if kpt is not None:
# if kpt is set, the hamiltonian will be complex (except at some special k points)
self.is_complex = True
else:
# NOTE: is this enough? VERY HACKY
kpt = np.zeros((3,))
self.is_complex = self.lb.is_complex
# reset expander with the correct kpoint such that the Hamiltonian is built correctly
self.lb._reset_expand(kpt)
# extract non-interacting components fo the Hamiltonian, namely the diagonal intra-uc
# contribution (atom) and inter-uc as well as off-diagonal intra-uc ones (bonds)
self.atom, bonds = BT.selectPart(self.lb.atom_h0, self.lb.bond_h0, None)
# keep real part only if Hamiltonian is real to save memory
self.atom = self.atom.toarray() if self.is_complex else self.atom.real.toarray()
self.atom -= mu * np.eye(self.atom.shape[0])
# these dictionaries will be used to save the expanders for the different bonds and flag
# whether the corresponding expander is trivial (i.e. identity matrix) or not. This way we
# can avoid carrying out the trivial dot products with identities, hence saving time.
self.Hx = {}
self.trivial_x = {}
self.Hy = {}
self.trivial_y = {}
self.Hz = {}
self.trivial_z = {}
self.Huc = {}
# for each bond, extract the corresponding expander. These expanders are what is then
# kronecker multiplied to expand each bond hamiltonian to the full matrix representation
# the is_complex flag here is used in the same way as for atom
for translation, bond in bonds.items():
tmpx, tmpy, tmpz = self.lb.expand["H0"][translation]
self.Hx[translation] = tmpx.toarray() if self.is_complex else tmpx.real.toarray()
self.trivial_x[translation] = np.allclose(self.Hx[translation], np.eye(self.Lx))
self.Hy[translation] = tmpy.toarray() if self.is_complex else tmpy.real.toarray()
self.trivial_y[translation] = np.allclose(self.Hy[translation], np.eye(self.Ly))
self.Hz[translation] = tmpz.toarray() if self.is_complex else tmpz.real.toarray()
self.trivial_z[translation] = np.allclose(self.Hz[translation], np.eye(self.Lz))
self.Huc[translation] = bond.toarray() if self.is_complex else bond.real.toarray()
[docs]
def dot(self, v: np.ndarray, out: np.ndarray) -> np.ndarray:
r"""Calculate the scalar product of H.
Args:
v (np.ndarray): vector to dot the Hamiltonian into.
out (np.ndarray): preallocated output array.
Returns:
np.ndarray: Dot product of H with v
"""
def _dot(
hx: csr_array,
tx: bool,
hy: csr_array,
ty: bool,
hz: csr_array,
tz: bool,
huc: csr_array,
vec: np.ndarray,
) -> np.ndarray:
"""Helper function for dot product.
Args:
hx (csr_array): x expander.
tx (bool): x trivial flag
hy (csr_array): y expander.
ty (bool): y trivial flag
hz (csr_array): z expander.
tz (bool): z trivial flag
huc (csr_array): inter-uc hamiltonian.
vec (np.ndarray): vector to dot into the H.
Returns:
np.ndarray: H dot vec.
"""
# simply dot each expander/hamiltonian into the vector. This vector is assumed to be in
# the tensor-like shape (Mx, My, Mz, Muc) where the M's are the different lengths in
# each dimension
tmp = np.einsum("uv, abcv -> abcu", huc, vec)
# the flags tell us if it's trivial or not. If so, no need to perform the product.
if not tz:
tmp = np.einsum("zc, abcu -> abzu", hz, tmp)
if not ty:
tmp = np.einsum("yb, abzu -> ayzu", hy, tmp)
if not tx:
tmp = np.einsum("xa, ayzu -> xyzu", hx, tmp)
return tmp
# initialize result vector (same "linear" shape as the input vector)
out[:] = 0.0
# reshape the input vector to the tensor-like shape needed for the dot
vv = v.reshape((self.Lx, self.Ly, self.Lz, self.lb.M_uc))
# loop over all the bonds
for key in self.Hx.keys():
# if the bond is (0, 0, 0), add the diagonal part (atom) to the intra-uc hamiltonian.
# Do so only in this first contribution and not in the second (herm conj) part, since
# being the diagonal part, it only has to be added in once
tmpUC = self.Huc[key] + self.atom if key == (0, 0, 0) else self.Huc[key]
# calculate the contribution of the dot product for this bond
out += _dot(
self.Hx[key],
self.trivial_x[key],
self.Hy[key],
self.trivial_y[key],
self.Hz[key],
self.trivial_z[key],
tmpUC,
vv,
).reshape(out.shape)
# calculate the contribution of the hc part of the bond-hamiltonian to the dot product
# NOTE: this has to be done separately because the bond-hamiltonian is directional,
# meaning that if, e.g. one particular bond expresses the connection between one uc
# (index i as reference) and the uc to its right (index i+1), the hc contribution will
# NOT be the one with the inter-uc hamiltoninan that has been hc, but the one with the
# transposed expanders (which will now connect the uc i+1 to the uc i).
out += _dot(
np.conj(self.Hx[key]).T,
self.trivial_x[key],
np.conj(self.Hy[key]).T,
self.trivial_y[key],
np.conj(self.Hz[key]).T,
self.trivial_z[key],
np.conj(self.Huc[key]).T,
vv,
).reshape(out.shape)
out += self.disLocPot.reshape(out.shape) * v
# return scaled result
return out
[docs]
class HamiltonianOperator(HamiltonianLB):
r"""Hamiltonian in direct operator representation.
This implementation describes the Hamiltonian by simply storing the recipes for all the bonds
defining the system. For example, in a 1D chain with a unit cell made of two atoms with an
intra-uc hopping of t_1 and an inter-uc hopping of t_2, the information stored in this
implementation would be:
.. math::
(0, 0, 0),
\begin{pmatrix}
0 & t_1 \\
t_1^* & 0
\end{pmatrix},
\quad (1, 0, 0),
\begin{pmatrix}
0 & t_2 \\
t_2^* & 0
\end{pmatrix}
Where the first tuple indicates the translation associated with the bond and the matrix gives
us information about the hopping amplitudes within unitcells.
Since the Hamiltonian matrix is never build, the memory impact of this implementation is
minimal.
"""
[docs]
def __init__(self, input_dict: dict, disDict: Optional[dict] = None) -> None:
r"""Initialization of the HamiltonianOperator class.
NOTE: This Hamiltonian is only partially tested at the moment, since it does not work for
periodic systems and making it fit into the general structure of tests would require more
work than necessary. As soon as the class is generalized, this should be added to the list
of tested Hamiltonians.
Args:
input_dict (dict): Input dictionary for the Lattice Builder to construct the
Hamiltonian
disDict (Optional[dict]): dictionary specifying the disorder.
Raises:
NotImplementedError: in case the input is periodic, spinful or BCS, these cases are
still not implemented.
"""
# Call the init of the parent class. This runs the whole interface with the LB. All that
# is left is to specify the Hamiltonian and/or all quantities needed to implement the dot
# product.
super().__init__(input_dict, disDict)
# NOTE: temporary, until a general implementation of the dot product is provided.
if any(x != "hard-wall" for x in self.lb.config["system"]["system_boundary_conditions"]):
raise NotImplementedError("HamiltonianOperator only implemented for hard-walls bc.")
if self.is_BCS:
raise NotImplementedError("Hamiltonian class not implemented yet for BCS systems.")
if self.is_spin_matrix:
raise NotImplementedError("Hamiltonian class not implemented yet for spinful case.")
self.set_Hamiltonian()
[docs]
def set_Hamiltonian(
self, kpt: Optional[np.ndarray] = None, mu: float = 0, Bz: float = 0
) -> None:
"""Hamiltonian setter method.
Raises:
NotImplementedError: if Bz != 0.
Args:
kpt (Optional[np.ndarray]): Bloch vector determining phase-shift at periodic
boundaries.
mu (float): additional chemical potential.
Bz (float): additional Bz field.
"""
if Bz != 0:
raise NotImplementedError("Nonzero Bz not implemented yet for Tensor Hamiltonian.")
# TODO: only for flake8 reasons. This will have to be properly generalized. This
# Hamiltonian is not implemented for anything different from hard-wall bc, so it does not
# matter anyways
if kpt is not None:
# if kpt is set, the hamiltonian will be complex (except at some special k points)
self.is_complex = True
else:
# NOTE: is this enough?
kpt = np.zeros((3,))
self.is_complex = self.lb.is_complex
# the try-except structure is needed because in the fake-spinful case (spinful case with
# equivalent spin channels), the spin-resolved quantities are available, not the spinless
# ones, meaning lb.atom_h0_uu has to be grabbed (equivalent to lb.atom_h0_dd). Same applies
# to bonds
try:
tmp_atom = self.lb.atom_h0
tmp_bonds = self.lb.bond_h0
except AttributeError:
tmp_atom = self.lb.atom_h0_uu
tmp_bonds = self.lb.bond_h0_uu
# extract non-interacting components fo the Hamiltonian, namely the diagonal intra-uc
# contribution (atom) and inter-uc as well as off-diagonal intra-uc ones (bonds)
# this is all I need to store for this representation.
# the try-except is needed to
self.atom, self.bonds = BT.selectPart(tmp_atom, tmp_bonds, None)
# keep real part only if Hamiltonian is real to save memory
self.atom = self.atom.toarray() if self.is_complex else self.atom.real.toarray()
self.atom -= mu * np.eye(self.atom.shape[0])
for key in self.bonds.keys():
if self.is_complex:
self.bonds[key] = self.bonds[key].toarray()
else:
self.bonds[key] = self.bonds[key].toarray().real
[docs]
def dot(self, v: np.ndarray, out: np.ndarray) -> np.ndarray:
r"""Calculate the scalar product of H.
A simple way of understanding how this dot product is performed is by thinking of a 1D
chain with nn hopping -1 and considering H*v. Here the translation of the bond is 1 (to the
right or to the left) and, given a certain vector v, we have
.. math::
H \cdot v =
\begin{pmatrix}
0 & -1 & 0 & 0 \\
-1 & 0 & -1 & 0 \\
0 & -1 & 0 & -1 \\
0 & 0 & -1 & 0 \\
\end{pmatrix}
\cdot
\begin{pmatrix}
v_1 \\
v_2 \\
v_3 \\
v_4
\end{pmatrix}
= -1
\begin{pmatrix}
v_2 \\
v_3 \\
v_4 \\
0
\end{pmatrix}
-1
\begin{pmatrix}
0 \\
v_1 \\
v_2 \\
v_3
\end{pmatrix}
where we can clearly see that the first vector is the contribution of the hopping to the
right and the second, the one of the hopping to the left. To determine the contribution of
each bond to the dot product, one has to simply apply the corresponding translation
(separately for back and forth hopping) and accumulate the result. This has to be done for
each dimension separately, but it is as simple as reshaping the vector v to have the shape
(x, y, z, uc) and applying the translations in the corresponding direction. The uc
dimension is kept as is and dotted directly into the bond associated matrices.
Args:
v (np.ndarray): vector to dot the Hamiltonian into.
out (np.ndarray): preallocated output array.
Returns:
np.ndarray: Dot product of H with v
"""
vec = v.reshape((self.Lx, self.Ly, self.Lz, self.lb.M_uc))
# add diagonal part from 'atom' and disorder potential
res = (
vec @ self.atom
+ (self.disLocPot.reshape((self.Lx, self.Ly, self.Lz, self.lb.M_uc))) * vec
)
for transl, bondtmp in self.bonds.items():
_transl = np.array(transl)
# I have to consider both contribution coming from applying the translation forward AND
# backward. The - sign in front of the first transl below is needed because of how
# I define the translation (within the determination of the indices for selection and
# assignement). In other words, this is just a convention for defining positive (or
# negative) translation to be the one goin to the right (or to the left). The ^dag in
# the second bond matrix (bondtmp) is needed because the bonds are directional, so we
# need to "reverse" their direction when considering the opposite translation.
for tr, bond in zip([-_transl, _transl], [bondtmp, np.conj(bondtmp.T)]):
# The following indices are determined in such a way that the correct slices are
# selected for both positive and negative translation. The trick is in the usage of
# max and of the fact that selecting a vector beyond it's length is allowed as
# long as we start within the vector, e.g. in the case above v[1:100] would return
# [v2, v3, v4] for a vector v = [v0, v1, v2, v3, v4] without errors.
# The same has to be applied to all directions and, for a non-trivial unit cell,
# the simple multiplication of the hopping into the shifted vectors becomes a dot
# product into the correct dimension (here the last dimension is the one of the uc)
# selection boundaries
sel_begx = max(0, tr[0])
sel_endx = self.Lx + tr[0]
sel_begy = max(0, tr[1])
sel_endy = self.Ly + tr[1]
sel_begz = max(0, tr[2])
sel_endz = self.Lz + tr[2]
# assignement boundaries
asm_begx = max(0, -tr[0])
asm_endx = self.Lx - tr[0]
asm_begy = max(0, -tr[1])
asm_endy = self.Ly - tr[1]
asm_begz = max(0, -tr[2])
asm_endz = self.Lz - tr[2]
# add contribution by selecting the correct slices and adding them to the correct
# slice of the result vector (see example in class docstring for more info)
res[asm_begx:asm_endx, asm_begy:asm_endy, asm_begz:asm_endz] += (
vec[sel_begx:sel_endx, sel_begy:sel_endy, sel_begz:sel_endz] @ bond
)
out[:] = res.reshape(out.shape)
# return the final result in the correct shape and normalization.
return out