"""Redfield Tensor
Class Details
-------------
"""
from __future__ import annotations
from typing import Any
# import time
import numpy
import scipy
from ... import COMPLEX, REAL
from ...core.managers import energy_units
from ...exceptions import QuantarheiError
from ...utils.logging import log_detail, log_quick
# from ...core.managers import BasisManaged
from ...utils.types import BasisManagedComplexArray
from ..hilbertspace.hamiltonian import Hamiltonian
from .relaxationtensor import RelaxationTensor
from .systembathinteraction import SystemBathInteraction
[docs]
class RedfieldRelaxationTensor(RelaxationTensor):
"""Redfield Relaxation Tensor
Unlike Foerster or Lindblad tensors (which operate entirely in the site
basis), the Redfield tensor is computed in the eigenbasis of the
Hamiltonian. The constructor diagonalizes the Hamiltonian internally and
stores the result in that eigenbasis.
If you obtain this tensor via ``OpenSystem.get_RelaxationTensor()``, it is
automatically back-transformed to the site basis before being returned —
no extra steps are needed.
If you construct it directly, you must back-transform explicitly::
RRT = RedfieldRelaxationTensor(ham, sbi, secular=True)
_, SS = numpy.linalg.eigh(ham._data)
S1 = numpy.linalg.inv(SS)
RRT.transform(S1, inv=SS)
Parameters
----------
ham : Hamiltonian
Hamiltonian of the system.
sbi : SystemBathInteraction
Object specifying the system-bath interaction
initialize : bool
If True, the tensor will be immediately calculated.
cutoff_time : float
Time in femtoseconds after which the integral kernel in the
definition of the relaxation tensor is assumed to be zero.
as_operators : bool
If True the tensor will not be constructed. Instead a set of
operators whose application is equal to the application of the
tensor will be defined and stored.
secular : bool
If True, the tensor is secularized after construction.
"""
# data = BasisManagedComplexArray("data")
Km = BasisManagedComplexArray("Km")
Lm = BasisManagedComplexArray("Lm")
Ld = BasisManagedComplexArray("Ld")
_has_cutoff_time: bool
cutoff_time: float | None
_is_initialized: bool
dim: int
_Km: numpy.ndarray
_Lm: numpy.ndarray
_Ld: numpy.ndarray
def __init__(
self,
ham: Hamiltonian,
sbi: SystemBathInteraction,
initialize: bool = True,
cutoff_time: float | None = None,
as_operators: bool = False,
name: str = "",
secular: bool = False,
) -> None:
self._initialize_basis()
#
# Check the types of the arguments
#
# Hamiltonian
if not isinstance(ham, Hamiltonian):
raise QuantarheiError("First argument must be a Hamiltonian")
# SystemBathInteraction
if not isinstance(sbi, SystemBathInteraction):
if sbi is not None:
raise QuantarheiError(
"Second argument must be of the SystemBathInteraction type"
)
self._secular_on_init = secular
self.Hamiltonian = ham
self.SystemBathInteraction: SystemBathInteraction = sbi
self.dim = self.Hamiltonian.dim
self.name = name
self._data_initialized = False
self._is_initialized = False
self._has_cutoff_time = False
self.as_operators = as_operators
if not self.as_operators:
self.data = numpy.zeros(
(self.dim, self.dim, self.dim, self.dim), dtype=numpy.complex128
)
# set cut-off time
if cutoff_time is not None:
self.cutoff_time = cutoff_time
self._has_cutoff_time = True
# initialize the tensor right at creation
if initialize:
try:
pass
except Exception:
pass
with energy_units("int"):
self._implementation(ham, sbi)
if self._secular_on_init:
self.secularize(legacy=False)
self.Iterm = None
self.has_Iterm = False
[docs]
def apply(self, oper: Any, target: Any = None, copy: bool = True) -> Any:
"""Applies the relaxation tensor on a superoperator"""
if self.as_operators:
print("Applying Relaxation tensor")
if copy:
import copy as _copy
oper_ven = _copy.copy(oper)
else:
oper_ven = oper
rho1 = oper.data
# apply tensor to data
Lm = self.Lm
Km = self.Km
Ld = self.Ld
Kd = numpy.zeros(Km.shape, dtype=numpy.float64)
Nm = Km.shape[0]
ven = numpy.zeros(oper.data.shape, dtype=numpy.complex128)
for mm in range(Nm):
Kd[mm, :, :] = numpy.transpose(Km[mm, :, :])
ven += (
numpy.dot(Km[mm, :, :], numpy.dot(rho1, Ld[:, :, mm]))
+ numpy.dot(Lm[:, :, mm], numpy.dot(rho1, Kd[mm, :, :]))
- numpy.dot(numpy.dot(Kd[mm, :, :], Lm[:, :, mm]), rho1)
- numpy.dot(rho1, numpy.dot(Ld[:, :, mm], Km[mm, :, :]))
)
oper_ven.data = ven
return oper_ven
return super().apply(oper, copy=copy)
[docs]
def get_population_rate(self, N: int, M: int) -> Any:
"""Returns the relaxation rate between states N -> M"""
if self.as_operators:
rt = 0.0
for ii in range(self.Lm.shape[2]):
rt += 2.0 * numpy.real(self.Lm[N, M, ii] * self.Km[ii, N, M])
return rt
return super().get_population_rate(N, M)
[docs]
def get_dephasing_rate(self, N: int, M: int) -> Any:
"""Returns the dephasing rate of a coherence between states N and M"""
if self.as_operators:
rt = 0.0
for ii in range(self.Lm.shape[2]):
A_a = numpy.dot(self.Km[ii, N, :], self.Lm[:, N, ii])
A_b = numpy.dot(self.Ld[M, :, ii], self.Km[ii, :, M])
KL = self.Km[ii, N, N] * self.Ld[M, M, ii]
LK = self.Lm[N, N, ii] * self.Km[ii, M, M]
rt += A_a + A_b - KL - LK
return rt
return super().get_dephasing_rate(N, M)
def _implementation(self, ham: Hamiltonian, sbi: SystemBathInteraction) -> None:
r"""Reference implementation, completely in Python
Implementation of Redfield relaxation tensor according to
V. May and O. Kuehn, Charge and Energy Transfer Dynamics in Molecular
System, Wiley-VCH, Berlin, 2000, 1st edition, chapter 3.8.2.
In particular we refer to Eq. (3.8.13) on page 132
We assume the system-bath interaction operator in form of Eq. (3.5.30)
with the bath part specified through two-point correlation functions.
We construct operators K_{m} introduced in Eq. (3.5.30) and
operators \Lambda_{m} of Eq. (3.8.11).
We do not delete the imaginary part of the tensor (as is done later
in Section 3.8.3 to get the so-called "Multi-level Redfield
Equations"). Such a deletion can be done later manually.
"""
log_detail("Reference time-independent Redfield tensor calculation")
# print("Reference Redfield implementation ...")
#
# dimension of the Hamiltonian (includes excitons
# with all multiplicities specified at its creation)
#
Na = ham.dim # data.shape[0]
# time axis
ta = sbi.TimeAxis
dt = ta.step
#
# is this beyond single excitation band?
#
multi_ex = False
# figure out if the aggregate or molecule specifies more than
# one exciton band
if sbi.aggregate is not None:
agg = sbi.aggregate
if agg.mult > 1:
multi_ex = True
elif sbi.molecule is not None:
mol = sbi.molecule
if mol.mult > 1:
multi_ex = True
#
# shorten the interval of integration if a cut-off time is set
#
if self._has_cutoff_time:
# index of the cut-off time on the time axis
tcut = ta.nearest(self.cutoff_time)
# select the section of the time axis up to the cut-off time
tm = ta.data[0:tcut]
# length of the section corresponds to the index of cut-off time
length = tcut
else:
# if cut-off time is not set then we take the whole time axis
tm = ta.data
# and the length corresponds to the length of the time axis
length = ta.length
#
# Get eigenenergies and transformation matrix of the Hamiltonian
#
# Use _data directly: _implementation must always diagonalize in the
# site basis so that Km is the correct eigenbasis coupling operator.
hD, SS = numpy.linalg.eigh(ham._data)
#
# Find all transition frequencies
#
Om = numpy.zeros((Na, Na))
for a in range(Na):
for b in range(Na):
Om[a, b] = hD[a] - hD[b]
# number of baths - one per monomer
Nb = sbi.N
#
# Site K_m operators
#
Nex = Nb
if multi_ex:
Nex = 2 * Nb
Km = numpy.zeros((Nex, Na, Na), dtype=numpy.float64)
# Transform site operators
S1 = scipy.linalg.inv(SS)
# FIXME: SBI should also be basis controlled
for ns in range(Nb):
Km[ns, :, :] = numpy.dot(S1, numpy.dot(sbi.KK[ns, :, :], SS))
if multi_ex:
system_any: Any = sbi.system
try:
ii = system_any.twoex_state[0, 0]
# print(ii)
# if ii >= 0:
# print("Something is wrong")
# raise QuantarheiError("twoex_state property ill-defined.")
except AttributeError:
raise QuantarheiError("System requires the twoex_state property.")
# There are also two-exciton bath which can also be converted into one per monomer
Nb2 = sbi.N
# projectors to two-exciton states
K2: numpy.ndarray = numpy.zeros((Nb2, Na, Na), dtype=REAL)
i_start = ham.rwa_indices[2] # here the two-exciton bands starts
ii = 0
for ms in range(i_start, ham.dim):
K2[ii, ms, ms] = 1.0
ii += 1
# here we do two-exciton projectors
for ns in range(Nb):
for ms in range(Nb):
if ms != ns:
n2ex = (
system_any.twoex_state[ms, ns] - Nb - 1
) # we number the K2 projectors from zero
Km[Nb + ns, :, :] += numpy.dot(
S1, numpy.dot(K2[n2ex, :, :], SS)
)
#
# \Lambda_m operator
#
# Integrals of correlation functions from the set
if not multi_ex:
Lm = numpy.zeros((Na, Na, Nb), dtype=numpy.complex128)
else:
Lm = numpy.zeros((Na, Na, 2 * Nb), dtype=numpy.complex128)
Nex = 2 * Nb
for ms in range(Nb):
log_quick("Calculating bath component", ms, "of", Nb, end="\r")
# print(ms, "of", Nb)
# for ns in range(Nb):
ns = ms
# correlation function of site ns (if ns == ms)
# or a cross-correlation function of sites ns and ms
# FIXME: reaching correct correlation function is a nightmare!!!
rc1 = sbi.CC.get_coft(ms, ns)
self._guts_Cmplx_Splines(ms, Lm, Km, Na, Om, length, rc1, tm)
log_quick()
eexp = numpy.exp(
-1.0j * Om[:, :, numpy.newaxis] * tm[numpy.newaxis, numpy.newaxis, :]
)
if multi_ex:
for ms in range(Nb):
ns = ms
# print("two-ex", ns, ":")
rc1 = sbi.get_coft(ms + 1, ns + 1)
rc = numpy.einsum("t,ijt->ijt", rc1[0:length], eexp)
cc_mnab = _integrate_last_axis(rc, dt)
Lm[:, :, Nb + ms] = numpy.einsum(
"ij,ij->ij", cc_mnab[:, :, -1], Km[Nb + ms, :, :]
)
# create the Hermite conjuged version of \Lamnda_m
if not multi_ex:
Ld = numpy.zeros((Na, Na, Nb), dtype=numpy.complex128)
Nex = Nb
else:
Ld = numpy.zeros((Na, Na, 2 * Nb), dtype=numpy.complex128)
Nex = 2 * Nb
for ms in range(Nex):
Ld[:, :, ms] += numpy.conj(numpy.transpose(Lm[:, :, ms]))
self._post_implementation(Km, Lm, Ld)
log_detail("... Redfield done")
def _post_implementation(
self, Km: numpy.ndarray, Lm: numpy.ndarray, Ld: numpy.ndarray
) -> None:
"""When components of the tensor are calculated, should they be
saved or converted into full tensor?
"""
if self.as_operators:
# save the operators - propagation methods must know about them
self.Km = Km
self.Lm = Lm
self.Ld = Ld
else:
# save the relaxation tensor
RR = self._convert_operators_2_tensor(Km, Lm, Ld)
self.data = RR
self._data_initialized = True
self._is_initialized = True
def _guts_Cmplx_Splines(
self,
ms: int,
Lm: numpy.ndarray,
Km: numpy.ndarray,
Na: int,
Om: numpy.ndarray,
length: int,
rc1: numpy.ndarray,
tm: numpy.ndarray,
) -> None:
for a in range(Na):
for b in range(Na):
# argument of the integration
eexp = numpy.exp(-1.0j * Om[a, b] * tm)
rc = rc1[0:length] * eexp
# spline integration instead of FFT
rr = numpy.real(rc)
ri = numpy.imag(rc)
sr = scipy.interpolate.UnivariateSpline(tm, rr, s=0).antiderivative()(
tm
)
si = scipy.interpolate.UnivariateSpline(tm, ri, s=0).antiderivative()(
tm
)
# we take the last value (integral to infinity)
cc_mnab = sr[length - 1] + 1.0j * si[length - 1]
# \Lambda_m operators
Lm[a, b, ms] += cc_mnab * Km[ms, a, b]
def _guts_Cmplx_Sum(
self,
ms: int,
Lm: numpy.ndarray,
Km: numpy.ndarray,
Na: int,
Om: numpy.ndarray,
length: int,
rc1: numpy.ndarray,
tm: numpy.ndarray,
) -> None:
import scipy
dt = tm[1] - tm[0]
for a in range(Na):
for b in range(Na):
# argument of the integration
eexp = numpy.exp(-1.0j * Om[a, b] * tm)
rc = rc1[0:length] * eexp
# we take integral to infinity
# cc_mnab = numpy.sum(rc)*dt
cc_mnab = scipy.integrate.trapz(rc, dx=dt)
# \Lambda_m operators
Lm[a, b, ms] += cc_mnab * Km[ms, a, b]
def _guts_Cmplx_FFT(
self,
ms: int,
Lm: numpy.ndarray,
Km: numpy.ndarray,
Na: int,
Om: numpy.ndarray,
length: int,
rc1: numpy.ndarray,
tm: numpy.ndarray,
) -> None:
for a in range(Na):
for b in range(Na):
# argument of the integration
eexp = numpy.exp(-1.0j * Om[a, b] * tm)
rc = rc1[0:length] * eexp
# spline integration instead of FFT
rr = numpy.real(rc)
ri = numpy.imag(rc)
sr = scipy.interpolate.UnivariateSpline(tm, rr, s=0).antiderivative()(
tm
)
si = scipy.interpolate.UnivariateSpline(tm, ri, s=0).antiderivative()(
tm
)
# we take the last value (integral to infinity)
cc_mnab = sr[length - 1] + 1.0j * si[length - 1]
# \Lambda_m operators
Lm[a, b, ms] += cc_mnab * Km[ms, a, b]
def _convert_operators_2_tensor(
self, Km: numpy.ndarray, Lm: numpy.ndarray, Ld: numpy.ndarray
) -> numpy.ndarray:
r"""Converts operator representation to the tensor one
Convertes operator representation of the Redfield tensor
into a truely tensor representation
Parameters
----------
Km : 3D array
K_m operators
Lm : 3D array
\Lambda_m operators
Ld : 3D array
Hermite conjuget \Lambda_m operators
"""
Na = self.Hamiltonian.data.shape[0]
Nb = Lm.shape[2] # self.SystemBathInteraction.N
RR = numpy.zeros((Na, Na, Na, Na), dtype=numpy.complex128)
# from ...implementations.cython.loopit import loopit
# tt1 = time.time()
for m in range(Nb):
Kd = numpy.transpose(Km[m, :, :])
# KdLm = numpy.dot(Kd,Lm[m,:,:])
# LdKm = numpy.dot(Ld[m,:,:],Km[m,:,:])
# for a in range(Na):
# print(a)
# for b in range(Na):
# for c in range(Na):
# for d in range(Na):
#
# RR[a,b,c,d] += (Km[m,a,c]*Ld[m,d,b]
# + Lm[m,a,c]*Kd[d,b])
# if b == d:
# RR[a,b,c,d] -= KdLm[a,c]
# if a == c:
# RR[a,b,c,d] -= LdKm[d,b]
_loopit(Km, Kd, Lm, Ld, Na, RR, m)
# tt2 = time.time()
# print(tt2-tt1)
return RR
[docs]
def initialize(self) -> None:
"""Initializes the Redfield tensor with values"""
self._implementation(self.Hamiltonian, self.SystemBathInteraction)
[docs]
def convert_2_tensor(self) -> None:
"""Converts internally the operator representation to a tensor one
Converst the representation of the relaxation tensor through a set
of operators into a tensor representation.
"""
if self.as_operators:
RR = self._convert_operators_2_tensor(self.Km, self.Lm, self.Ld)
self.data = RR
self._data_initialized = True
self.as_operators = False
def _loopit(
Km: numpy.ndarray,
Kd: numpy.ndarray,
Lm: numpy.ndarray,
Ld: numpy.ndarray,
Na: int,
RR: numpy.ndarray,
m: int,
) -> None:
KdLm = numpy.dot(Kd, Lm[:, :, m])
LdKm = numpy.dot(Ld[:, :, m], Km[m, :, :])
for a in range(Na):
for b in range(Na):
for c in range(Na):
for d in range(Na):
RR[a, b, c, d] += Km[m, a, c] * Ld[d, b, m] + Lm[a, c, m] * Kd[d, b]
if b == d:
RR[a, b, c, d] -= KdLm[a, c]
if a == c:
RR[a, b, c, d] -= LdKm[d, b]
def _integrate_last_axis(f: numpy.ndarray, dt: float) -> numpy.ndarray:
"""Cumulative Simpson integration over the last axis of a 2D or higher NumPy array.
Parameters:
f : ndarray
Input array with shape (..., N), where integration is over the last axis (time).
dt : float
Uniform time step between samples.
Returns:
F : ndarray
Cumulative integral of f along the last axis, with the same shape as f.
"""
if f.ndim < 2:
raise ValueError("Input array must be at least 2D (e.g. shape (i, t))")
N = f.shape[-1]
if N < 2:
raise ValueError("Need at least two points along the last axis")
F = numpy.zeros_like(f)
last = N - 1 if N % 2 == 1 else N - 2 # leave last interval for trapezoid if needed
idxs = numpy.arange(2, last + 1, 2)
for i3 in idxs:
F[..., i3] = F[..., i3 - 2] + (dt / 3) * (
f[..., i3 - 2] + 4 * f[..., i3 - 1] + f[..., i3]
)
F[..., i3 - 1] = (F[..., i3 - 2] + F[..., i3]) / 2 # midpoint
if N % 2 == 0:
# Final trapezoidal step
F[..., -1] = F[..., -2] + 0.5 * dt * (f[..., -2] + f[..., -1])
return F