"""Quantarhei package (http://www.github.com/quantarhei)
systembathinteraction module
"""
from __future__ import annotations
from typing import TYPE_CHECKING, Any
import numpy
from ... import REAL
from ...core.dfunction import DFunction
from ...core.saveable import Saveable
from ...core.time import TimeAxis
from ...exceptions import QuantarheiError
from ...qm.corfunctions.cfmatrix import CorrelationFunctionMatrix
from ...qm.corfunctions.correlationfunctions import c2g
from ...qm.corfunctions.functionstorage import FunctionStorage
if TYPE_CHECKING:
from ...builders.aggregates import Aggregate
from ...builders.molecules import Molecule
from ...builders.opensystem import OpenSystem
[docs]
class SystemBathInteraction(Saveable):
"""Describes interaction of an open quantum system with its environment
Stores the system--bath interaction operator in form of a set of operators
on the Hilbert space of the system and correlation functions of
the operator on the bath Hilbert space,
OR
It stores various relaxation and dephasing rates, to represent e.g. the
Linblad form.
Parameters
----------
sys_operators : list
List of the system part of the system-bath interaction Hamiltonian
components
bath_correlation_matrix: CorrelationFunctionMatrix
Object of the CorrelationFunctionMatrix type holding all correlation
functions needed for the description of system bath interaction
rates : list/tuple
List or tuple of rates. The total number of rates has to be
the same as the number of operators in the ``sys_operator`` list
drates : array
An array of dephasing rates. The dimension of the array must
correspond to the dimension of the treated system.
dtype : str
Type of the dehasing defined in `drates`. The types are "Lorentzian"
which is default, and correponds to a dephasing rate equation with
the term -\\gamma\rho_{ab} on the right hand side. The dephasing
is exponential. The type "Gaussian" results in a Gaussian dephasing
and corresponds to the term -\\gamma t \rho_{ab} on the right hand
side of the rate equation.
osites : list, array
List or array of site indices on which the oscillators reside. The
indices can repeat, indicating several modes on a single site.
orates : array
Oscillator decay rates. This rates corresponds to the dephasing
rate of the oscillations in a harmonic oscillator
system : {Molecule, Aggregate}
Molecule or Aggregate object in which the system--bath interaction
is specified
"""
def __init__(
self,
sys_operators: Any = None,
bath_correlation_matrix: Any = None,
rates: Any = None,
drates: Any = None,
dtype: str = "Lorentzian",
osites: Any = None,
orates: Any = None,
system: Any = None,
) -> None:
# information about aggregate is needed when dealing with
# multiple excitons
self.aggregate: Aggregate | None = None
self.molecule: Molecule | None = None
self.system: Aggregate | Molecule | OpenSystem | None = None
self.rates: Any = None
self.KK: numpy.ndarray | None = None
self.CC: CorrelationFunctionMatrix | None = None # correlation function matrix
self.GG: FunctionStorage | None = None # lineshape function storage
self.TimeAxis: TimeAxis | None = None
self.drates: Any = None
self.N = 0
self.osites: Any = None
self.orates: Any = None
self.sbitype = "Linear_Coupling"
self._has_gg_storage = False
self._gg_storage_config: dict | int | None = None
#
# version with bath correlation functions
#
if (sys_operators is not None) and (bath_correlation_matrix is not None):
self.sbitype = "Linear_Coupling"
# Find the length of the list of operators
if isinstance(sys_operators, list):
self.N = len(sys_operators)
else:
raise QuantarheiError("sys_operators argument has to a list")
# Second argument has to be a CorrelationFunctionMatrix
if not isinstance(bath_correlation_matrix, CorrelationFunctionMatrix):
raise QuantarheiError(
"ba_correlation_function argument has to a"
" CorrelationFunctionMatrix"
)
# Check that sys_operators and bath_correlation matrix has
# a compatible number of components
if bath_correlation_matrix.nob != self.N:
raise QuantarheiError(
f"Incompatile number of bath compoments: "
f"Correlation function matrix - {bath_correlation_matrix.nob:d} vs. operators {self.N:d}"
)
self.TimeAxis = bath_correlation_matrix.timeAxis
self.set_system(system)
if self.N > 0:
self._set_operators(sys_operators)
self.CC = bath_correlation_matrix
if rates is not None:
if len(rates) != self.N:
raise QuantarheiError("Wrong number of rates specified")
self.rates = rates
#
# version with system-bath operators and rates Lindblad form
#
elif (sys_operators is not None) and (rates is not None):
self.sbitype = "Lindblad_Form"
# Find the length of the list of operators
if isinstance(sys_operators, list):
self.N = len(sys_operators)
else:
raise QuantarheiError("First argument has to a list")
self.set_system(system)
if self.N > 0:
self._set_operators(sys_operators)
self.CC = None # bath_correlation_matrix
if len(rates) != self.N:
raise QuantarheiError("Wrong number of rates specified")
self.rates = rates
#
# version with phenomenological pure dephasing
#
elif (sys_operators is None) and (drates is not None):
self.sbitype = "Pure_Dephasing"
if len(drates.shape) != 2:
raise QuantarheiError(
"Pure dephasing rates must be defined by a matrix"
)
if drates.shape[0] != drates.shape[1]:
raise QuantarheiError(
"Pure dephasing rates must be defined by a square matrix"
)
self.N = drates.shape[0]
self.set_system(system)
self.CC = None
#
# version with Lindblad form for vibrational modes
#
elif (sys_operators is None) and (orates is not None):
self.sbitype = "Vibrational_Lindblad_Form"
if len(orates) != len(osites):
raise QuantarheiError(
"`orates` and `osites` arguments must have the same lengths"
)
self.set_system(system)
self.CC = None
self.orates = orates
self.osites = osites
[docs]
def set_system(self, system: Any) -> None:
"""Sets the system attribute"""
from ...builders.aggregates import Aggregate
from ...builders.molecules import Molecule
from ...builders.opensystem import OpenSystem
if system is not None:
if isinstance(system, Aggregate):
self.aggregate = system
self.molecule = None
self.system = self.aggregate
elif isinstance(system, Molecule):
self.aggregate = None
self.molecule = system
self.system = self.molecule
elif isinstance(system, OpenSystem):
self.aggregate = None
self.molecule = None
self.system = system
else:
raise QuantarheiError("Unknown system type")
def _set_operators(self, sys_operators: Any) -> None:
"""Sets the system part of the interaction"""
# First of the operators
KK = sys_operators[0]
# Get its dimension
dim = KK.data.shape[0]
self.KK = numpy.zeros((self.N, dim, dim), dtype=REAL)
self.KK[0, :, :] = numpy.real(KK.data)
# Save other operators and check their dimensions
for ii in range(1, self.N):
KK = sys_operators[ii]
if dim == KK.data.shape[0]:
self.KK[ii, :, :] = numpy.real(KK.data)
else:
raise QuantarheiError(
"Operators in the list are not of the same dimension"
)
[docs]
def get_time_axis(self) -> Any:
"""Returns the time axis of the storred correlation functions"""
return self.TimeAxis
[docs]
def get_correlation_function(self, where: Any) -> Any:
"""Returns the bath correlation function object defined by a pair of sites (tuple)"""
return self.CC.get_correlation_function(where[0], where[1])
[docs]
def get_coft(self, n: int, m: int) -> Any:
"""Returns bath correlation function corresponding to sites n and m"""
if self.sbitype != "Linear_Coupling":
raise QuantarheiError(
"Correlation functions only defined for "
"linear microscopic system-bath coupling"
)
# print("Returning coft" )
if self.system is None:
# print("Returning coft without the system" )
return self.CC.get_coft(n, m)
# FIXME: Molecule needs this method
bn = self.system.which_band[n]
bm = self.system.which_band[m]
# Ground state
if (bn == 0) and (bm == 0):
#
# This returns zero correlation function, which is consistent
#
assert self.CC is not None
return self.CC._cofts[0, :]
# First or higher excited state bands
if (bn >= 1) and (bm >= 1):
# print(bn,bm,"::",n-1,m-1)
#
# First band starts with n=1, but the correlation functions
# are stored by site index which starts with 0
#
return self.CC.get_coft(n - 1, m - 1)
# other bands return zero correlation function now
assert self.CC is not None
return self.CC._cofts[0, :]
[docs]
def get_coft_elsig(self, n_sig: Any, m_sig: Any) -> Any:
"""Returns bath correlation based on electronic signatures"""
if self.sbitype != "Linear_Coupling":
raise QuantarheiError(
"Correlation functions only defined for "
"linear microscopic system-bath coupling"
)
nb = numpy.sum(n_sig)
mb = numpy.sum(m_sig)
indices = []
if mb == nb:
ni = 0
for na in n_sig:
mi = 0
for ma in m_sig:
if (na == 1) and (ma == 1):
indices.append([ni, mi])
mi += 1
ni += 1
ret = numpy.zeros((self.TimeAxis.length), dtype=numpy.complex128)
for ind in indices:
# print(nb,":",ind[0],ind[1])
ret += self.get_coft(ind[0], ind[1])
return ret
assert self.CC is not None
return self.CC._cofts[0, :]
[docs]
def get_goft_storage(
self, config: dict | int | None = None, timeaxis: Any = None
) -> Any:
"""Returns a lineshape function storage based on correlation functions
The function calculates g(t) fuctions based on the correlation
functions specified in this object. This call fails if correlation
functions are not specified.
Parameters
----------
config : dict
Dictionary of the FunctionStorage configuration. If None submitted
it wil produced g(t) for a standard 3rd order response calculation.
"""
using_default_timeaxis = timeaxis is None
if (
using_default_timeaxis
and self._has_gg_storage
and (config is None or config == self._gg_storage_config)
):
return self.GG
# number of functions
Nf = self.CC.nof
# number of sites
Nb = self.CC.nob
# we define storage for lineshape functions with prescribed config
# FIXME: orinally I had Nb here, but that is wrong. Nf also does not work
if using_default_timeaxis:
timeaxis = [self.TimeAxis, self.TimeAxis]
gg = FunctionStorage(Nf, timeaxis=timeaxis, config=config)
# create functions to update storage
fcions = {}
for kk in range(Nb):
ii = self.CC.get_index_by_where((kk, kk))
if ii not in fcions:
# correlation function in discrete representation
cf = self.get_coft(kk + 1, kk + 1)
# integration into g(t)
gf = c2g(self.TimeAxis, cf)
# make it into spline function
df = DFunction(self.TimeAxis, gf)
gfunc = df.as_spline_function()
# save for later
fcions[ii] = gfunc
# set a function for a give site
gg.set_goft(kk, func=fcions[ii])
# make checks
if Nf != gg.get_number_of_functions():
raise QuantarheiError(
"Number of functions did not conserve in g(t) calculation"
)
if using_default_timeaxis:
self.GG = gg
self._has_gg_storage = True
self._gg_storage_config = config
return gg
[docs]
def has_temperature(self) -> bool:
"""Checkst if the Aggregate has a defined temperature"""
if (
self.sbitype == "Lindblad_Form"
or self.sbitype == "Vibrational_Lindblad_Form"
):
return False
try:
T = self.get_temperature()
if T >= 0.0:
return True
return False
except Exception:
return False
[docs]
def get_temperature(self) -> float:
"""Returns temperature associated with the bath"""
return self.CC.get_temperature()
[docs]
def get_reorganization_energy(self, i: int, j: int | None = None) -> Any:
"""Returns reorganization energy associated with a given site
If one index `i` is specified, the function returns reorganization
energy associated with i-th site. If two indices are specified,
the function returns crosscorrelation function (when i not equal j)
or a correlation function of the site (when i = j)
"""
if (
self.sbitype == "Lindblad_Form"
or self.sbitype == "Vibrational_Lindblad_Form"
):
return None
if j is None:
j = i
return self.CC.get_reorganization_energy(i, j)
[docs]
def get_correlation_time(self, i: int, j: int | None = None) -> Any:
if (
self.sbitype == "Lindblad_Form"
or self.sbitype == "Vibrational_Lindblad_Form"
):
return None
if j is None:
j = i
return self.CC.get_correlation_time(i, j)
[docs]
def get_sbitype(self) -> str:
"""Returns the type of SystemBathInteraction
Options are `Lindblad_Form`, `Vibrational_Lindblad_Form`
and `Linear_Coupling`
"""
return self.sbitype