"""Abstract Hamiltonian class interfacing directly to the Lattice Builder.
Copyright © 2019-2023 HQS Quantum Simulations GmbH. All Rights Reserved.
"""
from typing import Optional, Union, cast, Callable
import warnings
from abc import abstractmethod
from lattice_builder import Builder as LB
from lattice_validator import Validator as LV
from lattice_solver.chebyshev import calc_prefactors
from crystal_tools import bands
from .abstractHamiltonian import Hamiltonian
import qolossal.qolossal_tools as qt
import qolossal.qolossal_io as qio
import numpy as np
from scipy.sparse import csr_array
[docs]
class HamiltonianLB(Hamiltonian):
r"""Abstract Hamiltonian class that takes care of the interface with the LatticeBuilder."""
[docs]
def __init__(self, input_dict: dict, disDict: Optional[dict] = None) -> None:
r"""Initialization of the class.
This init takes care of validation and building of the system coresponding to the input
dictionary and will extract the general information needed by most Hamiltonians. The
Hamiltonian itself will have to be extracted by each the specific implementation. Note that
this class is still abstract as no dot product is defined.
Args:
input_dict (dict): Input dictionary for the Lattice Builder to construct the
Hamiltonian
disDict (Optional[dict]): dictionary specifying the disorder.
Raises:
ValueError: if the input_dict or dis dictionaries cannot be validated.
"""
lv = LV(profile="unlimited") # selection of qolossal profile
is_valid, dict_input_cluster = lv.validateConfiguration(input_dict)
if not is_valid:
raise ValueError("Lattice input can not be validated: " + str(dict_input_cluster))
# build
self.lb: LB = LB(dict_input_cluster)
# extract complex, spin and BCS flags. Keep a copy of them since, in the abstract class,
# they are independent of whether we build the Hamiltonian with the LB or in any other way
self.is_complex = self.lb.is_complex
self.is_spin_matrix = self.lb.is_spin_matrix
self.is_BCS = self.lb.is_BCS
# NOTE: this is only needed for the magnetic susceptibility calculation. Can be dropped
# out if AutoDiff is used or if linearResp superseeds getMagneticSusc
self.is_spinflip: bool
# ...and number of repetition of uc in each direction
self.Lx, self.Ly, self.Lz = self.lb.size_uc
# the distinction between dim and realDim
# get the dimension of the system.
self.realDim = np.prod(self.lb.size_uc) * self.lb.M_uc * (1 + self.is_spin_matrix)
# dim gets an extra factor of 2 if it's BCS.
self.dim = self.realDim * (1 + self.is_BCS)
# get the velocity, build the correct matrix depending on the case, then build the linear
# operator associated with it.
_V = (self.lb.get_velocity0(0), self.lb.get_velocity0(1), self.lb.get_velocity0(2))
# handle the velocity appropriately given the type of system (spinful and/or complex)
if self.is_spin_matrix:
# the velocity computed by the lb in the spin_matrix case will be a tuple for each
# direction, so it needs to be "expanded" for each tupe element.
_V = (qt.blocks2full(_V[0]), qt.blocks2full(_V[1]), qt.blocks2full(_V[2]))
else:
# in this case the input dictionary specifies the system as spinful, BUT the up and
# down spins are perfectly equivalent, so there is no need to keep both.
if isinstance(_V[0], tuple):
_V = (_V[0][0], _V[1][0], _V[2][0])
if self.is_BCS:
# if BCS simply repeat the velocity operator for the diagonal blocks (which correspond
# to c^dag c and c c^dag and leave an empty matrix for the off-diagonals)
_V = (
qt.blocks2full((_V[0], _V[0], csr_array(_V[0].shape))),
qt.blocks2full((_V[1], _V[1], csr_array(_V[1].shape))),
qt.blocks2full((_V[2], _V[2], csr_array(_V[2].shape))),
)
# get linear operators starting from the matrices
self.velocity_operator = (
qt.matrix_to_linop(_V[0]),
qt.matrix_to_linop(_V[1]),
qt.matrix_to_linop(_V[2]),
)
# get volume of unit cell intended as length if system is 1D, area if 2D and volume if 3D
lat_vecs = self.lb.config["unitcell"]["lattice_vectors"]
# in 1D the volume will be the length of this vector
tmpVolume = np.array(lat_vecs[0])
# in 2D it's |v1 x v2|
if len(lat_vecs) > 1:
tmpVolume = np.cross(tmpVolume, lat_vecs[1])
# in 3D it's |(v1 x v2) * v3|
if len(lat_vecs) > 2:
tmpVolume = np.dot(tmpVolume, lat_vecs[2])
# take the absolute value of the resulting vector (in 3D the result is a number, but the
# abs() must be taken anyways)
self.total_volume = np.sqrt(np.sum(tmpVolume**2)) * np.prod(self.lb.size_uc)
# FIXME: this is not yet general. Fix once it's clear what one should use.
# setup on-site disorder
if disDict is None:
disDict = {}
self.parseDisDict(disDict)
self.scrambleW()
# set the scaling
self._set_scaling_and_cheby()
[docs]
@abstractmethod
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.
Raises:
NotImplementedError: If the method is not overloaded in the children classes.
"""
raise NotImplementedError("The set_Hamiltonian function has not been defined.")
[docs]
def getH(self) -> csr_array:
"""Getter function for Hamiltonian in sparse format.
Returns:
csr_array: Hamiltonian matrix in sparse row format.
Raises:
RuntimeError: if the Hamiltonian is not available.
"""
tmp = self.lb.getH0()
if tmp is None:
raise RuntimeError("Hamiltonian not available. The system was not properly set up.")
if isinstance(tmp, tuple):
return qt.blocks2full((tmp[0], tmp[1], tmp[2]))
else:
return tmp
[docs]
def todense(self) -> np.ndarray:
"""Getter function for Hamiltonian in dense format.
Returns:
np.ndarray: Hamiltonian matrix in dense format.
"""
return self.getH().toarray()
[docs]
def generateSelector(self, what: str = "all") -> np.ndarray:
"""Generate selector given a certain input configuration.
Args:
what (str): input configuration.
Returns:
np.ndarray: array composed of 1's and 0's selecting the correct sublattice.
Raises:
ValueError: if the argument 'what' is invalid.
Warnings:
RuntimeWarning: if the selection is inconsistent with the system and a default choice
is made.
"""
# FIXME: this ONLY works for the sparse and dense implementations. But these are also the
# only ones to have a spinful implementation. Nonetheless this is not particularly safe.
# In general, every Hamiltonian will have to provide its selector, since the way different
# spin channels are taken into account depends on the implementation.
vec = np.zeros(self.dim, dtype=bool)
# will be all 1's if realDim == self.dim otherwise only the first half will be ones
vec[: self.realDim] = 1
if what in ["up", "down"]:
if self.is_spin_matrix:
if what == "up":
vec[self.realDim // 2 :] = 0
else:
vec[: self.realDim // 2] = 0
else:
message = "The system is spinless, selecting up or down spins will have no effect."
warnings.warn(message, category=RuntimeWarning, stacklevel=2)
elif what != "all":
raise ValueError(f"Invalid value {what!r} for argument 'what'.")
return vec
[docs]
def getSpectral(
self,
num_kpoints: float,
Order: int,
NOmegas: int,
NSamplesDis: int = 1,
verbose: int = 2,
eWindow: Optional[Union[list, np.ndarray]] = None,
) -> tuple:
r"""Calculate the spectral function.
This calculation is performed the same way as for the GF, the only difference being the
states used for evaluating the moments: in this case we directly use the Bloch functions.
The band structure is then simply :math:`-\operatorname{Im}[G_{kk}]/\pi` for k's along the
correct k path.
Args:
num_kpoints (float): number of kpoints for the spectral function calculation.
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: spectral function, frequencies and corresponding k-path information.
"""
qio.print_header(
verbose,
"spectral function",
{
"NSamplesDis": NSamplesDis,
"Order": Order,
"Nomegas": NOmegas,
"Number of unique k points": num_kpoints,
},
)
qio.warn_disorder(verbose, NSamplesDis, self.disDict)
qio.warn_scale(self.scaling.get_scale(), self.scaling.get_shift())
selector = self.generateSelector()
# save original disorder configuration to be reset at the end
original_disLocPot = self.disLocPot.copy()
# get energies given ewindow
energies = qt.get_energies_from_e_window(
scale=self.scaling.get_scale(),
shift=self.scaling.get_shift(),
n_omegas=NOmegas,
e_window=eWindow,
)
# save user-defined bc
user_bc = self.lb.BC
# change the bc to periodic for non-trivial directions only (dimension = len(lat_vec))
# NOTE: the strategy for this has to be streamlined, as it appears in the LB, too.
lat_vec = self.lb.config["unitcell"]["lattice_vectors"]
dimensions = len(lat_vec)
tmp = [1 if ii < dimensions else 0 for ii in range(3)]
self.lb.BC = tuple(tmp)
# get information needed to compute the Points along the correct kPath
kDict, kPath, primitiveVec, *_ = bands.returnLatticeAndPath(
self.lb.config["unitcell"]["lattice_vectors"]
)
# get points along standard kpath and corresponding special points)
kPathPoints, _, locationSpecialPoints, indexSpecialPoints = bands.constructPath(
kDict, kPath, primitiveVec, totalSamplePoints=num_kpoints
)
# the total number of kpoints can differ from num_kpoints since one could pass through
# special kpoints more than once. This is how band.constructpath works.
num_points_path = len(kPathPoints)
# initialize spectral function
spectral = np.zeros((NOmegas, num_points_path), dtype=float)
# before starting the calculation, force the chebyshev object to be complex, this is
# because the k-dependent Hamiltonian will be complex in any case
self._set_cheby(force_complex=True)
# stochastic calculation of moments
for dd in range(NSamplesDis):
qio.print_progress(verbose, "dis_head", dd, NSamplesDis)
# loop over selected indices
for k_idx, k in enumerate(kPathPoints):
# get bloch vectors corresponding to vector k
bloch = self.lb.get_bloch_states(k)
# the vector has to be repeated if spinful or if BCS
repeats = (1 + self.is_spin_matrix) * (1 + self.is_BCS)
# repeats is 1 if not spinful nor bcs, 2 or 4 if one or two are true. np.tile
# repeats the whole block of the vector, since we always first represent the up
# spins, then the down spins, or equivalent for creation/annihilation operators
# in Nambu space in the BCS case.
# NOTE: this is fine for now, but will depend on how the child Hamiltonian treats
# the spins! so not really general per-se.
bloch = np.tile(bloch, (repeats, 1))
# set Hamiltonian with boundary conditions that make the current k commensurate
self.set_Hamiltonian(kpt=k)
# initialize intermediate spectral function
tmp_spectral = np.zeros((NOmegas,), dtype=complex)
# loop over sites in the unit cell and spins (if necessary)
for imc in range(self.lb.M_uc * (1 + self.is_spin_matrix)):
# if spinful, the first M_uc will be the spin up, then the spin down. This way
# imc % M_uc gives us the site index.
if self.is_spin_matrix:
sel_ = "up" if imc // self.lb.M_uc == 0 else "down"
selector = self.generateSelector(sel_)
# selector is by default set to all sites, so it works in the spinless case
vec_ = bloch[:, imc % self.lb.M_uc] * selector
# NOTE: if not normalized again, the final spectral function will integrate
# to the number of particles (if the whole k space is spanned). This happens
# for all quantities at the moment, but should be consistenly changed once
# external feedback is given.
vec_ /= np.linalg.norm(vec_)
# calculate k-dependent green's function
_ = self.chebyshev.resolvent_greens_function(
energies,
1e-9,
Order + 1,
vec_,
vec_,
tmp_spectral,
resolvent_type="+",
advanced=False,
filter="jackson",
num_smooth=0,
)
# the final spectral function is then just -Im[GR_{kk}]
# accumulate result
spectral[:, k_idx] -= tmp_spectral.imag / np.pi
qio.print_progress(verbose - 1, "norm", k_idx, num_points_path)
# randomize the disorder configuration for next cycle
self.scrambleW()
qio.print_progress(verbose, "dis", dd, NSamplesDis)
# average calculated spectral function
spectral /= NSamplesDis
# reset disorder to original configuration
self.disLocPot = original_disLocPot.copy()
# restore user-defined boundary conditions
self.lb.BC = user_bc
# rebuild the original Hamiltonian
self.set_Hamiltonian()
# reset the Chebyshev object to the original one
self._set_cheby(force_complex=False)
return spectral, energies, locationSpecialPoints, indexSpecialPoints, kPath
[docs]
def getMu(
self,
T: float,
dos: Optional[np.ndarray] = None,
w: Optional[np.ndarray] = None,
NSamples: Optional[int] = None,
Order: Optional[int] = None,
NOmegas: Optional[int] = None,
NSamplesDis: int = 1,
verbose: int = 2,
) -> float:
r"""Calculate the chemical potential.
We calculate the chemical potential as the :math:`\mu` which satisfies
.. math::
N = \int\limits dE f(E-\mu) \rho(E)
where :math:`N` is the number of particles in the system, f is the fermi function and
:math:`\rho(E)` is the density of states.
In the case of BCS systems, the DOS depends on the chemical potential. This is therefore
calculated for every chemical potential and the corresponding number of particles is
obtained.
A bisection is used starting from the boundaries of the Hamiltonian and honing in on the
chemical potential resulting in a number of particles closest to the desired one. Starting
from the two values giving slightly more and slightly less particles than needed, a linear
extrapolation is finally performed to approximate the energy that results in the occupation
closest to the desired one.
TODO: it is probably better to split this calculation in the "one-shot" version and the
self-consistent one, suitable for BCS systems.
Args:
T (float): Temperature.
dos (Optional[np.ndarray]): precoumputed DoS. Will be the 1st one to be used in the
BCS case.
w (Optional[np.ndarray]): frequencies of the precomputed DoS.
NSamples (Optional[int]): Number of stochastic samples for the DoS calculation.
Order (Optional[int]): Order of the Chebyshev expansion for the DoS calculation.
NOmegas (Optional[int]): Number of frequencies for the DoS calculation.
NSamplesDis (int): Number of disorder configurations for the DoS calculation.
verbose (int): controls the verbosity of the method.
Raises:
ValueError: if the input is not consistent.
Returns:
float: Chemical potential.
"""
# check that input in consistent with system type (BCS or non BCS)
bcs_params = [NSamples, Order, NOmegas]
nonbcs_params = [dos, w]
# if BCS and not all bcs params are set or any nonbcs param is set, throw an error.
if self.is_BCS:
if any(x is None for x in bcs_params) or any(x is not None for x in nonbcs_params):
message = "The provided parameters are not compatible with the system. "
message += "For BCS systems, only provide the convergence parameters for the DoS "
message += "calculations needed for determining the chemical potential."
raise ValueError(message)
paramsDict = {"NSamples": NSamples, "Order": Order, "NOmegas": NOmegas}
# same applies for non bcs case
else:
if any(x is None for x in nonbcs_params) or any(x is not None for x in bcs_params):
message = "The provided parameters are not compatible with the system. "
message += "For non-BCS systems, only provide the DoS and corresponding "
message += "frequencies needed for determining the chemical potential."
raise ValueError(message)
# in principle I have already checked at this point that they are not none, but mypy
# does not understand this
if (dos is not None) and (w is not None) and (dos.shape != w.shape):
raise ValueError("The DOS and corresponing energies have incompatible shapes.")
paramsDict = {}
paramsDict["NSamplesDis"] = NSamplesDis
# at this point we know we have passed a correct set of arguments. This assignment and
# casting is only needed for mypy reasons.
_dos = cast(np.ndarray, dos)
_w = cast(np.ndarray, w)
_NSamples = cast(int, NSamples)
_Order = cast(int, Order)
_NOmegas = cast(int, NOmegas)
def _return_DOS(mu: float, seed: int) -> tuple[np.ndarray, np.ndarray]:
"""Function returning the DoS (recomputing it with the new chemical potential if BCS).
Args:
mu (float): Chemical potential.
seed (int): seed for the DoS calculation.
Returns:
tuple[np.ndarray, np.ndarray]: dos and corresponding frequencies.
"""
if self.is_BCS:
self.set_Hamiltonian(mu=mu)
# rescale Hamiltonian accordingly
self.scaleH_inv()
self.scaleH(1e-1)
# set seed and calculate DOS(mu)
np.random.seed(seed)
return self.getDOS(
NSamples=_NSamples,
Order=_Order,
NOmegas=_NOmegas,
NSamplesDis=NSamplesDis,
verbose=0,
)
return _dos, _w
qio.print_header(verbose, "chemical potential", paramsDict)
# FIXME: we will be able to remove this once fake_spinful will be removed
# for "fake spinful" systems, lb.N will refer to the spinful system, meaning a factor of
# two must be accounted for to get the correct number of particles. In all other cases
# realDim takes care of the normalization correctly.
_fake_spin_mult = 1
if (not self.is_spin_matrix) and (
self.lb.config["system"]["site_type"] == "spinful_fermions"
):
_fake_spin_mult = 2
# save shift and scale, we will reset them at the end
old_shift = self.scaling.get_shift()
old_scale = self.scaling.get_scale()
# determine the possible values for the bisection
if self.is_BCS:
# use the Hamiltonian boundaries as starting point in BCS systems
self.scaleH_inv()
mL, mR = self.getBoundsLanczos(100)
mus = np.linspace(mL, mR, _NOmegas)
else:
# use the frequencies in the non bcs case
mus = _w
_NOmegas = len(mus)
# indices of the extremes
iL = 0
iR = _NOmegas - 1
# random seed only needed for BCS dos calculations. Setting this to the same value every
# time we perform a DoS calculation ensures that the calculation is not affected by
# fluctuations in case of not enough stochastic sampling
seed = np.random.randint(2**32 - 1)
# this will continue until the left mu integrates to slightly less than N and the right one
# to slightly more (or equal)
finished = False
Niterations = np.ceil(np.log2(_NOmegas))
iiter = 0
while not finished:
# take middle index and calculate corresponding DOS and particle integral
iM = (iL + iR) // 2
# in BCS we need the DoS corresponding to the current chemical potential
_dos, _w = _return_DOS(mu=mus[iM], seed=seed)
# calculate particle number
dw = _w[1] - _w[0]
# the chemical potential for the N calculation should be 0 in the BCS case, since it
# is already 'included' in the Hamiltonian (hence the calculated DOS).
mu_for_fermi_integral = 0 if self.is_BCS else mus[iM]
fD = qt.fermi_distr(energies=_w, mu=mu_for_fermi_integral, T=T)
# calculate corresponding N
tmp_N = (_dos * fD).sum() * dw * self.realDim * _fake_spin_mult
# update extremes
iL = iL if tmp_N >= self.lb.N else iM
iR = iR if tmp_N < self.lb.N else iM
# stop if iR and iL are contiguous
if iR - iL == 1:
finished = True
qio.print_progress(verbose, "norm", iiter, Niterations)
iiter += 1
print("Finalizing calculation... ", end="", flush=True)
# once the two mus are found, a linear interpolation will be used to determine the mu that
# returns the number of particles closest to the desired one
partFin = {}
# this particle number has been calculated in the last iteration of the loop
partFin[iM] = tmp_N
# just find the other mu that has to be computed
iRemaining = iL if iL != iM else iR
# for BCS recompute the DOS with the new chemical potential
_dos, _w = _return_DOS(mu=mus[iRemaining], seed=seed)
# exactly as above
mu_for_fermi_integral = 0 if self.is_BCS else mus[iRemaining]
fD = qt.fermi_distr(energies=_w, mu=mu_for_fermi_integral, T=T)
# calculate corresponding N
partFin[iRemaining] = (_dos * fD).sum() * dw * self.realDim * _fake_spin_mult
# assume the number of particles to be a linear function between mus[iR] and mus[iL] and
# return the frequency that would best approximate self.lb.N with this assumption
mu_approx = mus[iL]
mu_approx += (self.lb.N - partFin[iL]) / (
(partFin[iR] - partFin[iL]) / (mus[iR] - mus[iL])
)
print("Done.", flush=True)
# reset Hamiltonian, shift and scale
self.set_Hamiltonian(mu=0)
self.scaling.set_scale(old_scale)
self.scaling.set_shift(old_shift)
return mu_approx
[docs]
def getCarrierDensity(self, dos: np.ndarray, w: np.ndarray, mu: float, T: float) -> float:
r"""Calculate the carrier density.
The calculation is performed by integrating the DoS time the occupation for frequencies
larger than the chemical potential and normalizing this by the system's volume. The formula
reads:
.. math::
n = \frac{1}{\Omega} \int\limits_{-\infty}^{\mu} f(E-\mu) \rho(E)
where :math:`\Omega` is the total volume, :math:`\mu` is the chemical potential, f is the
fermi function and :math:`\rho(E)` is the density of states.
Args:
dos (np.ndarray): Density of states.
w (np.ndarray): Energies for which the DOS is defined.
mu (float): Chemical potential.
T (float): Temperature.
Raises:
ValueError: If dos and w are not compatible.
Returns:
float: Carrier density in :math:`Å^{-D}` with D the dimensionality of the system.
"""
# TODO: This routine can become completely external. Just use the additional argument
# total_volume, since we multiply by realDim and then divide by vol_p_s * realDim
if dos.shape != w.shape:
raise ValueError("The DOS and corresponing energies have incompatible shapes.")
# can't use is_spin_matrix, see getMu
_spin_mult = 1 + int(self.lb.config["system"]["site_type"] == "spinful_fermions")
dw = w[1] - w[0]
# simply integrate dos * qt.fermi_distr for w > mu
wselector = w > mu
n_el = (dos[wselector] * qt.fermi_distr(w[wselector], mu, T)).sum() * dw * self.realDim
n_el *= _spin_mult
# n_el is the actual number in the system, not the density. We need to divide by
# the total volume to have a density
n_el /= self.total_volume
return n_el
[docs]
def getChargeSusceptibility(
self,
T: float,
mu: float,
dos: Optional[np.ndarray] = None,
w: Optional[np.ndarray] = None,
NSamples: Optional[int] = None,
Order: Optional[int] = None,
NOmegas: Optional[int] = None,
NSamplesDis: int = 1,
) -> float:
r"""Calculate the linear DC electric charge susceptibility.
This is equal to
.. math::
\chi_e = \frac{2 \pi}{\Omega} \frac{\partial N}{\partial \mu}
where :math:`\Omega` is the total volume, N is the number of particles and :math:`\mu` the
chemical potential.
The routine simply computes the numerical derivative of the number of particles with
respect to the chemical potential as an average of energies close to it.
TODO: add possibility to select sublattices, if deemed useful
Args:
T (float): Temperature.
mu (float): Chemical potential.
dos (Optional[np.ndarray]): precoumputed DoS.
w (Optional[np.ndarray]): frequencies of the precomputed DoS.
NSamples (Optional[int]): Number of stochastic samples for the DoS calculation.
Order (Optional[int]): Order of the Chebyshev expansion for the DoS calculation.
NOmegas (Optional[int]): Number of frequencies for the DoS calculation.
NSamplesDis (int): Number of disorder configurations for the DoS calculation.
Raises:
ValueError: if the input is not consistent.
Returns:
float: static charge susceptibility of the system.
"""
# check that input in consistent with system type (BCS or non BCS)
bcs_params = [NSamples, Order, NOmegas]
nonbcs_params = [dos, w]
# if BCS and not all bcs params are set or any nonbcs param is set, throw an error.
# TODO: this could go in the io module.
if self.is_BCS:
if any(x is None for x in bcs_params) or any(x is not None for x in nonbcs_params):
message = "The provided parameters are not compatible with the system. "
message += "For BCS systems, only provide the convergence parameters for the DoS "
message += "calculations needed for determining the charge susceptibility."
raise ValueError(message)
# same applies for non bcs case
else:
if any(x is None for x in nonbcs_params) or any(x is not None for x in bcs_params):
message = "The provided parameters are not compatible with the system. "
message += "For non-BCS systems, only provide the DoS and corresponding "
message += "frequencies needed for determining the charge susceptibility."
raise ValueError(message)
if (dos is not None) and (w is not None) and (dos.shape != w.shape):
raise ValueError("The DOS and corresponing energies have incompatible shapes.")
# at this point we know we have passed a correct set of arguments. This assignment and
# casting is only needed for mypy reasons.
_dos = cast(np.ndarray, dos)
_w = cast(np.ndarray, w)
_NSamples = cast(int, NSamples)
_Order = cast(int, Order)
_NOmegas = cast(int, NOmegas)
# to deal with fake spinful cases, see getMu
_fake_spin_mult = 1
if (not self.is_spin_matrix) and (
self.lb.config["system"]["site_type"] == "spinful_fermions"
):
_fake_spin_mult = 2
if self.is_BCS:
# use the Hamiltonian boundaries as starting point in BCS systems
mL, mR = self.scaling.get_energy_bounds()
mus = np.linspace(mL, mR, _NOmegas)
else:
# use the frequencies in the non bcs case
mus = _w
# find w that is the closest to the selected mu
imu = np.argmin(np.abs(mus - mu))
# take this and two adjacent values
mus = mus[[imu - 1, imu, imu + 1]]
# calculate the particle number for the selected mu and 2 adjacent values
N = np.zeros(3)
# use a seed for every dos calculation so that they do not depend on the fluctuations
# of the stochastic averaging
seed = np.random.randint(2**32 - 1)
for icount, mui in enumerate(mus):
if self.is_BCS:
self.set_Hamiltonian(mu=mui)
# setting the seed and calculate
np.random.seed(seed)
_dos, _w = self.getDOS(
NSamples=_NSamples,
Order=_Order,
NOmegas=_NOmegas,
NSamplesDis=NSamplesDis,
verbose=0,
)
dw = _w[1] - _w[0]
mu_for_fermi_integral = 0 if self.is_BCS else mui
fD = qt.fermi_distr(energies=_w, mu=mu_for_fermi_integral, T=T)
# N is just the integral of the dos times the fermi function
N[icount] = ((_dos * fD * dw)).sum() * self.realDim * _fake_spin_mult
dN = np.diff(N)
# the susceptibility is proportional to dN/dmu
chi_el = (dN / np.diff(mus)).mean()
# add prefactors
chi_el *= 2 * np.pi / (self.total_volume)
# reset Hamiltonian, shift and scale if it was changed (BCS case)
if self.is_BCS:
self.set_Hamiltonian(mu=0)
return chi_el
[docs]
def getMagneticSusceptibility(
self,
T: float,
mu: float,
dos_u: Optional[np.ndarray] = None,
dos_d: Optional[np.ndarray] = None,
w: Optional[np.ndarray] = None,
NSamples: Optional[int] = None,
Order: Optional[int] = None,
NOmegas: Optional[int] = None,
NSamplesDis: int = 1,
) -> float:
r"""Calculate the linear static magnetic susceptibility.
This is equal to
.. math::
\chi_m = \frac{2 \pi}{\Omega} \frac{\partial S_z}{\partial B_z}
where :math:`\Omega` is the total volume, :math:`S_z` is the spin operator in the z
direction and :math:`B_z` the magnetic field.
The routine simply computes the numerical derivative of the spin expectation value with
respect to the applied B field.
TODO: add possibility to select sublattices, if deemed useful
Args:
T (float): Temperature.
mu (float): Chemical potential.
dos_u (Optional[np.ndarray]): precoumputed DoS of spin up electrons.
dos_d (Optional[np.ndarray]): precoumputed DoS of spin down electrons.
w (Optional[np.ndarray]): frequencies of the precomputed DoS.
NSamples (Optional[int]): Number of stochastic samples for the DoS calculation.
Order (Optional[int]): Order of the Chebyshev expansion for the DoS calculation.
NOmegas (Optional[int]): Number of frequencies for the DoS calculation.
NSamplesDis (int): Number of disorder configurations for the DoS calculation.
Raises:
ValueError: if the input is not consistent.
Returns:
float: static magnetic susceptibility of the system.
"""
# The Bz field will induce (or strengthen) the energy splitting between spin up and spin
# down. The spin up energy will be shifted down by Bz/2 and, conversely, the spin down will
# be shifted up by the same amount. The DOS associated with either spin will change
# accordingly: shifting similarly for non-spinflip cases and shifting/mixing otherwise.
# This enters in the calculation of Sz since this is evaluated as (<N_u> - <N_d>) / 2.
# What is done here is to calculate <N_u/d> for the given chemical potential and for
# Bz=+-dw. This is equivalent to shifting the DOS's in the opposite direction (effect of
# Bz as mentioned above) in the non-spinflip case, while in the spinflip case to the same
# mu but with recomputed DoS(Bz).
sflip_params = [NSamples, Order, NOmegas]
nonsflip_params = [dos_u, dos_d, w]
# TODO: this could go in the io module.
if self.is_spinflip:
if any(x is None for x in sflip_params) or any(x is not None for x in nonsflip_params):
message = "The provided parameters are not compatible with the system. "
message += "For spinflip systems, only provide the convergence parameters for the "
message += "DoS calculations needed for determining the magnetic susceptibility."
raise ValueError(message)
# same applies for non-spinflip case
else:
if any(x is None for x in nonsflip_params) or any(x is not None for x in sflip_params):
message = "The provided parameters are not compatible with the system. For "
message += "non-spinflip systems, only provide the spinful DoS and corresponding "
message += "frequencies needed for determining the charge susceptibility."
raise ValueError(message)
if (
(dos_u is not None)
and (dos_d is not None)
and (w is not None)
and ((dos_u.shape != w.shape) or (dos_u.shape != dos_d.shape))
):
raise ValueError("The DOS and/or corresponing energies have incompatible shapes.")
# at this point we know we have passed a correct set of arguments. This assignment and
# casting is only needed for mypy reasons.
_dos_u = cast(np.ndarray, dos_u)
_dos_d = cast(np.ndarray, dos_d)
_w = cast(np.ndarray, w)
_NSamples = cast(int, NSamples)
_Order = cast(int, Order)
_NOmegas = cast(int, NOmegas)
# to deal with fake spinful cases, see getMu
_fake_spin_mult = 1
if (not self.is_spin_matrix) and (
self.lb.config["system"]["site_type"] == "spinful_fermions"
):
_fake_spin_mult = 2
if self.is_spinflip:
# use the Hamiltonian boundaries as starting point in BCS systems
self.scaleH_inv()
mL, mR = self.getBoundsLanczos(100)
_w = np.linspace(mL, mR, _NOmegas)
self.scaleH(0.2)
def _return_DOSs(Bz: float, seed: int) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Function returning the DoS up and down.
The DoS#s are recomputed with the new Bz if the system is spinflip.
Args:
Bz (float): Magnetic field in the z direction.
seed (int): seed for the DoS calculation.
Returns:
tuple[np.ndarray, np.ndarray, np.ndarray]: dos and corresponding frequencies.
"""
# TODO: this will have to be the spinflip flag
if self.is_spinflip:
self.set_Hamiltonian(Bz=Bz)
# rescale Hamiltonian accordingly
self.scaleH_inv()
self.scaleH(1e-1)
# set seed and calculate DOS(mu)
np.random.seed(seed)
recalc_dos_u, recalc_w = self.getDOS(
NSamples=_NSamples,
Order=_Order,
NOmegas=_NOmegas,
NSamplesDis=NSamplesDis,
selector=self.generateSelector("up"),
verbose=0,
)
np.random.seed(seed)
recalc_dos_d, recalc_w = self.getDOS(
NSamples=_NSamples,
Order=_Order,
NOmegas=_NOmegas,
NSamplesDis=NSamplesDis,
selector=self.generateSelector("down"),
verbose=0,
)
return recalc_dos_u, recalc_dos_d, recalc_w
return _dos_u, _dos_u, _w
# save shift and scale, we will reset them at the end
old_shift = self.scaling.get_shift()
old_scale = self.scaling.get_scale()
# find w that is the closest to the selected mu
imu = np.argmin(np.abs(_w - mu))
# I will perform the calulation for Bz - dBz, Bz, and Bz + dBz, with dBz = 2dw, where dw is
# the step for the frequencies passed (>0). This translates in the following calculations:
# Bz - dBz --> N_u with mu->mu - dBz/2 since the corresponding dos would be shifted UP by
# the smaller Bz field. The opposite happens for the spin down.
# Bz --> just the GS N_u/d calculation with mu->mu
# Bz + dBz --> the opposite of the first case, dos_up is shifted down, meaning equivalent
# higher mu, mu->mu + dBz/2
# NOTE: Bz is already in the Hamiltonian, there's no need to explicitly add its effect here
# initialize N, first column will be spin down, the second spin up.
N = np.zeros((3, 2))
seed = np.random.randint(2**32 - 1)
# NOTE: the reason for sweeping the values and not jsut do mu +- dmu is that 1) I allow
# for non uniform spacing, 2) even w generated via np.linspace are not perfectly uniform
# due to floating point approximations and can give wrong results for T=0 where the
# fermi function is a step function
list_mus = _w[[imu - 1, imu, imu + 1]]
for icount, iimu in enumerate(list_mus):
Bz = iimu - _w[imu]
# TODO: please rewrite this madness with iimu and imu and so on
_dos_u, _dos_d, _w = _return_DOSs(Bz=Bz, seed=seed)
mu_fermi_int_u = _w[imu]
mu_fermi_int_d = _w[imu]
# for a certain mu, e.g. mu = mu - dBz / 2, the corresponding calculation for the
# number of particles will be the one corresponding to Bz = Bz - dBz for spin up and
# Bz = Bz + dBz for spin down.
if not self.is_spinflip:
# this convoluted way is needed because of errors introduced in the calculation
# of the integrals because of roundoff errors when mu = mu +- Bz / 2 not aligning
# perfectly with the saved w, making one single contribution miss in the integral.
tmp = np.argmin(np.abs(_w - (_w[imu] + Bz / 2)))
mu_fermi_int_u = _w[tmp]
tmp = np.argmin(np.abs(_w - (_w[imu] - Bz / 2)))
mu_fermi_int_d = _w[tmp]
# N is just the integral of the dos times the fermi function with the proper mu
# NOTE: if the frequencies are not uniform, this will introduce an error depending on
# the difference between the steps in frequencies on either side.
dws = np.diff(_w, prepend=(2 * _w[0] - _w[1]))
N[icount, 1] = (_dos_u * qt.fermi_distr(_w, mu_fermi_int_u, T) * dws).sum()
N[icount, 1] *= self.realDim * _fake_spin_mult
N[icount, 0] = (_dos_d * qt.fermi_distr(_w, mu_fermi_int_d, T) * dws).sum()
N[icount, 0] *= self.realDim * _fake_spin_mult
# calculate Sz and dSz
Sz = np.diff(N, axis=1).squeeze() / 2
dSz = np.diff(Sz)
chi_m = (dSz / np.diff(list_mus)).mean()
# add prefactors
chi_m *= 2 * np.pi / (self.total_volume)
# reset Hamiltonian, shift and scale
self.set_Hamiltonian(Bz=0)
self.scaling.set_scale(old_scale)
self.scaling.set_shift(old_shift)
return chi_m
[docs]
def getLinearResponse(
self,
A: Callable,
P: 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 static linear response of an observable to a perturbation.
This is done by expanding
.. math::
\chi^A_P = \frac{2 \pi}{\Omega} \frac{\partial \langle \hat{A} \rangle}{\partial
\varepsilon}
where :math:`\hat{A}` is the observable, :math:`\hat{P}` the perturbation and
:math:`\varepsilon` its strength, i.e. the perturbation enters the Hamiltonian as
:math:`\hat{H}' = \varepsilon \hat{P}` with :math:`\varepsilon << 1`.
Args:
A (Callable): Linear operator associated with the observable.
P (Callable): Perturbation in linear operator form.
dtype (type): dtype of the linear operators.
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,
"linear response of provided observable",
{
"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()
# save original disorder configuration to be reset at the end
original_disLocPot = self.disLocPot.copy()
# use complex if either the Hamiltonian or the Operator are complex
is_anything_complex = complex in [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, same for random vector
mu_mn: np.ndarray = np.zeros((Order + 1, Order + 1), dtype=moments_dtype)
tmp_mu_mn: np.ndarray = np.zeros((Order + 1, Order + 1), dtype=moments_dtype)
x0 = np.zeros(self.dim, dtype=self.dtype)
# stochastic loop
# 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
x0 = qt.generateRandomVector(selector, x0, moments_dtype)
xR = np.zeros_like(x0)
xR = P(x0, out=xR)
# calculate and accumulate Chebyshev moments. According to the notes, we have
# <r| A Tn P Tm |r>. By using the cyclic property of the trace, we obtain
# <r| Tn P Tm A|r>, hence the calculation of xR above.
self.chebyshev.double_expansion_fast(
n_cheby_moments=Order + 1, linear_operator_V=A, bra=x0, ket=xR, res=tmp_mu_mn
)
mu_mn += tmp_mu_mn.astype(moments_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
mu_mn = mu_mn / (NSamples * NSamplesDis) * self.realDim
# reset disorder to original configuration and reset the cheby object
self.disLocPot = original_disLocPot.copy()
self._set_cheby(force_complex=False)
# get Fermi-Pade frequencies and residues to perform the matsubara summation
# TODO: this will have to be reworked once we have a tools package. Also KB should not
# be defined here.
izf, Rf = qt.getFermiPadeParams(NPadeFreq)
kB = 8.617333262145e-5
beta = 1 / (kB * T)
izf = izf / beta + mu
Rf /= beta
# get the Green's function coefficients
gm = calc_prefactors(self.scaling, izf, Order + 1)
# normalization coefficient based on the percentage of selected states.
perc_selected: float = selector.sum() / self.realDim
# evaluate the sum over the moments and over the pade frequencies
tmp: np.ndarray = np.einsum("zn, mn -> zm", gm, mu_mn)
tmp = np.einsum("zm, zm -> z", gm, tmp)
# normalization
lin_resp = (self.scaling.get_scale() ** 2) * tmp @ Rf
lin_resp *= perc_selected
lin_resp *= 2 * np.pi / self.total_volume
lin_resp = lin_resp.real if dtype is complex else lin_resp
# cast the value to float in the corresponding case
return lin_resp
[docs]
def getOperator(self, operator: str) -> tuple[Callable, type]:
"""Provide basic operators.
Operators available:
- ±number operator and its left version
- ±spin operator and its left version
- ±position and its left version
Args:
operator (str): operator to be returned.
Raises:
NotImplementedError: if the requested operator has not been implemented.
Returns:
tuple[Callable, type]: tuple contining the operator as a callable acting on an input
array and the dtype of the operator.
"""
sign = 2 * int(operator[0] != "-") - 1
operator = operator.strip("-")
# declaration of dtype
_dtype: type
if operator in ["N", "NDag"]:
_dtype = float
# define operator
def _operator(_x: np.ndarray, out: np.ndarray) -> np.ndarray:
"""Number operator.
Args:
_x (np.ndarray): input vector.
out (np.ndarray): output vector.
Returns:
np.ndarray: Number operator applied to the input vector.
"""
# simply +1 for particles and -1 for holes
out[: self.realDim] = _x[: self.realDim]
out[self.realDim :] = -_x[self.realDim :]
return sign * out
elif operator in ["Sz", "SzDag"]:
_dtype = float
# define operator
def _operator(_x: np.ndarray, out: np.ndarray) -> np.ndarray:
"""Spin operator.
Args:
_x (np.ndarray): input vector.
out (np.ndarray): output vector.
Returns:
np.ndarray: Spin operator applied to the input vector.
"""
# simply 1/2 for spin up (1st & 3rd quarters (in BCS)) and -1/2 otherwise
quarter = self.realDim // 2
out[:quarter] = _x[:quarter] / 2
out[quarter : self.realDim] = -_x[quarter : self.realDim] / 2
out[self.realDim : 3 * quarter] = _x[self.realDim : 3 * quarter] / 2
out[3 * quarter :] = -_x[3 * quarter :] / 2
return sign * out
elif operator in ["rx", "rxDag", "ry", "ryDag", "rz", "rzDag"]:
_dtype = float
idx_dir, _ = qt.dir2idx(operator[1] + operator[1])
r = self.lb.get_atom_positions()[:, idx_dir]
# define operator
def _operator(_x: np.ndarray, out: np.ndarray) -> np.ndarray:
"""Position operator.
Args:
_x (np.ndarray): input vector.
out (np.ndarray): output vector.
Returns:
np.ndarray: Position operator applied to the input vector.
"""
# simply r * _x
out[:] = r * _x
return sign * out
else:
raise NotImplementedError(f"Operator {operator} not implemented.")
return (_operator, _dtype)