"""Abstract Hamiltonian class for large systems.
Copyright © 2019-2023 HQS Quantum Simulations GmbH. All Rights Reserved.
"""
from typing import Callable, Union, Optional, Type
import warnings
from abc import ABC, abstractmethod
import qolossal.qolossal_tools as qt
import qolossal.qolossal_io as qio
import lattice_solver.chebyshev as lsc
import numpy as np
from numpy.typing import ArrayLike
from scipy.linalg import eigvalsh_tridiagonal
[docs]
class Hamiltonian(ABC):
r"""Abstract base class for a general Hamiltonian.
This class does not implement any Hamiltonian. Nonetheless it provides all methods for
calculating the properties of systems. This is possible since, within the linear method based
on the Chebyshev expansion, only the dot product of the Hamiltonian on a state is needed to
compute observables.
"""
[docs]
def __init__(self) -> None:
"""Empty constructor.
The attributes defined below are of general character, i.e. will have definite values in
all actual implementations of the Hamiltonian, and are needed for the various calculations
performed by the methods defined in this class.
disDict: dictionary defining disorder. Will be translated to a HQSchema.
disLocPot: local disorder potential.
scaling: chebyshev scaling object.
chebyshev: chebyshev expansion object.
dim: number of sites in the system. Can double in spinful OR BCS cases.
realDim: real number of sites.
velocity_operator: velocity operator in the x-, y- and z- directions.
total_volume: total volume of the system.
is_complex: flag for complex Hamiltonians
"""
self.disDict: dict
self.disLocPot: np.ndarray
self.scaling: lsc.chebyshev_scaling
self.chebyshev: lsc.chebyshev
self.dim: int
self.realDim: int
self.velocity_operator: Optional[tuple[Callable, Callable, Callable]]
self.total_volume: float
self.is_complex: bool
@property
def dtype(self) -> Type:
"""Dtype convenience function.
Returns:
Type: complex or float depending on the flag is_complex
"""
return complex if self.is_complex else float
[docs]
@abstractmethod
def dot(self, v: np.ndarray, out: np.ndarray) -> np.ndarray:
"""Empty dot product. Abstract method: it has to be overloaded.
Args:
v (np.ndarray): Vector to dot.
out (np.ndarray): preallocated output array.
Raises:
NotImplementedError: In any case, since this method HAS to be overwritten by children
classes.
"""
raise NotImplementedError("The dot product function has not been defined.")
[docs]
def _dot_scaling(self, v: np.ndarray, out: np.ndarray) -> np.ndarray:
"""Dot product including scaling.
This dot product is needed for evaluating moments of quantities that have not been yet
implemented in the lattice solver.
Args:
v (np.ndarray): Vector to dot.
out (np.ndarray): preallocated output array.
Returns:
np.ndarray: Dot product of H with v.
"""
out[:] = self.dot(v, out)
out[:] = (out - self.scaling.get_shift() * v) * self.scaling.get_scale()
return out
[docs]
def _set_scaling_and_cheby(
self,
lower_spectral_bound: float = -1,
upper_spectral_bound: float = 1,
delta: float = 0,
force_complex: bool = False,
) -> None:
"""Set the scaling and corresponding chebyshev objects.
Args:
lower_spectral_bound (float): Lower bound of the spectral range. Defaults to -1.
upper_spectral_bound (float): Upper bound of the spectral range. Defaults to 1.
delta (float): Spectral buffer for Hamiltonian scaling.
force_complex (bool): Flag for forcing complex chebyshev objects. Defaults to False.
"""
# simply set the scaling object with the provided parameters
self.scaling = lsc.chebyshev_scaling(lower_spectral_bound, upper_spectral_bound, delta)
# set the chebyshev object associated with the scaling object
self._set_cheby(force_complex)
[docs]
def _set_cheby(self, force_complex: bool = False) -> None:
"""Set the chebyshev object.
Args:
force_complex (bool): Flag for forcing complex chebyshev objects. Defaults to False.
"""
_dtype = "complex" if force_complex else self.dtype.__name__
# set the chebyshev object with the correct dtype
self.chebyshev = lsc.chebyshev(self.scaling, _dtype)
self.chebyshev.set_linear_operator(self.dot)
[docs]
def generateSelector(self, what: str = "all") -> np.ndarray: # noqa: ARG002
"""Trivial selector generator given a certain input configuration.
This method is just an empty placeholder, it will have to be specified for each
Hamiltonian implementation.
Args:
what (str): input configuration. Not used in this method.
Returns:
np.ndarray: array of 1's.
Warnings:
RuntimeWarning: warn all the time, since this function is always selecting everything.
"""
message = "Selector generation not implemented for current Hamiltonian. "
message += "All sites will be selected."
warnings.warn(message, category=RuntimeWarning, stacklevel=2)
return np.ones(self.dim, dtype=bool)
[docs]
def getBoundsLanczos(self, Order: int, tolerance: float = 1e-11) -> np.ndarray:
"""Get spectrum boundaries for the current Hamiltonian with the Lanczos method.
Args:
Order (int): Order of the Lanczos expansion
tolerance (float): Tolerance for the expansion. Defaults to 1e-11.
Returns:
np.ndarray: Boundaries of the spectrum.
"""
x0 = np.zeros(self.dim, dtype=self.dtype)
x0 = qt.generateRandomVector(self.generateSelector(), x0, self.dtype)
alphas, betas = qt.calcLanczos(self._dot_scaling, x0, Order, tolerance=tolerance)
evals = eigvalsh_tridiagonal(alphas, betas)
return evals[[0, -1]]
[docs]
def scaleH(
self,
delta: float = 0.1,
LanczOrder: Optional[int] = None,
energyBounds: Optional[Union[list, np.ndarray]] = None,
) -> None:
"""Rescale Hamiltonian to have spectrum between -1 and 1.
Args:
delta (float): Controls the extra compression of the band. This will be rescaled to be
(-1 + delta/2, 1 - delta/2). Defaults to 0.1.
LanczOrder (Optional[int]): Order of the Lanczos expansion used to determine the
bandwidth. Cannot be set together with energyBounds.
energyBounds (Optional[Union[list, np.ndarray]]): lower and upper energy boundaries
to determine scaling. Cannot be set together with LanczOrder.
Raises:
ValueError: If delta is not in (0, 2).
ValueError: If both LanczOrder and energyBounds are set.
"""
if delta <= 0 or delta >= 2:
raise ValueError("The scaling factor delta has to be in (0, 2)")
# the default behavior is an automatic scaling by estimating the energy boundaries through
# a Lanczos expansion. If the energy bounds are set, the LanczOrder should not be.
# Otherwise use automatic estimation of boundaries with the given LanczOrder (or the
# default 75)
if (energyBounds is not None) and (LanczOrder is not None):
message = "Both LanczOrder and energyBounds have been set. Select only one"
message += " to use the corresponding scaling method."
raise ValueError(message)
else:
# reset scaling before setting the new one
self.scaleH_inv()
# if the energy bounds are set, I don't have to do anything else but determine the
# scale otherwise estimate the boundaries through Lanczos expansion (at the given or
# default order)
if energyBounds is None:
# use default value if LanczOrder has not been set
if LanczOrder is None:
LanczOrder = 75
# estimate boundaries of the Hamiltonian via Lanczos expansion
evals = self.getBoundsLanczos(LanczOrder)
energyBounds = [evals[0], evals[-1]]
# set the scaling object
self._set_scaling_and_cheby(energyBounds[0], energyBounds[1], delta)
[docs]
def scaleH_inv(self) -> None:
"""Rescale Hamiltonian to have original spectrum."""
self._set_scaling_and_cheby()
[docs]
def scale2real(self, x: qt.fl_cmplx_or_ndarr) -> qt.fl_cmplx_or_ndarr:
"""Apply scaling to transform scaled frequencies to real ones.
Args:
x (qt.fl_cmplx_or_ndarr): array of scaled energies within (-1, 1).
Returns:
qt.fl_cmplx_or_ndarr: Real energies array.
"""
return qt._scale2real(x, self.scaling.get_scale(), self.scaling.get_shift())
[docs]
def parseDisDict(self, parse_dict: dict) -> None:
"""Parse disorder input dictionary.
The disorder dictionary specifies the properties of the disorder:
- type: kind of disorder, e.g. onsite (only type implemented for now).
- distribution: random distribution of disorder potential.
- sigma: standard deviation of the disorder distribution.
- avg: mean of distribution.
Note:
This will become obsolete with the introduction of a Schema for the disDict.
Args:
parse_dict (dict): Dictionary containing information about the disorder.
Raises:
ValueError: if the dictionary cannot be parsed
"""
self.disDict = {}
self.disDict["type"] = "onsite" if "type" not in parse_dict else parse_dict["type"]
self.disDict["distr"] = "uniform" if "distr" not in parse_dict else parse_dict["distr"]
self.disDict["sigma"] = 0.0 if "sigma" not in parse_dict else parse_dict["sigma"]
self.disDict["avg"] = 0.0 if "avg" not in parse_dict else parse_dict["avg"]
if self.disDict["type"] not in ["onsite"]:
raise ValueError("Invalid disorder type.")
if self.disDict["distr"] not in ["uniform", "gaussian", "laplace"]:
raise ValueError("Invalid disorder distribution.")
# NOTE: this might not be the best way to check.
if type(self.disDict["sigma"]) not in [float, int, np.float64]:
raise ValueError("Invalid disorder parameter 'sigma'.")
# NOTE: this might not be the best way to check.
if type(self.disDict["avg"]) not in [float, int, np.float64]:
raise ValueError("Invalid disorder parameter 'avg'.")
[docs]
def scrambleW(self) -> None:
"""Function that scrambles the disorder vector."""
# FIXME: this is not done properly. Fix once it's clear what one should use.
if self.disDict["type"] == "onsite":
if self.disDict["sigma"] == 0 and self.disDict["avg"] == 0:
self.disLocPot = np.zeros(self.dim, dtype=float)
elif self.disDict["distr"] == "uniform":
# np.sqrt(12) is the factor between the std of a uniform distribution and its width
self.disLocPot = (
np.sqrt(12) * self.disDict["sigma"] * (np.random.random(self.dim) - 0.5)
+ self.disDict["avg"]
)
elif self.disDict["distr"] == "gaussian":
self.disLocPot = np.random.normal(
loc=self.disDict["avg"], scale=self.disDict["sigma"], size=self.dim
)
self.disLocPot = np.clip(
self.disLocPot,
a_min=-5 * self.disDict["sigma"] + self.disDict["avg"],
a_max=5 * self.disDict["sigma"] + self.disDict["avg"],
)
elif self.disDict["distr"] == "laplace":
# 1/sqrt(2) is the factor between the std and the parameter λ of a laplace distr.
self.disLocPot = np.random.laplace(
loc=self.disDict["avg"],
scale=self.disDict["sigma"] / np.sqrt(2),
size=self.dim,
)
self.disLocPot = np.clip(
self.disLocPot,
a_min=-7 * self.disDict["sigma"] + self.disDict["avg"],
a_max=7 * self.disDict["sigma"] + self.disDict["avg"],
)
[docs]
def getTrResolvent(
self,
NSamples: int,
Order: int,
rescale: bool = True,
use_KPM: bool = True,
verbose: int = 2,
selector: Optional[np.ndarray] = None,
) -> Callable:
r"""Returns a function evaluating the trace of the resolvent.
This method computes the Chebyshev moments for a given order and given number of stochastic
samples to evaluate the trace. Can be used to evaluate the trace of the resolvent,
.. math::
\mathrm{Tr} \left[ \frac{1}{z - \hat{H}} \right]
for a given 'frequency' :math:`z` without the need to recompute the Chebyshev expansion.
This is useful if it is not a priori clear on how many frequencies the trace of the
resolvent needs to be evaluated.
The trace of the resolvent will be evaluated at the given frequencies. For real frequencies
the retarded resolvent is chosen.
NOTE: The trace is estimated as the *average* over all (normalized) stochastic
samples. Hence, multiply the result by the *linear* size of the Hamiltonian to get
the "true" estimate for the trace.
NOTE: If the number of samples is larger that the size of the selected subspace,
the trace is computed exactly (apart from the normalization, see above).
Args:
NSamples (int): Number of stochastic samples to approximate the trace.
Order (int): Order of the Chebyshev expanbsion to approximate the resolvent.
rescale (bool): Whether generated function transforms input into interval (-1, 1). \
Defaults to True.
use_KPM (bool): Flag controlling the use of KPM.
verbose (int): Controls verbosity of output and warnings. Defaults to 2.
selector (Optional[np.ndarray]): selects a sublattice of interest for the calculation.
Returns:
Callable: Function evaluating the trace of the resolvent for a given frequency.
"""
qio.print_header(
verbose,
"Chebyshev expansion of resolvent.",
{"Nsamples": NSamples, "Order": Order, "NSamplesDis": 1},
)
qio.warn_scale(self.scaling.get_scale(), self.scaling.get_shift())
if selector is None:
selector = self.generateSelector()
# initialize moments.
moments = np.zeros((Order + 1,), dtype=self.dtype)
tmp_moments = np.zeros((Order + 1,), dtype=self.dtype)
# get linear size of selected space
Nselected = selector.sum()
# stochastic calculation of moments for DOS
if Nselected <= NSamples:
if verbose > 0:
print("Calculating exact trace.")
# use exact trace
for rr, idx in enumerate(np.nonzero(selector)):
# setup exact trace in the natural basis of the Hamiltonian
x0 = np.zeros((selector.size), dtype=self.dtype)
x0[idx] = 1
# calculate and accumulate Chebyshev moments for the vector x0
self.chebyshev.calc_moments(Order + 1, bra=x0, ket=x0, res=tmp_moments)
moments += tmp_moments
# completion percentage printout
qio.print_progress(verbose, "norm", rr, Nselected)
# average calculated moments (for consistency with stochastic method)
moments /= Nselected
else:
if verbose > 0:
print("Calculating stochastic trace.")
# use stochastic trace
for rr in range(NSamples):
# initialize random vector
x0 = np.zeros(self.dim, dtype=self.dtype)
x0 = qt.generateRandomVector(selector, x0, dtype=self.dtype)
# calculate and accumulate Chebyshev moments for the random vector x0
self.chebyshev.calc_moments(Order + 1, bra=x0, ket=x0, res=tmp_moments)
moments += tmp_moments
# completion percentage printout
qio.print_progress(verbose, "norm", rr, NSamples)
# average calculated moments
moments /= NSamples
# set rescaling function
if rescale:
# apply rescaling
_z_rescale = self.scale2transformed
C = self.scaling.get_scale()
else:
# trivial function -> take bare input frequencies
def _id(z: qt.fl_cmplx_or_ndarr) -> qt.fl_cmplx_or_ndarr:
"""Trivial function.
Args:
z (qt.fl_cmplx_or_ndarr): input.
Returns:
qt.fl_cmplx_or_ndarr: input.
"""
return z
_z_rescale = _id
C = 1.0
# construct resolvent function using get_resolvent from qolossalTools module
def _get_tr_resolvent(z: qt.fl_cmplx_or_ndarr) -> qt.fl_cmplx_or_ndarr:
# apply rescaling function
z = _z_rescale(z)
# get function returning the trace of the resolvent
tr_resolvent = qt.get_tr_resolvent(moments, use_KPM=use_KPM)
return tr_resolvent(z) * C
return _get_tr_resolvent
[docs]
def _get_DOS_moments(
self,
NSamples: int,
Order: int,
NSamplesDis: int,
verbose: int,
selector: np.ndarray,
) -> np.ndarray:
"""Calculates the moments of the Chebyshev expansion of the density of states.
Only meant for internal and expert usage. This method does not implement any failsafe
mechanism or consistency checks.
Args:
NSamples (int): Number of stochastic samples for the calculation of the DoS.
Order (int): Order of the Chebyshev expansion used for the calculation
NSamplesDis (int): Number of disorder configurations to average over.
verbose (int): controls verbosity of output and warnings.
selector (np.ndarray): selects a sublattice of interest for the calculation.
Returns:
np.ndarray: Moments associated with the density of states
"""
# save original disorder configuration to be reset at the end
original_disLocPot = self.disLocPot.copy()
# initialize moments. These will be real even if the Hamiltonian is complex -- the DOS is
# a real quantity!
moments = np.zeros((Order + 1,))
tmp = np.zeros((Order + 1,), dtype=self.dtype)
# initialize random vector, I'll always use the same memory.
x0 = np.zeros(self.dim, dtype=self.dtype)
# stochastic calculation of moments
for dd in range(NSamplesDis):
qio.print_progress(verbose, "dis_head", dd, NSamplesDis)
for rr in range(NSamples):
# initialize random vector and normalize
x0 = qt.generateRandomVector(selector, x0, dtype=self.dtype)
# calculate and accumulate Chebyshev moments for the random vector x0 (pass it both
# as right and left vector)
self.chebyshev.calc_moments(Order + 1, bra=x0, ket=x0, res=tmp)
moments += tmp.real
# completion percentage printout
qio.print_progress(verbose - 1, "norm", rr, NSamples)
# randomize the disorder configuration for next cycle
self.scrambleW()
qio.print_progress(verbose, "dis", dd, NSamplesDis)
# average calculated moments
moments /= NSamples * NSamplesDis
# reset disorder to original configuration
self.disLocPot = original_disLocPot.copy()
return moments
[docs]
def getDOS(
self,
NSamples: int,
Order: int,
NOmegas: int = 101,
NSamplesDis: int = 1,
verbose: int = 2,
selector: Optional[np.ndarray] = None,
eWindow: Optional[Union[list, np.ndarray]] = None,
) -> tuple[np.ndarray, np.ndarray]:
r"""Calculates DoS through stochastic averaging of Chebyshev expansion.
This density of states integrates to 1, i.e. it is defined as:
.. math::
\rho(E) = \frac{1}{N} \sum\limits_i \delta(E - E_i) = \frac{1}{N} {\rm Tr}
\left[ \delta(E - \hat{H}) \right]
with :math:`E_i` the eigenvalues of the Hamiltonian and N the number of orbitals.
More details on the calculation are readily available in the literature. For example:
https://doi.org/10.1016/j.physrep.2020.12.001
Args:
NSamples (int): Number of stochastic samples for the calculation of the DoS.
Order (int): Order of the Chebyshev expansion used for the calculation
NOmegas (int): Number of frequencies to evaluate the DoS at. Defaults to 101.
NSamplesDis (int): Number of disorder configurations to average over.
verbose (int): controls verbosity of output and warnings.
selector (Optional[np.ndarray]): selects a sublattice of interest for the calculation.
eWindow (Optional[Union[list, np.ndarray]]): energy window for the result.
Returns:
tuple: DoS and corresponding frequencies.
"""
qio.print_header(
verbose,
"DoS",
{"NSamplesDis": NSamplesDis, "Nsamples": NSamples, "Order": Order, "Nomegas": NOmegas},
)
qio.warn_disorder(verbose, NSamplesDis, self.disDict)
qio.warn_scale(self.scaling.get_scale(), self.scaling.get_shift())
# TODO: the selector passed here should be a vector of ones and zeros. This means that the
# call to this function will most times have H.generateSelector() as the argument.
# Although extremely flexible, since one can pass custom selectors here, this is not
# particularly elegant. Once we have the need, one could pass only the string or expression
# to then call inside here the self.generateSelector with that string.
if selector is None:
selector = self.generateSelector()
# calculate moments
moments = self._get_DOS_moments(NSamples, Order, NSamplesDis, verbose, selector)
# get energies at which to evaluate the DOS.
energies = qt.get_energies_from_e_window(
scale=self.scaling.get_scale(),
shift=self.scaling.get_shift(),
n_omegas=NOmegas,
e_window=eWindow,
)
# get dos
DOS = qt.get_DOS_from_moments(
moments=moments,
frequencies=energies,
scale=self.scaling.get_scale(),
shift=self.scaling.get_shift(),
use_jackson_kernel=True,
)
# normalize by the selector integral
DOS *= selector.sum() / self.realDim
return DOS, energies
[docs]
def getGF(
self,
selector: Union[str, ArrayLike],
Order: int,
NOmegas: int,
NSamplesDis: int = 1,
verbose: int = 2,
eWindow: Optional[Union[list, np.ndarray]] = None,
) -> tuple[np.ndarray, np.ndarray]:
"""Calculate the retarded Green's Function in the real space basis.
This calculation is performed the same way as for the DOS, the only difference being the
states used for evaluating the moments: in this case we directly use the real space basis.
Args:
selector (Union[str, ArrayLike]): specifies what to compute ('diag', 'full', or array
of row-column pairs)
Order (int): Order of the Chebyshev expansion.
NOmegas (int): Number of frequencies to evaluate.
NSamplesDis (int): Number of disorder configurations to average over.
verbose (int): verbosity of routine.
eWindow (Optional[Union[list, np.ndarray]]): energy window for the result.
Returns:
tuple[np.ndarray, np.ndarray]: GF and corresponding frequencies.
"""
qio.print_header(
verbose,
"Green's function",
{"NSamplesDis": NSamplesDis, "Order": Order, "Nomegas": NOmegas},
)
qio.warn_disorder(verbose, NSamplesDis, self.disDict)
qio.warn_scale(self.scaling.get_scale(), self.scaling.get_shift())
# interpret selector to obtain the list of indices to loop through to calculate the GF
# use realDim to take care of the BCS case
idxList = qt.idxListFromSelector(selector, self.realDim)
lenList = idxList.shape[0]
# save original disorder configuration to be reset at the end
original_disLocPot = self.disLocPot.copy()
# get energies at which to evaluate the green's function
energies = qt.get_energies_from_e_window(
scale=self.scaling.get_scale(),
shift=self.scaling.get_shift(),
n_omegas=NOmegas,
e_window=eWindow,
)
# initialize GF and tmp array to store intermediate results
GF = np.zeros((lenList, NOmegas), dtype=complex)
GF_tmp = np.zeros((NOmegas), dtype=complex)
# stochastic calculation of moments
for dd in range(NSamplesDis):
qio.print_progress(verbose, "dis_head", dd, NSamplesDis)
# loop over selected indices
for ii, idx in enumerate(idxList):
# initialize vectors for the expansion
vL = np.zeros(self.dim)
vR = np.zeros(self.dim)
# idx contains the indices of the GF to be calculated.
vL[idx[0]] = 1
vR[idx[1]] = 1
# simply call the chebyshev object resolvent gf with the selected vector
# TODO: add the control over eta and jackson kernel
_ = self.chebyshev.resolvent_greens_function(
energies,
1e-9,
Order + 1,
vL,
vR,
GF_tmp,
resolvent_type="+",
advanced=False,
filter="jackson",
num_smooth=0,
)
# we need to accumulate the result, as for different disorder configurations, we
# need to re-calculate the GF
GF[ii] += GF_tmp
qio.print_progress(verbose - 1, "norm", ii, lenList)
# randomize the disorder configuration for next cycle
self.scrambleW()
qio.print_progress(verbose, "dis", dd, NSamplesDis)
GF /= NSamplesDis
# reset disorder to original configuration
self.disLocPot = original_disLocPot.copy()
# reshape output if the full GF has been calculated since it will only contain the upper
# triangular part of the GF matrix
if isinstance(selector, str) and selector == "full":
GF = GF.reshape(self.realDim, self.realDim, NOmegas)
return GF, energies
[docs]
def _get_conductivity_moments(
self,
tensor_element: str,
NSamples: int,
Order: int,
NSamplesDis: int,
verbose: int,
selector: np.ndarray,
) -> np.ndarray:
"""Calculates the moments of the Chebyshev expansion of the conductivity.
Only meant for internal and expert usage. This method does not implement any failsafe
mechanism or consistency checks.
Args:
tensor_element (str): specifies which element of the conductivity tensor to calculate,
as literal, e.g. 'xx'.
NSamples (int): Number of stochastic samples for the calculation of the conductivity.
Order (int): Order of the Chebyshev expansion used for the calculation.
NSamplesDis (int): Number of disorder configurations to average over.
verbose (int): verbosity of routine. Defaults to 1.
selector (Optional[np.ndarray]): selects a sublattice of interest for the calculation.
Raises:
ValueError: if the velocity operator is not set for the current object.
Returns:
np.ndarray: Moments associated with the conductivity
"""
if self.velocity_operator is None:
raise ValueError(qio._no_velocity_operator_error_message)
# save original disorder configuration to be reset at the end
original_disLocPot = self.disLocPot.copy()
# obtain correct index for velocity selection
directionL, directionR = qt.dir2idx(tensor_element)
# initialize moments and allocate random vector
Cmn: np.ndarray = np.zeros((Order + 1, Order + 1), dtype=complex)
tmp_Cmn: np.ndarray = np.zeros((Order + 1, Order + 1), dtype=complex)
# force vectors to be complex -- if H is real, velocity is imaginary, so it's never a
# purely real case!
xR = np.zeros(self.dim, dtype=complex)
# also force the cheby object to be complex!
self._set_cheby(force_complex=True)
# stochastic calculation of moments
for dd in range(NSamplesDis):
qio.print_progress(verbose, "dis_head", dd, NSamplesDis)
for rr in range(NSamples):
# initialize random vector and normalize it
xR = qt.generateRandomVector(selector, xR, complex)
xL = np.zeros_like(xR)
xL = self.velocity_operator[directionL](xR, xL)
# calculate Chebyshev moments of the conductivity. Passing (V^† @ x0) = (V @ x0)
# as left vector since it has to be the ket form (since the Hamiltonian will be
# applied from the left). The routine will take care of conjugating at the end.
# Because this will always return complex results (see comment for xR), explicitly
# take it astype(dtype)
self.chebyshev.double_expansion_fast(
n_cheby_moments=Order + 1,
linear_operator_V=self.velocity_operator[directionR],
bra=xL,
ket=xR,
res=tmp_Cmn,
)
Cmn += tmp_Cmn
qio.print_progress(verbose - 1, "norm", rr, NSamples)
# randomize the disorder configuration for next cycle
self.scrambleW()
qio.print_progress(verbose, "dis", dd, NSamplesDis)
# average the moments
Cmn = Cmn / (NSamples * NSamplesDis)
# reset disorder to original configuration
self.disLocPot = original_disLocPot.copy()
# reset the chebyshev object to the original dtype
self._set_cheby(force_complex=False)
to_return = Cmn if self.dtype is complex else Cmn.real
return to_return
[docs]
def getConductivity(
self,
tensor_element: str,
T: float,
NSamples: int,
Order: int,
NOmegas: int = 101,
NSamplesDis: int = 1,
verbose: int = 2,
selector: Optional[np.ndarray] = None,
eWindow: Optional[Union[list, np.ndarray]] = None,
) -> tuple[np.ndarray, np.ndarray]:
r"""Calculates the conductivity of the system from the Kubo-Greenwood formula.
This routine follows the method presented in
https://doi.org/10.1016/j.physrep.2020.12.001
In short: a Chebyshev expansion of the zero temperature Kubo formula is performed. This
translates into the expansion of the function :math:`\delta(E - \hat{H}) \hat{V}_{\alpha}
\delta(E - \hat{H})` :math:`\hat{V}_{\beta}`, where :math:`E` is the energy,
:math:`\hat{H}` is the Hamiltonian and :math:`\hat{V}` is the velocity operator in the
direction :math:`\alpha\,{\rm and }\,\beta`. The conductivity reads
.. math::
\sigma(E) = \frac{\pi \hbar e^2}{\Omega} {\rm Tr}\left[ \delta(E - \hat{H})
\hat{V}_{\alpha} \delta(E - \hat{H}) \right]
The temperature dependent result is obtained by simply convolving the zero temperature
result by a 'Fermi window', i.e. minus the derivative of the Fermi function with respect to
the energy.
NOTE: the output of this routine is in units of :math:`e^2/(h Å^{(D-2)})` with D the
dimensionality of the system and does NOT include any spin- degeneracy coefficients for
spinless systems.
Args:
tensor_element (str): specifies which element of the conductivity tensor to calculate,
as literal, e.g. 'xx'.
T (float): Temperature in Kelvins.
NSamples (int): Number of stochastic samples for the calculation of the conductivity.
Order (int): Order of the Chebyshev expansion used for the calculation
NOmegas (int): Number of energies to evaluate the conductivity at. Defaults to 101.
NSamplesDis (int): Number of disorder configurations to average over.
verbose (int): verbosity of routine. Defaults to 1.
selector (Optional[np.ndarray]): selects a sublattice of interest for the calculation.
eWindow (Optional[Union[list, np.ndarray]]): energy window for the result.
Returns:
tuple: Conductivity and corresponding frequencies.
"""
qio.print_header(
verbose,
"conductivity",
{"NSamplesDis": NSamplesDis, "Nsamples": NSamples, "Order": Order, "Nomegas": NOmegas},
)
qio.warn_disorder(verbose, NSamplesDis, self.disDict)
qio.warn_scale(self.scaling.get_scale(), self.scaling.get_shift())
if selector is None:
selector = self.generateSelector()
# get the moments associated with the conductivity
mu_mn = self._get_conductivity_moments(
tensor_element, NSamples, Order, NSamplesDis, verbose, selector
)
# get chemical potentials at which to evaluate the conductivity
chemical_potentials = qt.get_energies_from_e_window(
scale=self.scaling.get_scale(),
shift=self.scaling.get_shift(),
n_omegas=NOmegas,
e_window=eWindow,
)
# evaluate the zero temperature Kubo-Greenwood conductivity given the moments.
conductivity = qt.get_KG_conductivity_from_moments(
moments=mu_mn,
chemical_potentials=chemical_potentials,
scale=self.scaling.get_scale(),
shift=self.scaling.get_shift(),
total_volume=self.total_volume,
system_linear_size=self.realDim,
use_jackson_kernel=True,
)
# get the temperature dependent result from the zero temperature one
conductivity = qt.get_T_dependent_conductivity(conductivity, chemical_potentials, T)
# cast to real if the system is real, this makes things simpler since in real cases, the
# imaginary part could still be different from 0 for convergence reasons.
conductivity = conductivity if self.dtype is complex else conductivity.real
return conductivity, chemical_potentials
[docs]
def getConductivityMu(
self,
tensor_element: str,
NSamples: int,
Order: int,
mu: Union[float, np.ndarray],
NSamplesDis: int = 1,
verbose: int = 2,
selector: Optional[np.ndarray] = None,
) -> Union[float, complex, np.ndarray]:
r"""Calculates the conductivity of the system defined by H at a single chemical potential.
Given a Hamiltonian, the corresponding conductivity can be calculated as
.. math::
\sigma(\mu) = \frac{\pi \Omega}{\Delta E^2} \braket{\phi^L|\phi^R}
where :math:`\Omega` is the volume of the system, :math:`\Delta E` is the scaling applied
to the Hamiltonian to ensure the spectrum lies within (-1, 1) and the vectors :math:`\phi`
are
.. math::
\ket{\phi^L} = \sum_n Im[G^+_n(\mu)] T_n(\tilde{H}) V \ket{r}
\ket{\phi^R} = \sum_n Im[G^+_n(\mu)] V T_n(\tilde{H}) \ket{r}
where :math:`\ket{r}` is the random vector, :math:`V` is the velocity operator,
:math:`G^+_n` is the n-th GF coefficient and :math:`T_n(\tilde{H})` is the m-th Chebyshev
polynomial.
NOTE: the output of this routine is in units of :math:`e^2/(h Å^{(D-2)})` with D the
dimensionality of the system and does NOT include any spin- degeneracy coefficients for
spinless systems.
Args:
tensor_element (str): specifies which element of the conductivity tensor to calculate,
as literal, e.g. 'xx'.
NSamples (int): Number of stochastic samples for the calculation of the conductivity.
Order (int): Order of the Chebyshev expansion used for the calculation
mu (Union[float, np.ndarray]): Energy(ies) (as in chemical potential) at which to
evaluate the conductivity.
NSamplesDis (int): Number of disorder configurations to average over.
verbose (int): verbosity of routine. Defaults to 1.
selector (Optional[np.ndarray]): selects a sublattice of interest for the DOS.
Raises:
ValueError: if the velocity operator is not set for the current object.
Returns:
qt.fl_cmplx_or_ndarr: Conductivity at the specified chemical potential(s).
"""
qio.print_header(
verbose,
"conductivity",
{"NSamplesDis": NSamplesDis, "Nsamples": NSamples, "Order": Order, "mu": mu},
)
qio.warn_disorder(verbose, NSamplesDis, self.disDict)
qio.warn_scale(self.scaling.get_scale(), self.scaling.get_shift())
if self.velocity_operator is None:
raise ValueError(qio._no_velocity_operator_error_message)
if selector is None:
selector = self.generateSelector()
# obtain correct index for velocity selection
directionL, directionR = qt.dir2idx(tensor_element)
# save original disorder configuration to be reset at the end
original_disLocPot = self.disLocPot.copy()
# mu can be a single value or multiple, so make mu at least 1D. The gm coefficients
# will then have the correct shape, compatible with the moments calculating method.
mus = np.atleast_1d(mu)
# calculate prefactors (fixed 1e-9 eta for now since I always apply the kernel)
gm = lsc.calc_prefactors(self.scaling, mus, 1e-9, Order + 1)
# take the imaginary part and transpose, to not have to change too much afterwards.
# NOTE: this will anyway change since the routine for calculating the accumulated double
# moments will be included in the chebyshev module
gm = lsc.apply_jackson(gm).imag.T
# initialize conductivity. Needs shape (len(mus),) to be compatible with cheby method
cond = np.zeros(mus.shape, dtype=complex)
# allocate random vector
xR = np.zeros(self.dim, dtype=complex)
# stochastic calculation of moments
for dd in range(NSamplesDis):
qio.print_progress(verbose, "dis_head", dd, NSamplesDis)
for rr in range(NSamples):
# initialize random vector
xR = qt.generateRandomVector(selector, xR, complex)
xL = np.zeros_like(xR)
xL = self.velocity_operator[directionL](xR, xL)
# calculate Chebyshev moments of the conductivity. Passing (V^† @ x0) = (V @ x0)
# as left vector since it has to be the ket form (since the Hamiltonian will be
# applied from the left). The routine will take care of conjugating at the end.
cond += qt.ChebyDoubleMomentsAccumulate(
self._dot_scaling,
gm,
xL,
xR,
self.velocity_operator[directionR],
Order,
dtype=complex,
)
qio.print_progress(verbose - 1, "norm", rr, NSamples)
# randomize the disorder configuration for next cycle
self.scrambleW()
qio.print_progress(verbose, "dis", dd, NSamplesDis)
# average conductivity contributions
cond = cond / (NSamples * NSamplesDis)
# reset disorder to original configuration
self.disLocPot = original_disLocPot.copy()
# normalize result. See qolossal_tools::get_KG_conductivity_from_moments for more details.
cond *= 2 * self.scaling.get_scale() ** 2 / (self.total_volume / self.realDim)
# cast to real if the system is real, this makes things simpler since in real cases, the
# imaginary part could still be different from 0 for convergence reasons.
to_return = cond if self.dtype is complex else cond.real
# make sure I return a single value (not an array) if a single value was requested
if isinstance(mu, float):
to_return = to_return[0]
return to_return
[docs]
def getConductivity_KB(
self,
tensor_element: str,
T: float,
NSamples: int,
Order: int,
NOmegas: int = 101,
NSamplesDis: int = 1,
verbose: int = 2,
selector: Optional[np.ndarray] = None,
eWindow: Optional[Union[list, np.ndarray]] = None,
NOmegasInt: Optional[int] = None,
) -> tuple[np.ndarray, np.ndarray]:
r"""Calculates the conductivity of the system evaluating the Kubo-Bastin formula.
This routine follows the method presented in
https://doi.org/10.1103/PhysRevLett.114.116602
In short: a Chebyshev expansion of the Kubo-Bastin formula is performed. This translates
into the expansion of the function:
.. math::
v_\alpha \delta(\varepsilon - \hat{H}) \hat{v}_\beta \frac{dG^+(\varepsilon)}
{d\varepsilon} - \hat{v}_\alpha \frac{dG^-(\varepsilon)}{d\varepsilon} \hat{v}_\beta
\delta(\varepsilon - \hat{H})
where :math:`\varepsilon` is the energy, :math:`\hat{H}` is the Hamiltonian,
:math:`\hat{v}` is the velocity operator and :math:`\alpha, \beta` indicate the direction.
The modus operandi is exactly the same as for the other conductivity calculations based on
the Chebyshev expansion.
NOTE: the output of this routine is in units of :math:`e^2/(h Å^{(D-2)})` with D the
dimensionality of the system and does NOT include any spin- degeneracy coefficients for
spinless systems.
Args:
tensor_element (str): specifies which element of the conductivity tensor to calculate,
as literal, e.g. 'xx'.
T (float): Temperature in Kelvins.
NSamples (int): Number of stochastic samples for the calculation of the conductivity.
Order (int): Order of the Chebyshev expansion used for the calculation
NOmegas (int): Number of energies to evaluate the conductivity at. Defaults to 101.
NSamplesDis (int): Number of disorder configurations to average over.
verbose (int): verbosity of routine. Defaults to 1.
selector (Optional[np.ndarray]): selects a sublattice of interest for the calculation.
eWindow (Optional[Union[list, np.ndarray]]): energy window for the result.
NOmegasInt (Optional[int]): number of frequencies to use to calculate the E integral.
Raises:
ValueError: if the velocity operator is not set for the current object.
Returns:
tuple: Conductivity and corresponding frequencies.
"""
qio.print_header(
verbose,
"Kubo-Bastin conductivity",
{
"NSamplesDis": NSamplesDis,
"Nsamples": NSamples,
"Order": Order,
"NOmegas": NOmegas,
"Temperature (K)": T,
},
)
qio.warn_disorder(verbose, NSamplesDis, self.disDict)
qio.warn_scale(self.scaling.get_scale(), self.scaling.get_shift())
if self.velocity_operator is None:
message = "The velocity operator has not been set. This calculation can therefore not "
message += "be performed."
raise ValueError(message)
if selector is None:
selector = self.generateSelector()
# get the moments associated with the conductivity
mu_mn = self._get_conductivity_moments(
tensor_element, NSamples, Order, NSamplesDis, verbose, selector
)
# get the chemical potentials at which to calculate the conductivity
chemical_potentials = qt.get_energies_from_e_window(
scale=self.scaling.get_scale(),
shift=self.scaling.get_shift(),
n_omegas=NOmegas,
e_window=eWindow,
)
# default NOmegasInt to NOmegas if not provided
NOmegasInt = NOmegas if NOmegasInt is None else NOmegasInt
# get the kubo-bastin conductivity from the moments
conductivity = qt.get_KB_conductivity_from_moments(
mu_mn,
chemical_potentials,
self.scaling.get_scale(),
self.scaling.get_shift(),
self.total_volume,
self.realDim,
T,
NOmegasInt,
)
# cast to real if the system is real, this makes things simpler since in real cases, the
# imaginary part could still be different from 0 for convergence reasons.
conductivity = conductivity if self.dtype is complex else conductivity.real
return conductivity, chemical_potentials
[docs]
def _get_expectation_value_moments(
self,
A: Callable,
expectation_value_dtype: type,
NSamples: int,
Order: int,
NSamplesDis: int,
verbose: int,
selector: np.ndarray,
) -> np.ndarray:
"""Calculates the moments of the Chebyshev expansion of an operator's expectation value.
Only meant for internal and expert usage. This method does not implement any failsafe
mechanism or consistency checks.
Args:
A (Callable): Linear operator.
expectation_value_dtype (type): dtype of the linear operator.
NSamples (int): Number of stochastic samples.
Order (int): Order of the Chebyshev expansion used for the calculation
NSamplesDis (int): Number of disorder configurations to average over.
verbose (int): controls verbosity of output and warnings.
selector (np.ndarray): selects a sublattice of interest for the calculation.
Returns:
np.ndarray: Moments associated with the expectation value of operator A.
"""
# save original disorder configuration to be reset at the end
original_disLocPot = self.disLocPot.copy()
is_anything_complex = complex in [expectation_value_dtype, self.dtype]
moments_dtype = complex if is_anything_complex else float
self._set_cheby(force_complex=is_anything_complex)
# initialize moments to the correct dtype
moments: np.ndarray = np.zeros((Order + 1,), dtype=expectation_value_dtype)
tmp_moments: np.ndarray = np.zeros((Order + 1,), dtype=moments_dtype)
xL = np.zeros(self.dim, dtype=self.dtype)
# stochastic calculation of moments
for dd in range(NSamplesDis):
qio.print_progress(verbose, "dis_head", dd, NSamplesDis)
for rr in range(NSamples):
# initialize random vector
xL = qt.generateRandomVector(selector, xL, expectation_value_dtype)
xR = np.zeros_like(xL)
xR = A(xL, xR)
# calculate and accumulate Chebyshev moments
self.chebyshev.calc_moments(Order + 1, bra=xL, ket=xR, res=tmp_moments)
moments += tmp_moments.astype(expectation_value_dtype)
# completion percentage printout
qio.print_progress(verbose - 1, "norm", rr, NSamples)
# randomize the disorder configuration for next cycle
self.scrambleW()
qio.print_progress(verbose, "dis", dd, NSamplesDis)
# average calculated moments
moments /= NSamples * NSamplesDis
# reset disorder to original configuration and reset the cheby object
self.disLocPot = original_disLocPot.copy()
self._set_cheby(force_complex=False)
return moments
[docs]
def getExpectationValue(
self,
A: Callable,
dtype: type,
T: float,
mu: float,
NSamples: int,
Order: int,
NPadeFreq: int = 101,
NSamplesDis: int = 1,
verbose: int = 2,
selector: Optional[np.ndarray] = None,
) -> Union[float, complex]:
r"""Calculates the expectation value of an operator through Chebyshev expansion.
This is done by performing the Chebyshev expansion of
.. math::
\langle \hat{A} \rangle = {\rm Tr} \left[ \hat{A} f(\hat{H}) \right]
where :math:`\hat{A}` is the operator, :math:`\hat{H}` is the hamiltonian and f is the
fermi function. This expression is expanded in terms of Chebyshev polynomials. The integral
over the energies is performed via the resummation of the Fermi-Pade' expansion of the
Fermi function.
Args:
A (Callable): Linear operator.
dtype (type): dtype of the linear operator.
T (float): temperature in Kelvins. Must be greater than 0.
mu (float): chemical potential.
NSamples (int): Number of stochastic samples.
Order (int): Order of the Chebyshev expansion used for the calculation
NPadeFreq (int): Number of Pade frequencies to evaluate the approximation to
the Fermi function.
NSamplesDis (int): Number of disorder configurations to average over.
verbose (int): controls verbosity of output and warnings.
selector (Optional[np.ndarray]): selects a sublattice of interest for the calculation.
Raises:
ValueError: if temperature is too small.
Returns:
Union[float, complex]: expectation value of operator A.
"""
if T < 1e-10:
raise ValueError("Temperature too low, the Fermi-Pade method breaks down for T -> 0.")
qio.print_header(
verbose,
"expectation value of provided operator",
{
"NSamplesDis": NSamplesDis,
"Temperature (K)": T,
"Nsamples": NSamples,
"Order": Order,
"NPadeFreq": NPadeFreq,
},
)
qio.warn_disorder(verbose, NSamplesDis, self.disDict)
qio.warn_scale(self.scaling.get_scale(), self.scaling.get_shift())
if selector is None:
selector = self.generateSelector()
# use complex if either the Hamiltonian or the Operator are complex
expV_dtype = complex if ((self.dtype is complex) or (dtype is complex)) else float
# get moments associated with the expectation value of operator A
moments = self._get_expectation_value_moments(
A=A,
expectation_value_dtype=expV_dtype,
NSamples=NSamples,
Order=Order,
NSamplesDis=NSamplesDis,
verbose=verbose,
selector=selector,
)
# get expectation value from the moments. normalize by the selector to scale it up to the
# linear size of the system.
expV = selector.sum() * qt.get_expectation_value_from_moments(
moments=moments,
T=T,
mu=mu,
scale=self.scaling.get_scale(),
shift=self.scaling.get_shift(),
n_pade_freq=NPadeFreq,
)
# cast the value to float in the corresponding case, since the gm contributions will
# always be complex
toreturn = expV if expV_dtype is complex else expV.real
return toreturn
[docs]
def getConductivityOptical(
self,
tensor_element: str,
T: float,
mu: float,
NSamples: int,
Order: int,
NOmegas: int = 101,
NSamplesDis: int = 1,
verbose: int = 2,
selector: Optional[np.ndarray] = None,
eWindow: Optional[Union[list, np.ndarray]] = None,
NOmegasInt: Optional[int] = None,
) -> tuple[np.ndarray, np.ndarray]:
r"""Calculates the complex optical conductivity of the system.
In short: a Chebyshev expansion of the Kubo formula is performed. The procedure is similar
to the one presented in
https://journals.aps.org/prb/abstract/10.1103/PhysRevB.94.235405
Details on the routine can be found in the internal notes. In particular, eq. 32 of the
notes is implemented.
NOTE: the output of this routine is in units of :math:`e^2/(h Å^{(D-2)})` with D the
dimensionality of the system and does NOT include any spin- degeneracy coefficients for
spinless systems.
Args:
tensor_element (str): specifies which element of the conductivity tensor to calculate,
as literal, e.g. 'xx'.
T (float): Temperature in Kelvins.
mu (float): Chemical potential.
NSamples (int): Number of stochastic samples for the calculation of the conductivity.
Order (int): Order of the Chebyshev expansion used for the calculation
NOmegas (int): Number of energies to evaluate the conductivity at. Defaults to 101.
NSamplesDis (int): Number of disorder configurations to average over.
verbose (int): verbosity of routine. Defaults to 1.
selector (Optional[np.ndarray]): selects a sublattice of interest for the calculation.
eWindow (Optional[Union[list, np.ndarray]]): energy window for the result.
NOmegasInt (Optional[int]): number of frequencies to use to calculate the E integral.
Returns:
tuple: Conductivity and corresponding frequencies.
"""
qio.print_header(
verbose,
"optical conductivity",
{
"Nsamples": NSamples,
"Order": Order,
"NOmegas": NOmegas,
"mu": mu,
"Temperature (K)": T,
"NSamplesDis": NSamplesDis,
},
)
qio.warn_disorder(verbose, NSamplesDis, self.disDict)
qio.warn_scale(self.scaling.get_scale(), self.scaling.get_shift())
if selector is None:
selector = self.generateSelector()
# default NOmegasInt to NOmegas if not provided
NOmegasInt = NOmegas if NOmegasInt is None else NOmegasInt
# get energies at which to evaluate the optical conductivity
# NOTE: eWindow specifies the frequencies of the perturbation! meaning it's not related to
# any system's property, hence only the scale is relevant for the transformation, not the
# shift!
energies = qt.get_energies_from_e_window(
scale=self.scaling.get_scale(), shift=0, n_omegas=NOmegas, e_window=eWindow
)
# get the moments associated with the conductivity
moments = self._get_conductivity_moments(
tensor_element, NSamples, Order, NSamplesDis, verbose, selector
)
# expand the moments to obtain the optical conductivity
optical_conductivity = qt.get_optical_conductivity_from_moments(
moments=moments,
frequencies=energies,
scale=self.scaling.get_scale(),
shift=self.scaling.get_shift(),
total_volume=self.total_volume,
system_linear_size=self.realDim,
T=T,
chemical_potential=mu,
NOmegasInt=NOmegasInt,
verbose=verbose,
)
return optical_conductivity, energies