"""Provides typical Bath correlation function types.
Most important types of bath or energy gap correlation functions are
provided. Where possible, the correlation function is calculated
from the parameters from analytical formulae. Where such formulae are
not available, correlation function is calculated by transformation
of the spectral density.
Types of correlation function provided
--------------------------------------
OverdampedBrownian-HighTemperature :
OverdampedBrownian oscillator in high temperature limit
OverdampedBrownian :
General overdampedBrownian oscillator
Examples
--------
>>> from quantarhei import TimeAxis
>>> params = dict(ftype="OverdampedBrownian", cortime=100, reorg=20, T=300)
>>> time = TimeAxis(0.0,1000,1.0)
>>> with energy_units("1/cm"):
... cf = CorrelationFunction(time,params)
>>> with energy_units("1/cm"):
... print(cf.get_reorganization_energy())
20.0
Reorganization energy of a correlation function can be calculated from the
shape of the spectral density by integrating over it. The accuracy
of such estimation depends on numerics, hence the relative tolerance of
only 1.0e-4 below
>>> lamb_definition = cf.get_reorganization_energy()
>>> lamb_measured = cf.measure_reorganization_energy()
>>> print(numpy.allclose(lamb_definition, lamb_measured, rtol=1.0e-4))
True
Details of Classes Provided
---------------------------
"""
from __future__ import annotations
from copy import deepcopy
from typing import Any, cast
import numpy
import numpy.typing
import scipy.interpolate as interp
from ...core.dfunction import DFunction
from ...core.frequency import FrequencyAxis
from ...core.managers import UnitsManaged, energy_units
from ...core.time import TimeAxis
from ...core.units import kB_intK
from ...core.wrappers import enforce_energy_units_context
from ...exceptions import QuantarheiError
[docs]
class CorrelationFunction(DFunction, UnitsManaged):
"""Bath (energy-gap) correlation function.
Provides several standard spectral-density models as analytical or
numerically evaluated correlation functions. The type is selected via
the ``'ftype'`` key of ``params``.
Parameters
----------
axis : TimeAxis
Time grid on which the correlation function is evaluated.
params : dict or list of dict
Parameter dictionary (or list of dictionaries for a composite
correlation function). Must contain at least ``'ftype'``, ``'reorg'``,
and ``'T'`` (temperature in Kelvin).
values : numpy.ndarray, optional
Pre-computed complex values of the correlation function at the
time points of ``axis``. When supplied, the analytical evaluation
is skipped.
"""
allowed_types = (
"OverdampedBrownian-HighTemperature",
"OverdampedBrownian",
"OverdampedBrownian_from_Specdens",
"UnderdampedBrownian",
"Underdamped",
"B777",
"CP29",
"Value-defined",
)
analytical_types = ("OverdampedBrownian-HighTemperature", "OverdampedBrownian")
energy_params = (
"reorg",
"omega",
"freq",
"fcp",
"g_FWHM",
"l_FWHM",
"freq1",
"freq2",
"gamma",
)
energy_units: str
lamb: float | numpy.ndarray
temperature: float | numpy.ndarray
axis: TimeAxis
@enforce_energy_units_context
def __init__(
self, axis: Any = None, params: Any = None, values: Any = None
) -> None:
super().__init__()
if (axis is not None) and (params is not None):
# FIXME: values might also need handling according to specified
# energy units
self.energy_units = self.manager.get_current_units("energy")
# handle axis (it can be TimeAxis or FrequencyAxis)
if not isinstance(axis, TimeAxis):
taxis = axis.get_TimeAxis()
self.axis = taxis
else:
self.axis = axis
# handle params
self.params = [] # this will always be a list of components
p2calc = []
try:
# if this passes, we assume params is a dictionary
params.keys()
self._is_composed = False
p2calc.append(params)
except AttributeError:
# othewise we assume it is a list of dictionaries
self._is_composed = True
for p in params:
p2calc.append(p)
self.lamb = 0.0
self.temperature = -1.0
self.cutoff_time = 0.0
for params in p2calc:
try:
ftype = params["ftype"]
if ftype not in CorrelationFunction.allowed_types:
raise QuantarheiError("Unknown CorrelationFunction type")
# we mutate the parameters into internal units
prms = {}
for key in params.keys():
if key in self.energy_params:
prms[key] = self.convert_energy_2_internal_u(params[key])
else:
prms[key] = params[key]
except KeyError:
raise QuantarheiError(
"Dictionary of parameters does not contain `ftype` key"
)
self.params.append(prms)
if values is None:
#
# loop over parameter sets
#
for prms in self.params:
# try:
# ftype = params["ftype"]
#
# if ftype not in CorrelationFunction.allowed_types:
# raise QuantarheiError("Unknown CorrelationFunction type")
#
# # we mutate the parameters into internal units
# prms = {}
# for key in params.keys():
# if key in self.energy_params:
# prms[key] = self.convert_energy_2_internal_u(params[key])
# else:
# prms[key] = params[key]
#
# except:
# raise QuantarheiError("Dictionary of parameters does not contain "
# +" `ftype` key")
ftype = prms["ftype"]
if ftype == "OverdampedBrownian-HighTemperature":
self._make_overdamped_brownian_ht(prms) # , values=values)
elif ftype == "OverdampedBrownian":
self._make_overdamped_brownian(prms) # , values=values)
elif ftype == "UnderdampedBrownian":
self._make_underdamped_brownian(prms) # , values=values)
elif ftype == "Underdamped":
self._make_underdamped(params, values=values)
elif ftype == "B777":
self._make_B777(params, values=values)
elif ftype == "CP29":
self._make_CP29_spectral_density(params, values=values)
elif ftype == "Value-defined":
self._make_value_defined(prms, values)
else:
raise QuantarheiError(
"Unknown correlation function type or"
"type domain combination."
)
# self.params.append(prms)
else:
self._add_me(self.axis, values)
# update reorganization energy
self.lamb = 0.0
self.temperature = self.params[0]["T"]
for prms in self.params:
self.lamb += prms["reorg"]
if self.temperature != prms["T"]:
raise QuantarheiError(
"Inconsistent temperature! "
"Temperatures of all "
"components have to be the same"
)
# FIXME: set cut-off time and temperature
# self._set_temperature_and_cutoff_time(self.params[0])
def __repr__(self) -> str:
ftype = self.params[0]["ftype"] if self.params else "unknown"
temp = self.temperature if hasattr(self, "temperature") else None
return f"CorrelationFunction(ftype={ftype!r}, T={temp})"
def _matsubara(self, kBT: float, ctime: float, nof: int) -> numpy.ndarray:
"""Matsubara frequency part of the Brownian correlation function"""
msf: float | numpy.ndarray = 0.0
nut = 2.0 * numpy.pi * kBT
time = self.axis.data
for i in range(0, nof):
n = i + 1
msf += (
nut
* n
* numpy.exp(-nut * n * time)
/ ((nut * n) ** 2 - (1.0 / ctime) ** 2)
)
return numpy.asarray(msf)
def _set_temperature_and_cutoff_time(
self, temperature: float, ctime: float
) -> None:
"""Sets the temperature and cutoff time of for the component"""
# Temperatures of all components have to be the same
# is this the first time that temperature is assigned?
if self.temperature == -1.0:
self.temperature = temperature
elif self.temperature != temperature:
raise QuantarheiError(
"Inconsistent temperature! Temperatures of all "
"components have to be the same"
)
# longest cortime has to be preserved
new_cutoff_time = 5.0 * ctime
if new_cutoff_time > self.cutoff_time:
self.cutoff_time = new_cutoff_time
def _make_overdamped_brownian(self, params: dict) -> None: # , values=None):
"""Creates the overdamped Brownian oscillator component
of the correlation function
"""
temperature = params["T"]
ctime = params["cortime"]
lamb = params["reorg"]
if "matsubara" in params.keys():
nmatsu = params["matsubara"]
else:
nmatsu = 10
kBT = kB_intK * temperature
time = self.axis.data
cfce = (lamb / ctime) * (
1.0 / numpy.tan(1.0 / (2.0 * kBT * ctime))
) * numpy.exp(-time / ctime) - 1.0j * (lamb / ctime) * numpy.exp(-time / ctime)
cfce += (4.0 * lamb * kBT / ctime) * self._matsubara(kBT, ctime, nmatsu)
# this is a call to the function inherited from DFunction class
self._add_me(self.axis, cfce)
# update reorganization energy
self.lamb += lamb
# check temperature and update cutoff time
self._set_temperature_and_cutoff_time(temperature, 5.0 * ctime)
def _make_overdamped_brownian_ht(self, params: dict) -> None: # , values=None):
"""Creates the high temperature overdamped Brownian oscillator
component of the correlation function
"""
temperature = params["T"]
ctime = params["cortime"]
lamb = params["reorg"]
kBT = kB_intK * temperature
time = self.axis.data
cfce = (2.0 * lamb * kBT) * numpy.exp(-time / ctime) - 1.0j * (
lamb / ctime
) * numpy.exp(-time / ctime)
# this is a call to the function inherited from DFunction class
self._add_me(self.axis, cfce)
# update reorganization energy
self.lamb += lamb
# check temperature and update cutoff time
self._set_temperature_and_cutoff_time(temperature, 5.0 * ctime)
def _make_underdamped_brownian(self, params: dict) -> None: # , values=None):
"""Creates underdamped Brownian oscillator component of the correlation
function
"""
from .spectraldensities import SpectralDensity
temperature = params["T"]
ctime = params["gamma"]
# omega = params["freq"]
lamb = params["reorg"]
# kBT = kB_intK*temperature
time = self.axis # .data
with energy_units("int"):
# Make it via SpectralDensity
fa = SpectralDensity(time, params)
cf = fa.get_CorrelationFunction(temperature=temperature)
cfce = cf.data
# this is a call to the function inherited from DFunction class
self._add_me(self.axis, cfce)
# update reorganization energy
self.lamb += lamb
# check temperature and update cutoff time
self._set_temperature_and_cutoff_time(temperature, 5.0 / ctime)
def _make_underdamped(self, params: dict, values: Any = None) -> None:
from .spectraldensities import SpectralDensity
temperature = params["T"]
ctime = params["gamma"]
# use the units in which params was defined
lamb = params["reorg"]
time = self.axis # .data
if values is not None:
cfce = values
else:
# Make it via SpectralDensity
fa = SpectralDensity(time, params)
cf = fa.get_CorrelationFunction(temperature=temperature)
cfce = cf.data
# this is a call to the function inherited from DFunction class
self._add_me(self.axis, cfce)
# update reorganization energy
self.lamb += lamb
# check temperature and update cutoff time
self._set_temperature_and_cutoff_time(temperature, 5.0 * ctime)
def _make_B777(self, params: dict, values: Any = None) -> None:
from .spectraldensities import SpectralDensity
temperature = params["T"]
ctime = params["gamma"]
# use the units in which params was defined
lamb = self.manager.iu_energy(params["reorg"], units=self.energy_units)
time = self.axis # .data
if values is not None:
cfce = values
else:
# Make it via SpectralDensity
fa = SpectralDensity(time, params)
cf = fa.get_CorrelationFunction(temperature=temperature)
cfce = cf.data
# this is a call to the function inherited from DFunction class
self._add_me(self.axis, cfce)
# update reorganization energy
self.lamb += lamb
# check temperature and update cutoff time
self._set_temperature_and_cutoff_time(temperature, 5.0 * ctime)
def _make_CP29_spectral_density(self, params: dict, values: Any = None) -> None:
from .spectraldensities import SpectralDensity
temperature = params["T"]
ctime = params["gamma"]
# omega = params["freq"]
# use the units in which params was defined
lamb = self.manager.iu_energy(params["reorg"], units=self.energy_units)
time = self.axis # .data
if values is not None:
cfce = values
else:
# Make it via SpectralDensity
fa = SpectralDensity(time, params)
cf = fa.get_CorrelationFunction(temperature=temperature)
cfce = cf.data
# this is a call to the function inherited from DFunction class
self._add_me(self.axis, cfce)
# update reorganization energy
self.lamb += lamb
# check temperature and update cutoff time
self._set_temperature_and_cutoff_time(temperature, 5.0 * ctime)
def _make_value_defined(self, params: dict, values: Any) -> None:
lamb = params["reorg"]
temperature = params["T"]
if "cutoff-time" in params.keys():
ctime = params["cutoff-time"]
else:
ctime = self.axis.max
if values is not None:
if len(values) == self.axis.length:
cfce = values
else:
raise QuantarheiError("Incompatible values")
else:
raise QuantarheiError("Valued-defined correlation function without values")
# this is a call to the function inherited from DFunction class
self._add_me(self.axis, cfce)
# update reorganization energy
self.lamb += lamb
# check temperature and update cutoff time
self._set_temperature_and_cutoff_time(temperature, ctime)
#
# Aritmetic operations
#
def __add__(self, other: DFunction) -> CorrelationFunction:
"""Addition of two correlation functions"""
if not isinstance(other, CorrelationFunction):
raise TypeError("Can only add CorrelationFunction to CorrelationFunction")
t1 = self.axis
t2 = other.axis
if t1 == t2:
# f = CorrelationFunction(t1, params=self.params)
# f.add_to_data(other)
f = deepcopy(self) # Works with value defined functions
f.add_to_data2(other)
# with energy_units("int"):
# f = CorrelationFunction(t1, params=self.params)
# #f = self.deepcopy()
# f.add_to_data(other)
else:
raise QuantarheiError(
"In addition, functions have to share the same TimeAxis object"
)
return f
def __iadd__(self, other: DFunction) -> CorrelationFunction:
"""Inplace addition of two correlation functions"""
if not isinstance(other, CorrelationFunction):
raise TypeError("Can only add CorrelationFunction to CorrelationFunction")
self.add_to_data2(other)
return self
[docs]
def add_to_data(self, other: CorrelationFunction) -> None:
"""Addition of data from a specified CorrelationFunction to this object"""
t1 = self.axis
t2 = other.axis
if t1 == t2:
self.data += other.data
self.lamb += other.lamb # reorganization energy is additive
if other.cutoff_time > self.cutoff_time:
self.cutoff_time = other.cutoff_time
if self.temperature != other.temperature:
raise QuantarheiError(
"Cannot add two correlation functions on different temperatures"
)
for p in other.params:
self.params.append(p)
self._is_composed = True
self._is_empty = False
else:
raise QuantarheiError(
"In addition, functions have to share the same TimeAxis object"
)
[docs]
def add_to_data2(self, other: CorrelationFunction) -> None:
"""Addition of data from a specified CorrelationFunction to this object"""
if self == other:
# print(other.energy_units)
# print(other.params)
with energy_units("int"): # other.energy_units):
ocor = CorrelationFunction(other.axis, other.params)
else:
ocor = other
t1 = self.axis
t2 = ocor.axis
if t1 == t2:
self.data += ocor.data
self.lamb += ocor.lamb # reorganization energy is additive
if ocor.cutoff_time > self.cutoff_time:
self.cutoff_time = ocor.cutoff_time
if self.temperature != ocor.temperature:
raise QuantarheiError(
"Cannot add two correlation functions on different temperatures"
)
for p in ocor.params:
self.params.append(p)
self._is_composed = True
self._is_empty = False
else:
raise QuantarheiError(
"In addition, functions have to share the same TimeAxis object"
)
[docs]
def reorganization_energy_consistent(self, rtol: float = 1.0e-3) -> bool:
"""Checks if the reorganization energy is consistent with the data
Calculates reorganization energy from the data and checks if it
is within specified tolerance from the expected value
"""
lamb1 = self.measure_reorganization_energy()
lamb2 = self.convert_energy_2_current_u(self.lamb)
if (abs(lamb1 - lamb2) / (lamb1 + lamb2)) < rtol:
return True
return False
[docs]
def is_analytical(self) -> bool:
"""Returns `True` if analytical
Returns `True` if the CorrelationFunction object is constructed
by analytical formula. Returns `False` if the object was constructed
by numerical transformation from spectral density.
"""
return bool(cast(dict[str, Any], self.params)["ftype"] in self.analytical_types)
[docs]
def get_temperature(self) -> float | numpy.ndarray:
"""Returns the temperature of the correlation function"""
return self.temperature
[docs]
def get_reorganization_energy(self) -> float | numpy.ndarray:
"""Returns the reorganization energy of the correlation function"""
return self.convert_energy_2_current_u(self.lamb)
[docs]
def get_correlation_time(self) -> float | numpy.ndarray:
"""Returns correlation time associated with the first component
of the bath correlation function
"""
return self.params[0]["cortime"]
[docs]
def measure_reorganization_energy(self) -> float | numpy.ndarray:
"""Calculates the reorganization energy of the correlation function
Calculates the reorganization energy of the correlation function by
integrating its imaginary part.
"""
# with energy_units("int"):
primitive = c2h(self.axis, self.data)
lamb = -numpy.imag(primitive[self.axis.length - 1])
return self.convert_energy_2_current_u(lamb)
[docs]
def copy(self) -> CorrelationFunction:
"""Creates a copy of the current correlation function"""
with energy_units("int"):
cfce = CorrelationFunction(self.axis, self.params)
return cfce
[docs]
def get_SpectralDensity(self, fa: FrequencyAxis | None = None) -> Any:
"""Returns a SpectralDensity corresponding to this CorrelationFunction.
If a FrequencyAxis object is included, the SpectralDensity
object will be returned with that FrequencyAxis instance as its
frequency axis.
"""
from .spectraldensities import SpectralDensity
# protect this from external units
with energy_units("int"):
frequencies = self.axis.get_FrequencyAxis()
vals = self.get_OddFTCorrelationFunction().data
if fa is not None:
if numpy.all(numpy.isclose(fa.data, frequencies.data, 1e-5)):
# time = ta
# pass
frequencies = fa
else:
raise QuantarheiError(
"The provided FrequencyAxis does not "
"have the same data as the Fourier "
"transformed axis"
)
# FIXME: how to set the limit of SpectralDensity at w->0
spectd = SpectralDensity(frequencies, self.params, values=vals)
return spectd
[docs]
def get_FTCorrelationFunction(self) -> FTCorrelationFunction:
"""Returns a Fourier transform of the correlation function
Returns a Fourier transform of the correlation function in form
of an instance of a special class ``FTCorrelationFunction``
"""
with energy_units("int"):
ftcf = FTCorrelationFunction(self.axis, self.params)
return ftcf
[docs]
def get_OddFTCorrelationFunction(self) -> OddFTCorrelationFunction:
"""Returns a odd part of the Fourier transform of correlation function
Returns the odd part of a Fourier transform of the correlation
function in form of an instance of a special class
``OddFTCorrelationFunction``
"""
with energy_units("int"):
oftcf = OddFTCorrelationFunction(self.axis, self.params)
return oftcf
[docs]
def get_EvenFTCorrelationFunction(self) -> EvenFTCorrelationFunction:
"""Returns a even part of the Fourier transform of correlation function
Returns the even part of a Fourier transform of the correlation
function in form of an instance of a special class
``EvenFTCorrelationFunction``
"""
with energy_units("int"):
eftcf = EvenFTCorrelationFunction(self.axis, self.params)
return eftcf
[docs]
class LineshapeFunction(DFunction, UnitsManaged):
"""Lineshape function obtained by double integration of a correlation function.
Computes ``g(t) = int_0^t dt' int_0^{t'} dt'' C(t'')`` numerically,
where ``C`` is a :class:`CorrelationFunction` built from ``params``.
Parameters
----------
axis : TimeAxis
Time grid on which the correlation function is defined.
params : dict
Parameter dictionary for the underlying :class:`CorrelationFunction`.
Must contain at least ``"ftype"``, ``"reorg"``, and ``"T"``.
values : numpy.ndarray, optional
Pre-computed correlation function values; passed through to
:class:`CorrelationFunction`. Default is ``None``.
lfactor : int, optional
Length multiplication factor applied to the time axis before
integration. Default is ``1``.
Raises
------
Exception
If ``params`` is ``None``.
"""
@enforce_energy_units_context
def __init__(
self, axis: Any = None, params: Any = None, values: Any = None, lfactor: int = 1
) -> None:
self.lfactor = lfactor
self.axis = axis
self.params = params
self.values = values
if params is None:
raise QuantarheiError("Argument 'params' must not be None")
# FIXME: This DOES NOT WORK!!!, we need something else to distinguish
try:
N = len(params)
except TypeError:
N = 1
# print(N, params)
N = 1
if N == 1:
tt = TimeAxis(axis.start, lfactor * axis.length, axis.step)
cf = CorrelationFunction(axis=tt, params=params, values=values)
gg = c2g(tt, cf.data)
else:
k = 0
for prm in params:
tt = TimeAxis(axis.start, lfactor * axis.length, axis.step)
cf = CorrelationFunction(axis=tt, params=prm, values=values)
if k == 0:
cf_out = cf
else:
cf_out.__iadd__(cf)
k += 1
gg = c2g(tt, cf_out.data)
super().__init__(tt, gg)
[docs]
class FTCorrelationFunction(DFunction, UnitsManaged):
"""Fourier transform of the correlation function
Numerically calculated Fourier transform of the correlation function
Parameters
----------
axis: TimeAxis
Time interval from which the frequency interval is calculated
params: dictionary
Dictionary of the correlation function parameters
"""
energy_params = ("reorg", "omega", "freq")
def __init__(self, axis: Any, params: Any, values: Any = None) -> None:
super().__init__()
if not isinstance(axis, FrequencyAxis):
faxis = axis.get_FrequencyAxis()
self.axis = faxis
else:
self.axis = axis
# handle params
self.params: (
list[Any] | dict[Any, Any]
) = [] # this will always be a list of components
p2calc = []
try:
# if this passes, we assume params is a dictionary
params.keys()
self._is_composed = False
p2calc.append(params)
except AttributeError:
# othewise we assume it is a list of dictionaries
self._is_composed = True
for p in params:
p2calc.append(p)
for params in p2calc:
try:
ftype = params["ftype"]
if ftype not in CorrelationFunction.allowed_types:
raise QuantarheiError("Unknown CorrelationFunction type")
# we mutate the parameters into internal units
prms = {}
for key in params.keys():
if key in self.energy_params:
prms[key] = self.convert_energy_2_internal_u(params[key])
else:
prms[key] = params[key]
except KeyError:
raise QuantarheiError(
"Dictionary of parameters does not contain `ftype` key"
)
self.params = prms
if values is None:
# data have to be protected from change of units
with energy_units("int"):
cfce = CorrelationFunction(axis, self.params)
ftvals = cfce.get_Fourier_transform()
self.data = ftvals.data
# self.axis = ftvals.axis
else:
# This is not protected from change of units!!!!
self.data = values
# self.axis = cfce.axis.get_FrequencyAxis()
[docs]
class OddFTCorrelationFunction(DFunction, UnitsManaged):
"""Odd part of the Fourier transform of the correlation function
Numerically calculated odd part Fourier transform of the correlation
function. Calculated as Fourier transform of the imaginary part of the
correlation function.
Parameters
----------
axis: TimeAxis
Time interval from which the frequency interval is calculated
params: dictionary
Dictionary of the correlation function parameter
Examples
--------
>>> ta = TimeAxis(0.0,1000,1.0)
>>> params = dict(ftype="OverdampedBrownian",reorg=20,cortime=100,T=300)
>>> with energy_units("1/cm"):
... ocf = OddFTCorrelationFunction(ta,params)
... print(numpy.allclose(ocf.at(-100), -ocf.at(100)))
True
"""
def __init__(self, axis: Any, params: Any, values: Any = None) -> None:
super().__init__()
if not isinstance(axis, FrequencyAxis):
faxis = axis.get_FrequencyAxis()
self.axis = faxis
else:
self.axis = axis
# handle params
self.params = [] # this will always be a list of components
p2calc = []
try:
# if this passes, we assume params is a dictionary
params.keys()
self._is_composed = False
p2calc.append(params)
except AttributeError:
# othewise we assume it is a list of dictionaries
self._is_composed = True
for p in params:
p2calc.append(p)
for params in p2calc:
ftype = params["ftype"]
if ftype not in CorrelationFunction.allowed_types:
raise QuantarheiError("Unknown Correlation Function Type")
self.params.append(params)
# We create CorrelationFunction and FTT it
if params["ftype"] == "Value-defined":
if values is None:
raise QuantarheiError()
else:
cfce = CorrelationFunction(axis, params, values=values)
else:
cfce = CorrelationFunction(axis, params)
cfce.data = 1j * numpy.imag(cfce.data)
# data have to be protected from change of units
with energy_units("int"):
ftvals = cfce.get_Fourier_transform()
ndata = numpy.real(ftvals.data)
self._add_me(self.axis, ndata)
[docs]
class EvenFTCorrelationFunction(DFunction, UnitsManaged):
"""Even part of the Fourier transform of the correlation function
Numerically calculated even part Fourier transform of the correlation
function. Calculated as Fourier transform of the real part of the
correlation function.
Parameters
----------
axis: TimeAxis
Time interval from which the frequency interval is calculated
params: dictionary
Dictionary of the correlation function parameter
Examples
--------
>>> ta = TimeAxis(0.0,1000,1.0)
>>> params = dict(ftype="OverdampedBrownian",reorg=20,cortime=100,T=300)
>>> with energy_units("1/cm"):
... ecf = EvenFTCorrelationFunction(ta,params)
... print(numpy.allclose(ecf.at(-100), ecf.at(100)))
True
"""
def __init__(self, axis: Any, params: Any, values: Any = None) -> None:
super().__init__()
if not isinstance(axis, FrequencyAxis):
faxis = axis.get_FrequencyAxis()
self.axis = faxis
else:
self.axis = axis
# handle params
self.params = [] # this will always be a list of components
p2calc = []
try:
# if this passes, we assume params is a dictionary
params.keys()
self._is_composed = False
p2calc.append(params)
except AttributeError:
# othewise we assume it is a list of dictionaries
self._is_composed = True
for p in params:
p2calc.append(p)
for params in p2calc:
ftype = params["ftype"]
if ftype not in CorrelationFunction.allowed_types:
raise QuantarheiError("Unknown Correlation Function Type: " + ftype)
self.params.append(params)
# We create CorrelationFunction and FTT it
if params["ftype"] == "Value-defined":
if values is None:
raise QuantarheiError()
else:
cfce = CorrelationFunction(axis, params, values=values)
else:
cfce = CorrelationFunction(axis, params)
cfce.data = numpy.real(cfce.data)
# data have to be protected from change of units
with energy_units("int"):
ftvals = cfce.get_Fourier_transform()
ndata = numpy.real(ftvals.data)
self._add_me(self.axis, ndata)
# FIXME: these functions can go to DFunction
[docs]
def c2g(timeaxis: TimeAxis, coft: numpy.ndarray) -> numpy.ndarray:
"""Converts correlation function to lineshape function
Explicit numerical double integration of the correlation
function to form a lineshape function.
Parameters
----------
timeaxis : cu.oqs.time.TimeAxis
TimeAxis of the correlation function
coft : complex numpy array
Values of correlation function given at points specified
in the TimeAxis object
"""
time = timeaxis
preal = numpy.real(coft)
pimag = numpy.imag(coft)
splr = interp.UnivariateSpline(time.data, preal, s=0).antiderivative()(time.data)
splr = interp.UnivariateSpline(time.data, splr, s=0).antiderivative()(time.data)
spli = interp.UnivariateSpline(time.data, pimag, s=0).antiderivative()(time.data)
spli = interp.UnivariateSpline(time.data, spli, s=0).antiderivative()(time.data)
goft = splr + 1j * spli
return goft
[docs]
def c2h(timeaxis: TimeAxis, coft: numpy.ndarray) -> numpy.ndarray:
"""Integrates correlation function in time with an open upper limit
Explicit numerical integration of the correlation
function to form a precursor to the lineshape function.
Parameters
----------
timeaxis : TimeAxis
TimeAxis of the correlation function
coft : complex numpy array
Values of correlation function given at points specified
in the TimeAxis object
"""
time = timeaxis
preal = numpy.real(coft)
pimag = numpy.imag(coft)
splr = interp.UnivariateSpline(time.data, preal, s=0).antiderivative()(time.data)
spli = interp.UnivariateSpline(time.data, pimag, s=0).antiderivative()(time.data)
hoft = splr + 1j * spli
return hoft
[docs]
def h2g(timeaxis: TimeAxis, coft: numpy.ndarray) -> numpy.ndarray:
"""Integrates and integrated correlation function
Explicit numerical integration of the correlation
function to form a precursor to the lineshape function.
Parameters
----------
timeaxis : TimeAxis
TimeAxis of the correlation function
coft : complex numpy array
Values of correlation function given at points specified
in the TimeAxis object
"""
return c2h(timeaxis, coft)
[docs]
def oscillator_scalled_CorrelationFunction(
time: TimeAxis,
params: dict,
omega: float,
target_time: float,
Nmax: int = 5,
HR: float = 0.01,
silent: bool = True,
) -> CorrelationFunction:
"""Scales reorganization energy of the system-bath interaction
Returns a bath correlation function with reorganization energy scalled
such that it achieves a required relaxation time between the first
vibrationally excited and the ground state of an oscillator with a
given frequency.
Parameters
----------
time: TimeAxis
The time axis on which the correlation function is defined
params: dictionary
A dictionary of the correlation function parameters. Requires
the "reorg" item.
omega: float
Frequency/energy of the oscillator
target_time: float
The requested relaxation time in fs
Nmax: int
The mumber of harmonic oscillator levels used for the calculation
Default is 5.
HR: float
Huang-Rhys factor to be used in the calculation. Default is 0.1.
silent: bool
If silent is not True, some info about internal calculations
is printed.
Examples
--------
>>> omega = 200.0
>>> target_time = 100.0
Defining correltion function and time axis
>>> time = TimeAxis(0.0, 1000, 1.0)
>>> params = dict(ftype="OverdampedBrownian", T=300, cortime=30.0,
... reorg=100.0, matsubara=30)
>>> with energy_units("1/cm"):
... cf = oscillator_scalled_CorrelationFunction(time, params, omega,
... target_time, Nmax=7, HR=0.3)
>>> print("%10.5f" % params["reorg"])
32.97806
"""
# from ..qm.corfunctions import CorrelationFunction
from ...builders.modes import Mode
from ...builders.molecules import Molecule
from ...core.managers import eigenbasis_of
Nmx = Nmax
cf = CorrelationFunction(time, params)
with energy_units("1/cm"):
mol = Molecule([0.0, 10000.0])
mol.set_electronic_rwa([0, 1])
md = Mode(omega)
mol.add_Mode(md)
md.set_HR(1, HR)
md.set_nmax(0, Nmx)
md.set_nmax(1, Nmx)
# mol.set_mode_environment(mode=0, elstate=0, corfunc=cf)
# mol.set_mode_environment(mode=0, elstate=1, corfunc=cf)
mol.set_mode_environment(mode=0, elstate="ALL", corfunc=cf)
mol.set_transition_environment((0, 1), cf)
# HH = mol.get_Hamiltonian()
# sbi = mol.get_SystemBathInteraction()
time = cf.axis
RR, ham = mol.get_RelaxationTensor(
time, relaxation_theory="stR", as_operators=False
)
ref_time = 1.0 / numpy.real(RR.data[0, 0, 1, 1])
# print("Ref. time:", ref_time)
ratio = ref_time / target_time
# print("Ratio:", ratio)
orig_reorg = params["reorg"]
reorg = ratio * orig_reorg
params["reorg"] = reorg
cfn = CorrelationFunction(time, params)
if not silent:
with eigenbasis_of(ham):
print("Original time (S0):", 1.0 / numpy.real(RR.data[0, 0, 1, 1]))
print(
"Original time (S1):",
1.0 / numpy.real(RR.data[Nmx, Nmx, Nmx + 1, Nmx + 1]),
)
return cfn