"""Linear absorption spectrum
Linear absorption spectrum of a molecule or an aggregate of molecules.
"""
from __future__ import annotations
import copy
from typing import Any
import numpy
import scipy
from .. import COMPLEX
from ..builders import Aggregate, Molecule, OpenSystem
from ..core.frequency import FrequencyAxis
from ..core.managers import EnergyUnitsManaged, eigenbasis_of, energy_units
from ..core.time import TimeAxis, TimeDependent
from ..core.units import kB_intK
from ..core.wrappers import prevent_basis_context
from ..exceptions import ImplementationError, QuantarheiError
from ..qm.hilbertspace.operators import ReducedDensityMatrix
from ..utils import derived_type
from .abs2 import AbsSpectrum
from .linear_spectra import LinSpectrum
# TODO: move _c2g from the calculator and use _c2g instead self._c2g
[docs]
class LinSpectrumCalculator(EnergyUnitsManaged):
"""Linear absorption spectrum
Parameters
----------
timeaxis : quantarhei.TimeAxis
TimeAxis on which the calculation will be performed
system : quantarhei.Molecule or quantathei.Aggregate
System for which the absorption spectrum will be calculated.
Examples
--------
Calcutor has to be created using a Molecule, Aggregate or OpenSystem
>>> time = TimeAxis(0.0, 1000, 1.0)
>>> absc = AbsSpectrumCalculator(time)
Traceback (most recent call last):
...
TypeError: system must be of type [<class 'quantarhei.builders.molecules.Molecule'>, <class 'quantarhei.builders.aggregates.Aggregate'>, <class 'quantarhei.builders.opensystem.OpenSystem'>]
>>> with energy_units("1/cm"):
... mol = Molecule([0.0, 10000.0])
>>> absc = AbsSpectrumCalculator(time, mol)
The method bootstrap() must be called before calculate()
>>> abs = absc.calculate()
Traceback (most recent call last):
...
quantarhei.exceptions.QuantarheiError: Calculator must be bootstrapped first: call bootstrap() method of this object.
Molecule does not provide RWA information automatically
>>> HH = mol.get_Hamiltonian()
>>> HH.has_rwa
False
Setting the RWA frequency can be therefore done explicitely
>>> absc.bootstrap(rwa=1.0)
>>> print(absc.rwa)
1.0
When RWA is specified for the Hamiltonian
>>> HH.set_rwa([0, 1])
setting RWA explicitely in bootstrap() is ignored. The value from the
Molecule is used.
>>> absc.bootstrap(rwa=1.1)
>>> print("%5.5f" % absc.rwa)
1.88365
RWA can be set by the Molecule
>>> with energy_units("1/cm"):
... mol = Molecule([0.0, 10000.0])
>>> absc = AbsSpectrumCalculator(time, mol)
>>> mol.set_electronic_rwa([0, 1])
>>> absc.bootstrap()
"""
TimeAxis = derived_type("TimeAxis", TimeAxis)
system = derived_type("system", [Molecule, Aggregate, OpenSystem])
def __init__(
self,
timeaxis: Any,
system: Any = None,
dynamics: str = "secular",
relaxation_tensor: Any = None,
rate_matrix: Any = None,
effective_hamiltonian: Any = None,
) -> None:
# protected properties
self.TimeAxis = timeaxis
self.system = system
if dynamics in ["secular", "non-secular"]:
self.dynamics = dynamics
else:
raise QuantarheiError("Unknown type of dynamics: " + str(dynamics))
# unprotected properties
self._relaxation_tensor = None
self._rate_matrix = None
self._relaxation_hamiltonian = None
self._has_relaxation_tensor = False
self._gass_lineshape = False
self._gauss_broad = False
self._is_adiabatic = False
self._adiabatic_noBath = False
if relaxation_tensor is not None:
self._relaxation_tensor = relaxation_tensor
self._has_relaxation_tensor = True
if effective_hamiltonian is not None:
self._relaxation_hamiltonian = effective_hamiltonian
if rate_matrix is not None:
self._rate_matrix = rate_matrix
self._has_rate_matrix = True
self._pop_fluor = "vertical"
self.rwa: float | numpy.ndarray = 0.0
self.prop_has_rwa = False
self.bootstrapped = False
[docs]
def bootstrap(
self,
rwa: float = 0.0,
prop: Any = None,
lab: Any = None,
HWHH: Any = None,
axis: Any = None,
gauss: Any = None,
adiabatic: Any = None,
fluor: Any = None,
**kwargs: Any,
) -> None:
"""This function sets some additional information before calculation
>>> time = TimeAxis(0.0, 1000, 1.0)
>>> with energy_units("1/cm"):
... mol = Molecule([0.0, 10000.0])
>>> absc = AbsSpectrumCalculator(time, mol)
>>> absc.bootstrap()
Traceback (most recent call last):
...
quantarhei.exceptions.QuantarheiError: RWA not set by system nor explicitely.
"""
self.prop = prop
HH = self.system.get_Hamiltonian()
if HH.has_rwa:
# if HH has RWA, the 'rwa' argument is ignored
HR = HH.get_RWA_skeleton()
Ng = HH.rwa_indices[0]
Ne = HH.rwa_indices[1]
self.rwa = self.convert_2_internal_u(HR[Ne] - HR[Ng])
self.prop_has_rwa = True # Propagator is in RWA
else:
if rwa > 0.0:
self.rwa = self.convert_2_internal_u(rwa)
else:
raise QuantarheiError("RWA not set by system nor explicitely.")
# self.rwa = self.convert_2_internal_u(rwa)
with energy_units("int"):
# sets the frequency axis for plottig
self.frequencyAxis = self.TimeAxis.get_FrequencyAxis()
# it must be shifted by rwa value
self.frequencyAxis.data += self.rwa
if HWHH is not None:
self._gass_lineshape = True
self.HWHH = self.convert_2_internal_u(HWHH)
if gauss is not None:
self._gauss_broad = True
self.gauss = self.convert_2_internal_u(gauss)
if adiabatic is not None:
if adiabatic == False:
self._is_adiabatic = False
else:
self._is_adiabatic = True
if (adiabatic == "SubtractBath") or (adiabatic == "NoBath"):
self._adiabatic_noBath = True
if fluor is not None:
if fluor == "vertical":
self._pop_fluor = "vertical"
elif fluor == "adiabatic":
self._pop_fluor = "adiabatic"
else:
raise Warning(
"fluor keyword value unknown. Allowed values are vertical or adiabatic."
)
if axis is not None:
if axis == "x":
self.ld_axis = numpy.array([1.0, 0.0, 0.0], dtype="f8")
elif axis == "y":
self.axis = numpy.array([0.0, 1.0, 0.0], dtype="f8")
elif axis == "z":
self.ld_axis = numpy.array([0.0, 0.0, 1.0], dtype="f8")
else:
self.ld_axis = axis / numpy.linalg.norm(axis)
else:
self.ld_axis = numpy.array([0.0, 0.0, 1.0], dtype="f8")
# if isinstance(self.system, Aggregate):
# self.system.diagonalize()
self.bootstrapped = True
[docs]
@prevent_basis_context
def calculate(
self, raw: bool = False, from_dynamics: bool = False, alt: bool = False
) -> Any:
"""Calculates the absorption spectrum"""
if not self.bootstrapped:
raise QuantarheiError(
"Calculator must be bootstrapped first: "
"call bootstrap() method of this object."
)
with energy_units("int"):
if self.system is not None:
if from_dynamics:
# alt = True is for testing only
spect = self._calculate_abs_from_dynamics(raw=raw, alt=alt)
elif isinstance(self.system, Aggregate):
spect = self._calculate_aggregate(
relaxation_tensor=self._relaxation_tensor,
rate_matrix=self._rate_matrix,
relaxation_hamiltonian=self._relaxation_hamiltonian,
raw=raw,
)
elif isinstance(self.system, (Molecule, OpenSystem)):
# self._calculate_Molecule(rwa)
spect = self._calculate_monomer(raw=raw)
else:
raise QuantarheiError("System to calculate spectrum for not defined")
return spect
def _calculateMolecule(self, rwa: float) -> None:
if self.system._has_system_bath_coupling:
raise ImplementationError("Not yet implemented")
else:
# calculating stick spectra
stick_width = 1.0 / 0.1
def _equilibrium_excit_populations(
self, AG: Any, temperature: float = 300, relaxation_hamiltonian: Any = None
) -> Any:
if relaxation_hamiltonian:
H = relaxation_hamiltonian
else:
H = AG.get_Hamiltonian()
with eigenbasis_of(H):
rho0 = AG.get_DensityMatrix(
condition_type="thermal_excited_state",
relaxation_theory_limit="weak_coupling",
temperature=temperature,
relaxation_hamiltonian=relaxation_hamiltonian,
)
return rho0
[docs]
def one_transition_spectrum_abs(self, tr: dict[str, Any]) -> numpy.ndarray:
"""Calculates spectrum of one transition"""
ta = tr["ta"] # TimeAxis
dd = tr["dd"] # transition dipole strength
om = tr["om"] # frequency - rwa
gg = tr["gg"] # natural broadening (constant or time dependent)
fwhm = tr["fwhm"] # Additional gaussian broadening of the spectra
sgm = fwhm / (2 * numpy.sqrt(2 * numpy.log(2)))
# CD and fluorescence can be calculated in this step
# TODO if rotatory strength defined calculate also circular dichroism spectra
# TOOD calculate fluorescence spectra (for fluorescence there should be a switch because it should be calculated only for the first transition)
if self.system._has_system_bath_coupling:
# ct = tr["ct"] # correlation function
# convert correlation function to lineshape function
gt = tr["gt"]
# gt = _c2g(ta,ct.data)
# calculate time dependent response
at = numpy.exp(-gt - 1j * om * ta.data)
else:
# calculate time dependent response
at = numpy.exp(-1j * om * ta.data)
# plt.figure()
# plt.title("Absorption")
# plt.plot(ta.data,numpy.real(at))
# plt.plot(ta.data,numpy.imag(at))
if len(gg) == 1:
gam = gg[0]
rt = numpy.exp(gam * ta.data)
at *= rt
# print("Constant: ", rt[20], len(at))
else:
rt = numpy.exp((gg) * ta.data)
at *= rt
# print("Time dependent: len = ", rt[20], len(rt))
if fwhm != 0.0:
gauss = numpy.exp(-2 * (numpy.pi**2) * (sgm**2) * (ta.data**2))
at *= gauss
# Fourier transform the result
ft = dd * numpy.fft.hfft(at) * ta.step
ft = numpy.fft.fftshift(ft)
# invert the order because hfft is a transform with -i
ft = numpy.flipud(ft)
# cut the center of the spectrum
Nt = ta.length # len(ta.data)
return ft[Nt // 2 : Nt + Nt // 2]
[docs]
def one_transition_spectrum_ld(self, tr: dict[str, Any]) -> numpy.ndarray:
"""Calculates spectrum of one transition"""
ta = tr["ta"] # TimeAxis
ld = tr["ld"] # linear dichroism strength
om = tr["om"] # frequency - rwa
gg = tr["gg"] # natural broadening (constant or time dependent)
fwhm = tr["fwhm"] # Additional gaussian broadening of the spectra
sgm = fwhm / (2 * numpy.sqrt(2 * numpy.log(2)))
if self.system._has_system_bath_coupling:
# ct = tr["ct"] # correlation function
# convert correlation function to lineshape function
# gt = _c2g(ta,ct.data)
gt = tr["gt"]
# calculate time dependent response
at = numpy.exp(-gt - 1j * om * ta.data)
else:
# calculate time dependent response
at = numpy.exp(-1j * om * ta.data)
if len(gg) == 1:
gam = gg[0]
rt = numpy.exp(gam * ta.data)
at *= rt
# print("Constant: ", rt[20], len(at))
else:
rt = numpy.exp((gg) * ta.data)
at *= rt
# print("Time dependent: len = ", rt[20], len(rt))
if fwhm != 0.0:
gauss = numpy.exp(-2 * (numpy.pi**2) * (sgm**2) * (ta.data**2))
at *= gauss
# Fourier transform the result
ft = ld * numpy.fft.hfft(at) * ta.step
ft = numpy.fft.fftshift(ft)
# invert the order because hfft is a transform with -i
ft = numpy.flipud(ft)
# cut the center of the spectrum
Nt = ta.length # len(ta.data)
return ft[Nt // 2 : Nt + Nt // 2]
[docs]
def one_transition_spectrum_gauss(
self, tr: dict[str, Any]
) -> tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]:
"""Calculates spectrum of one transition using gaussian broadening
of the stick spectra. The definition is the same as for exat tools:
https://doi.org/10.1002/jcc.25118
"""
fa = tr["fa"] # Frequency axis
HWHH = tr["HWHH"] # Half width at the half hight (maximum)
dd = tr["dd"] # transition dipole strength
rr = tr["rr"] # transition dipole strength
ld = tr["ld"] # linear dichroism strength
om = tr["om"] + self.rwa # frequency
# LineShape = lambda p, x: (x/(p[1]*np.sqrt(2*m.pi))*np.exp(-0.5*((x-p[0])/p[1])**2))
# broad = broad/np.sqrt(2*np.log(2))
sigma = HWHH / numpy.sqrt(2 * numpy.log(2))
# x = ta.data
data = (
fa.data
/ (sigma * numpy.sqrt(2 * numpy.pi))
* numpy.exp(-0.5 * ((fa.data - om) / sigma) ** 2)
)
data_abs = dd * data
data_CD = rr * data
data_LD = ld * data
return data_abs, data_CD, data_LD
[docs]
def one_transition_spectrum_fluor(self, tr: dict[str, Any]) -> numpy.ndarray:
"""Calculates spectrum of one transition"""
ta = tr["ta"] # TimeAxis
dd = tr["dd"] # transition dipole strength
om = tr["om"] # frequency - rwa
gg = tr["gg"] # natural broadening (constant or time dependent)
fwhm = tr["fwhm"] # Additional gaussian broadening of the spectra
sgm = fwhm / (2 * numpy.sqrt(2 * numpy.log(2)))
# CD and fluorescence can be calculated in this step
# TODO if rotatory strength defined calculate also circular dichroism spectra
# TOOD calculate fluorescence spectra (for fluorescence there should be a switch because it should be calculated only for the first transition)
if self.system._has_system_bath_coupling:
# ct = tr["ct"] # correlation function
re = tr["re"] # reorganisation energy
# convert correlation function to lineshape function
# gt = _c2g(ta,ct.data)
gt = tr["gt"]
# calculate time dependent response
at = numpy.exp(-numpy.conjugate(gt) - 1j * om * ta.data + 2j * re * ta.data)
else:
# calculate time dependent response
at = numpy.exp(-1j * om * ta.data)
# plt.figure()
# plt.title("Absorption")
# plt.plot(ta.data,numpy.real(at))
# plt.plot(ta.data,numpy.imag(at))
if len(gg) == 1:
gam = gg[0]
rt = numpy.exp(gam * ta.data)
at *= rt
# print("Constant: ", rt[20], len(at))
else:
rt = numpy.exp((gg) * ta.data)
at *= rt
# print("Time dependent: len = ", rt[20], len(rt))
if fwhm != 0.0:
gauss = numpy.exp(-2 * (numpy.pi**2) * (sgm**2) * (ta.data**2))
at *= gauss
# Fourier transform the result
ft = dd * numpy.fft.hfft(at) * ta.step
ft = numpy.fft.fftshift(ft)
# invert the order because hfft is a transform with -i
ft = numpy.flipud(ft)
# cut the center of the spectrum
Nt = ta.length # len(ta.data)
return ft[Nt // 2 : Nt + Nt // 2]
[docs]
def one_transition_spectrum_cd(self, tr: dict[str, Any]) -> numpy.ndarray:
"""Calculates spectrum of one transition"""
ta = tr["ta"] # TimeAxis
rr = tr["rr"] # transition dipole strength
om = tr["om"] # frequency - rwa
gg = tr["gg"] # natural broadening (constant or time dependent)
fwhm = tr["fwhm"] # Additional gaussian broadening of the spectra
sgm = fwhm / (2 * numpy.sqrt(2 * numpy.log(2)))
# CD and fluorescence can be calculated in this step
# TODO if rotatory strength defined calculate also circular dichroism spectra
# TOOD calculate fluorescence spectra (for fluorescence there should be a switch because it should be calculated only for the first transition)
if self.system._has_system_bath_coupling:
# ct = tr["ct"] # correlation function
# convert correlation function to lineshape function
# gt = _c2g(ta,ct.data)
gt = tr["gt"]
# calculate time dependent response
at = numpy.exp(-gt - 1j * om * ta.data)
else:
# calculate time dependent response
at = numpy.exp(-1j * om * ta.data)
# plt.figure()
# plt.title("Absorption")
# plt.plot(ta.data,numpy.real(at))
# plt.plot(ta.data,numpy.imag(at))
if len(gg) == 1:
gam = gg[0]
rt = numpy.exp(gam * ta.data)
at *= rt
# print("Constant: ", rt[20], len(at))
else:
rt = numpy.exp((gg) * ta.data)
at *= rt
# print("Time dependent: len = ", rt[20], len(rt))
if fwhm != 0.0:
gauss = numpy.exp(-2 * (numpy.pi**2) * (sgm**2) * (ta.data**2))
at *= gauss
# Fourier transform the result
ft = rr * numpy.fft.hfft(at) * ta.step
ft = numpy.fft.fftshift(ft)
# invert the order because hfft is a transform with -i
ft = numpy.flipud(ft)
# cut the center of the spectrum
Nt = ta.length # len(ta.data)
return ft[Nt // 2 : Nt + Nt // 2]
def _excitonic_coft_old(self, SS: numpy.ndarray, AG: Any, n: int) -> numpy.ndarray:
"""Returns energy gap correlation function data of an exciton state"""
# FIXME: works only for 2 level molecules
c0 = AG.monomers[0].get_transition_environment((0, 1)).data # get_egcf((0,1))
Nt = len(c0)
# SystemBathInteraction
sbi = AG.get_SystemBathInteraction()
# CorrelationFunctionMatrix
cfm = sbi.CC
# get number of monomeric basis states
Na = 0
for monomer in AG.monomers:
Na += monomer.nel - 1
ct = numpy.zeros((Nt), dtype=numpy.complex128)
# Na = AG.nmono
for kk in range(Na):
# nkk = AG.monomers[kk].egcf_mapping[0]
for ll in range(Na):
# nll = AG.monomers[ll].egcf_mapping[0]
ct += (
(SS[kk + 1, n + 1] ** 2)
* (SS[ll + 1, n + 1] ** 2)
* cfm.get_coft(kk, ll)
)
# *AG.egcf_matrix.get_coft(nkk,nll))
return ct
def _excitonic_coft(self, SS: numpy.ndarray, AG: Any, n: int) -> numpy.ndarray:
"""Returns energy gap correlation function data of an exciton state n"""
# SystemBathInteraction
sbi = AG.get_SystemBathInteraction()
# CorrelationFunctionMatrix
cfm = sbi.CC
c0 = AG.monomers[0].get_egcf((0, 1))
Nt = len(c0)
ct = numpy.zeros((Nt), dtype=numpy.complex128)
# electronic states corresponding to single excited states
elst = numpy.where(AG.which_band == 1)[0]
for el1 in elst:
for el2 in elst:
if cfm.cpointer[el1 - 1, el2 - 1] == 0:
continue
coft = cfm.get_coft(el1 - 1, el2 - 1)
for kk in AG.vibindices[el1]:
for ll in AG.vibindices[el2]:
ct += (SS[kk, n] ** 2) * (SS[ll, n] ** 2) * coft
return ct
def _excitonic_coft_all(self, SS: numpy.ndarray, AG: Any) -> numpy.ndarray:
"""Returns energy gap correlation function data of an exciton state n"""
# SystemBathInteraction
sbi = AG.get_SystemBathInteraction()
# CorrelationFunctionMatrix
cfm = sbi.CC
c0 = AG.monomers[0].get_egcf((0, 1))
Nt = len(c0)
Nst = AG.HamOp.dim
ct = numpy.zeros((Nst, Nt), dtype=numpy.complex128)
# electronic states corresponding to single excited states
import time
timecount: float = 0.0
elst = numpy.where(AG.which_band == 1)[0]
start = time.time()
for el1 in elst:
for el2 in elst:
coft = cfm.get_coft(el1 - 1, el2 - 1)
start2 = time.time()
for kk in AG.vibindices[el1]:
for ll in AG.vibindices[el2]:
ct[:, :] += numpy.dot(
numpy.expand_dims(
(SS[kk, :] ** 2) * (SS[ll, :] ** 2), axis=1
),
numpy.expand_dims(coft, axis=0),
)
stop2 = time.time()
timecount += stop2 - start2
stop = time.time()
print(stop - start, stop - start - timecount)
return ct
def _excitonic_reorg_energy(self, SS: numpy.ndarray, AG: Any, n: int) -> float:
"""Returns the reorganisation energy of an exciton state"""
# SystemBathInteraction
sbi = AG.get_SystemBathInteraction()
# CorrelationFunctionMatrix
cfm = sbi.CC
rg = 0.0
# electronic states corresponding to single excited states
elst = numpy.where(AG.which_band == 1)[0]
for el1 in elst:
reorg = cfm.get_reorganization_energy(el1 - 1, el1 - 1)
for kk in AG.vibindices[el1]:
rg += (SS[kk, n] ** 2) * (SS[kk, n] ** 2) * reorg
return rg
def _excitonic_rotatory_strength(
self, SS: numpy.ndarray, AG: Any, energy: numpy.ndarray, n: int
) -> float:
# Initialize rotatory strength
Rot_n = 0
for ii in range(AG.Ntot):
for jj in range(ii + 1, AG.Ntot):
Rot_n += SS[ii, n] * SS[jj, n] * AG.RR[ii, jj]
# print(ii,jj,SS[ii,n],SS[jj,n],AG.RR[ii,jj])
#
#
# # electronic states corresponding to single excited states
# elst = numpy.where(AG.which_band == 1)[0]
# for el1 in elst:
# # get monomer number
# mon_indx = numpy.nonzero(AG.elsigs[el1])[0][0]
# mon1 = AG.monomers[mon_indx]
# Ri = mon1.position
# for el2 in elst:
# if el2>el1:
# mon_indx = numpy.nonzero(AG.elsigs[el2])[0][0]
# mon2 = AG.monomers[mon_indx]
# Rj = mon2.position
#
# for ii in AG.vibindices[el1]:
# di = AG.vibdipoles[0,ii]
# for jj in AG.vibindices[el2]:
# dj = AG.vibdipoles[0,jj]
# Scale by excitation energy:
Rot_n *= energy[n]
# print(Rot_n,energy,n)
return Rot_n
def _excitonic_rotatory_strength_fullv(
self, SS: numpy.ndarray, AG: Any, energy: numpy.ndarray, n: int
) -> Any:
# Initialize rotatory strength
Rot_n = 0
# Rot_nm = 0
#
# DD_vel = self.system.DD[0].copy()
# for ii in range(1:DD_vel.shape[0]):
# DD_vel[ii] *= -self.system.HH[ii,ii]
#
# print("Rot1:",numpy.dot(SS[:,n],numpy.dot(AG.RRm,SS[:,n]))/energy[n])
# print("Rot2:",numpy.dot(SS[:,n],numpy.dot(AG.RRv,SS[:,n])))
if AG._has_velocity_dipoles:
Rot_n = numpy.dot(SS[:, n], numpy.dot(AG.RRm, SS[:, n])) / energy[n]
else:
# Rot_n = energy[n]*numpy.dot(SS[:,n],numpy.dot(AG.RR+AG.RRm,SS[:,n]))
Rot_n = numpy.dot(SS[:, n], numpy.dot(AG.RRv, SS[:, n]))
# Rot_n = energy[n]*numpy.dot(SS[:,n],numpy.dot(AG.RR+AG.RRm,SS[:,n]))
# Rot_n = numpy.dot(SS[:,n],numpy.dot(AG.RRm,SS[:,n]))/energy[n]
# for ii in range(AG.Ntot):
# for jj in range(AG.Ntot):
# Rot_n += SS[ii,n]*SS[jj,n]*(AG.RRv[ii,jj]+AG.RRm[ii,jj])
#
# for ii in range(AG.Ntot):
# for jj in range(AG.Ntot):
# Rot_nm += SS[ii,n]*SS[jj,n]*AG.RRm[ii,jj]
return Rot_n # ,Rot_nm
def _subtract_site_reorg(
self, AG: Any, Hin: Any, subtract_bath: bool = True
) -> tuple[Any, numpy.ndarray]:
# Subtract the reorganisation energy from the sites
sbi = AG.get_SystemBathInteraction()
# CorrelationFunctionMatrix
cfm = sbi.CC
reorg_site = numpy.zeros(Hin.dim)
HH = copy.deepcopy(Hin)
# electronic states corresponding to single excited states
elst = numpy.where(AG.which_band == 1)[0]
for el1 in elst:
reorg = cfm.get_reorganization_energy(el1 - 1, el1 - 1)
reorg_bath = 0.0
if subtract_bath:
coft = cfm.cfuncs[cfm.get_index_by_where((el1 - 1, el1 - 1))]
for parm in coft.params:
if parm["ftype"] == "OverdampedBrownian":
reorg_bath += parm["reorg"]
for kk in AG.vibindices[el1]:
reorg_site[kk] = reorg - reorg_bath
HH._data[kk, kk] -= reorg_site[kk]
return HH, reorg_site
def _site2excit_reorg(
self, reorg_site: numpy.ndarray, SS: numpy.ndarray
) -> numpy.ndarray:
# Get exciton reorganization energy from the site one
reorg_exct = numpy.zeros(SS.shape[1])
for n in range(SS.shape[1]):
reorg_exct[n] = numpy.dot(SS[:, n] ** 4, reorg_site)
return reorg_exct
def _calculate_monomer(self, raw: bool = False) -> Any:
"""Calculates the absorption spectrum of a monomer"""
ta = self.TimeAxis
# transition frequency
om = self.system.elenergies[1] - self.system.elenergies[0]
# transition dipole moment
dm = self.system.dmoments[0, 1, :]
# dipole^2
dd = numpy.dot(dm, dm)
# natural life-time from the dipole moment
gama = [0.0] # [-1.0/self.system.get_electronic_natural_lifetime(1)]
sbi = self.system.get_SystemBathInteraction()
reorg = sbi.CC.get_reorganization_energy(0, 0)
if self.system._has_system_bath_coupling:
# correlation function
ct = self.system.get_egcf((0, 1))
gt = _c2g(ta, ct.data)
tr = {
"ta": ta,
"dd": dd,
"om": om - self.rwa,
"ct": ct,
"gt": gt,
"gg": gama,
"fwhm": 0.0,
}
else:
tr = {"ta": ta, "dd": dd, "om": om - self.rwa, "gg": gama, "fwhm": 0.0}
if self._gauss_broad:
tr["fwhm"] = self.gauss
tr["re"] = reorg
if self._gauss_broad:
tr["fwhm"] = self.gauss
# calculates the one transition of the monomer
data = numpy.real(self.one_transition_spectrum_abs(tr))
data_fl = numpy.real(self.one_transition_spectrum_fluor(tr))
for ii in range(2, self.system.Nb[1] + 1):
# transition frequency
om = self.system.elenergies[ii] - self.system.elenergies[0]
# transition dipole moment
dm = self.system.dmoments[0, ii, :]
# dipole^2
dd = numpy.dot(dm, dm)
# natural life-time from the dipole moment
gama = [0.0] # [-1.0/self.system.get_electronic_natural_lifetime(ii)]
if self.system._has_system_bath_coupling:
# correlation function
ct = self.system.get_egcf((0, ii))
gt = _c2g(ta, ct.data)
tr = {
"ta": ta,
"dd": dd,
"om": om - self.rwa,
"ct": ct,
"gt": gt,
"gg": gama,
"fwhm": 0.0,
}
else:
tr = {"ta": ta, "dd": dd, "om": om - self.rwa, "gg": gama, "fwhm": 0.0}
if self._gauss_broad:
tr["fwhm"] = self.gauss
if self._gauss_broad:
tr["fwhm"] = self.gauss
data += numpy.real(self.one_transition_spectrum_abs(tr))
# we only want to retain the upper half of the spectrum
Nt = len(self.frequencyAxis.data) // 2
do = self.frequencyAxis.data[1] - self.frequencyAxis.data[0]
st = self.frequencyAxis.data[Nt // 2]
# we represent the Frequency axis anew
axis = FrequencyAxis(st, Nt, do)
# multiply the spectrum by frequency (compulsory prefactor)
if not raw:
data = axis.data * data
data_fl = (axis.data**3) * data_fl
spect = AbsSpectrum(axis=axis, data=data)
return spect
def _calculate_abs_from_dynamics(self, raw: bool = False, alt: bool = False) -> Any:
"""Calculates the absorption spectrum of a molecule from its dynamics"""
#
# Frequency axis
#
# we only want to retain the upper half of the spectrum
Nt = len(self.frequencyAxis.data) // 2
do = self.frequencyAxis.data[1] - self.frequencyAxis.data[0]
st = self.frequencyAxis.data[Nt // 2]
# we represent the Frequency axis anew
axis = FrequencyAxis(st, Nt, do)
# the spectrum will be calculate for this system
system = self.system
HH = system.get_Hamiltonian()
if not HH.has_rwa:
raise QuantarheiError(
"Hamiltonian has to define Rotating Wave Approximation"
)
# with energy_units("1/cm"):
# print("Skeleton")
# print(numpy.diag(HH.get_RWA_skeleton()))
# print(HH.get_RWA_data())
# transition dipole moment
DD = system.get_TransitionDipoleMoment()
rhoeq = system.get_thermal_ReducedDensityMatrix()
# the propagator to propagate optical coherences
prop = self.prop
# FIXME: time axis must be consistent with other usage of the calculator
time = prop.TimeAxis
secular = False
if alt:
at = _spect_from_dyn(time, HH, DD, prop, rhoeq, secular)
else:
at = _spect_from_dyn_single(time, HH, DD, prop, rhoeq, secular)
#
# Fourier transform of the time-dependent result
#
ft = numpy.fft.hfft(at) * time.step
ft = numpy.fft.fftshift(ft)
# invert the order because hfft is a transform with -i
ft = numpy.flipud(ft)
# cut the center of the spectrum
Nt = time.length # len(ta.data)
data = ft[Nt // 2 : Nt + Nt // 2]
# FIXME: Fluorescence spectrum from dynamics not implemented
data_fl = numpy.zeros_like(data)
#
# multiply the response by \omega
#
if not raw:
data = axis.data * data
spect_abs = LinSpectrum(axis=axis, data=data)
fluor_spect = LinSpectrum(axis=axis, data=data_fl)
return {"abs": spect_abs, "fluor": fluor_spect}
# return spect
def _calculate_aggregate(
self,
relaxation_tensor: Any = None,
relaxation_hamiltonian: Any = None,
rate_matrix: Any = None,
raw: bool = False,
) -> dict[str, Any]:
"""Calculates the absorption spectrum of a molecular aggregate"""
ta = self.TimeAxis
# Hamiltonian of the system
if relaxation_hamiltonian is None:
HH = self.system.get_Hamiltonian()
else:
HH = relaxation_hamiltonian
if self._is_adiabatic:
# subtract site reorganisation energy from the sites
HH, reorg_site = self._subtract_site_reorg(
self.system, HH, subtract_bath=self._adiabatic_noBath
)
SS = HH.diagonalize() # transformed into eigenbasis
energy = numpy.diag(HH.data)
if self._is_adiabatic:
# get exciton reorganization energy
reorg_excit = self._site2excit_reorg(reorg_site, SS)
# shift the diagonal of the exciton hamiltonian
for ii in range(HH.dim):
HH._data[ii, ii] += reorg_excit[ii]
# SS = HH.diagonalize() # transformed into eigenbasis
# Transition dipole moment operator
DD = self.system.get_TransitionDipoleMoment()
# transformed into the basis of Hamiltonian eigenstates
DD.transform(SS)
# Frequency axis
# we only want to retain the upper half of the spectrum
Nt = len(self.frequencyAxis.data) // 2
do = self.frequencyAxis.data[1] - self.frequencyAxis.data[0]
st = self.frequencyAxis.data[Nt // 2]
# we represent the Frequency axis anew
axis = FrequencyAxis(st, Nt, do)
# TimeAxis
tr = {"ta": ta, "fa": axis, "fwhm": 0}
# If gaussian groadening is specified calculate also spectra
if self._gass_lineshape:
tr["HWHH"] = self.HWHH
if self._gauss_broad:
tr["fwhm"] = self.gauss
if relaxation_tensor is not None:
RR: Any = relaxation_tensor
RR.transform(SS)
gg = []
RR_any: Any = RR
if isinstance(RR, TimeDependent):
for ii in range(HH.dim):
gg.append(RR_any.data[:, ii, ii, ii, ii])
else:
for ii in range(HH.dim):
gg.append([RR_any.data[ii, ii, ii, ii]])
tr["gg"] = gg[1]
elif rate_matrix is not None:
RR = rate_matrix # rate matrix is in excitonic basis
gg = []
RR_any = RR
if isinstance(RR, TimeDependent):
for ii in range(HH.dim):
gg.append(RR_any.data[:, ii, ii] / 2.0)
else:
for ii in range(HH.dim):
gg.append([RR_any.data[ii, ii] / 2.0])
tr["gg"] = gg[1]
else:
tr["gg"] = [0.0]
# get square of transition dipole moment here
tr["dd"] = DD.dipole_strength(0, 1)
# first transition energy
tr["om"] = HH.data[1, 1] - HH.data[0, 0] - self.rwa
# get a transformed ct here
ct = self._excitonic_coft(SS, self.system, 1)
tr["ct"] = ct
# calculate g(t)
tr["gt"] = _c2g(tr["ta"], tr["ct"].data)
# get reorganization energy
tr["re"] = self._excitonic_reorg_energy(SS, self.system, 1)
# get rotatory strength
tr["rr"] = self._excitonic_rotatory_strength_fullv(SS, self.system, energy, 1)
# tr["rr"] = self._excitonic_rotatory_strength(SS,self.system,energy,1)
# print(1,convert(tr["rr"],"int","1/cm")*numpy.pi*1e-4)
dip = DD.data[0, 1]
tr["ld"] = 3 * (numpy.dot(dip, self.ld_axis)) ** 2 - tr["dd"] # *3/2
# tr["ld"] = 3*(DD. # *3/2
self.system._has_system_bath_coupling = True
temperature = self.system.sbi.get_temperature()
if self._pop_fluor == "vertical":
rho_eq_exct = self._equilibrium_excit_populations(
self.system, temperature=temperature
)
elif self._pop_fluor == "adiabatic":
eng = numpy.zeros(HH.dim - 1)
rho_eq_exct = numpy.zeros(HH.dim)
for ii in range(1, HH.dim):
reorg = self._excitonic_reorg_energy(SS, self.system, ii)
eng[ii - 1] = HH.data[ii, ii] - HH.data[0, 0] - self.rwa - reorg
# now get equlibrium population on energy levels
for ii in range(1, min(14, HH.dim)):
rho_eq_exct[ii] = numpy.exp(-eng[ii - 1] / kB_intK / temperature)
rho_eq_exct /= numpy.sum(rho_eq_exct)
rho_eq_exct = numpy.diag(rho_eq_exct)
else:
raise Warning(
"Unknown method for computing excited state equilibrium populations"
)
# print(tr["ct"])
# print(max(tr["ct"]))
# print(self.one_transition_spectrum(tr))
# print(max(self.one_transition_spectrum(tr)))
#
# Calculates spectrum of a single transition
#
data = numpy.real(self.one_transition_spectrum_abs(tr))
data_fl = numpy.real(
rho_eq_exct.data[1, 1] * self.one_transition_spectrum_fluor(tr)
)
data_cd = numpy.real(self.one_transition_spectrum_cd(tr))
data_ld = numpy.real(self.one_transition_spectrum_ld(tr))
if self._gass_lineshape:
data_gauss, data_cd_gauss, data_ld_gauss = numpy.real(
self.one_transition_spectrum_gauss(tr)
)
# print("Population mine:",rho_eq_exct.data[1, 1])
# FOR THE VIBRONIC SYSTEM THE SPECTRA HAVE TO BE SUMED THROUGH THE GROUND STATES (VIBRONIC)
for jj in range(1, min(1, self.system.Nb[0])): # sum over the ground states
tr["dd"] = DD.dipole_strength(jj, 1)
# first transition energy
tr["om"] = HH.data[1, 1] - HH.data[jj, jj] - self.rwa
data_fl += numpy.real(
rho_eq_exct.data[1, 1] * self.one_transition_spectrum_fluor(tr)
)
for ii in range(2, HH.dim):
if relaxation_tensor is not None or rate_matrix is not None:
tr["gg"] = gg[ii]
else:
tr["gg"] = [0.0]
# print(tr["gg"])
# tr[1] = DD.dipole_strength(0,ii) # update transition dipole moment
tr["dd"] = DD.dipole_strength(0, ii)
# tr[2] = HH.data[ii,ii]-HH.data[0,0]-rwa
tr["om"] = HH.data[ii, ii] - HH.data[0, 0] - self.rwa
# tr[3] = self._excitonic_coft(SS,self.system,ii-1) # update ct here
tr["ct"] = self._excitonic_coft(SS, self.system, ii)
tr["gt"] = _c2g(tr["ta"], tr["ct"].data)
tr["re"] = self._excitonic_reorg_energy(SS, self.system, ii)
tr["rr"] = self._excitonic_rotatory_strength_fullv(
SS, self.system, energy, ii
)
dip = DD.data[0, ii]
tr["ld"] = 3 * (numpy.dot(dip, self.ld_axis)) ** 2 - tr["dd"] # *3/2
# tr["rr"] = self._excitonic_rotatory_strength(SS,self.system,energy,ii)
# print(ii,convert(HH.data[ii,ii]-HH.data[0,0]-tr["re"],"int","1/cm"),convert(HH.data[ii,ii]-HH.data[0,0]-2*tr["re"],"int","1/cm"),convert(tr["re"],"int","1/cm"))
# print(ii,convert(tr["rr"],"int","1/cm")*numpy.pi*1e-4)
# conversion factor is convert rotatory strength to inverse centimeters and multiply *numpy.pi*1e-4
#
# Calculates spectrum of a single transition
#
data += numpy.real(self.one_transition_spectrum_abs(tr))
data_fl += numpy.real(
rho_eq_exct.data[ii, ii] * self.one_transition_spectrum_fluor(tr)
)
data_cd += numpy.real(self.one_transition_spectrum_cd(tr))
data_ld += numpy.real(self.one_transition_spectrum_ld(tr))
if self._gass_lineshape:
data_gauss_tmp, data_cd_gauss_tmp, data_ld_gauss_tmp = numpy.real(
self.one_transition_spectrum_gauss(tr)
)
data_gauss += data_gauss_tmp
data_cd_gauss += data_cd_gauss_tmp
data_ld_gauss += data_ld_gauss_tmp
# print("Population mine:",rho_eq_exct.data[ii, ii])
# FOR THE VIBRONIC SYSTEM THE SPECTRA HAVE TO BE SUMED THROUGH THE GROUND STATES (VIBRONIC)
for jj in range(
1, min(ii, self.system.Nb[0])
): # sum over the ground states
tr["dd"] = DD.dipole_strength(jj, ii)
# first transition energy
tr["om"] = HH.data[ii, ii] - HH.data[jj, jj] - self.rwa
data_fl += numpy.real(
rho_eq_exct.data[ii, ii] * self.one_transition_spectrum_fluor(tr)
)
# multiply the spectrum by frequency (compulsory prefactor)
if not raw:
data = axis.data * data
data_fl = axis.data * data_fl
data_cd = axis.data * data_cd
data_ld = axis.data * data_ld
# data_gauss = axis.data*data_gauss
# data_cd_gauss = axis.data*data_cd_gauss
# transform all quantities back
S1 = numpy.linalg.inv(SS)
HH.transform(S1)
DD.transform(S1)
if relaxation_tensor is not None:
RR.transform(S1)
abs_spect = LinSpectrum(axis=axis, data=data)
fluor_spect = LinSpectrum(axis=axis, data=data_fl)
CD_spect = LinSpectrum(axis=axis, data=data_cd)
LD_spect = LinSpectrum(axis=axis, data=data_ld)
if self._gass_lineshape:
abs_spect_gauss = LinSpectrum(axis=axis, data=data_gauss)
CD_spect_gauss = LinSpectrum(axis=axis, data=data_cd_gauss)
LD_spect_gauss = LinSpectrum(axis=axis, data=data_ld_gauss)
return {
"abs": abs_spect,
"fluor": fluor_spect,
"CD": CD_spect,
"LD": LD_spect,
"LD gauss": LD_spect_gauss,
"abs gauss": abs_spect_gauss,
"CD gauss": CD_spect_gauss,
}
return {"abs": abs_spect, "fluor": fluor_spect, "CD": CD_spect, "LD": LD_spect}
[docs]
class AbsSpectrumCalculator(LinSpectrumCalculator):
def __init__(
self,
timeaxis: Any,
system: Any = None,
dynamics: str = "secular",
relaxation_tensor: Any = None,
rate_matrix: Any = None,
effective_hamiltonian: Any = None,
) -> None:
super().__init__(
timeaxis,
system=system,
dynamics=dynamics,
relaxation_tensor=relaxation_tensor,
rate_matrix=rate_matrix,
effective_hamiltonian=effective_hamiltonian,
)
[docs]
@prevent_basis_context
def calculate(
self, raw: bool = False, from_dynamics: bool = False, alt: bool = False
) -> Any:
if not self.bootstrapped:
raise QuantarheiError(
"Calculator must be bootstrapped first: "
"call bootstrap() method of this object."
)
spect: Any = None
with energy_units("int"):
if self.system is not None:
if from_dynamics:
# alt = True is for testing only
spect = self._calculate_abs_from_dynamics(raw=raw, alt=alt)
if isinstance(self.system, Molecule):
# self._calculate_Molecule(rwa)
spect = self._calculate_monomer(raw=raw)
elif isinstance(self.system, Aggregate):
spect = self._calculate_aggregate(
relaxation_tensor=self._relaxation_tensor,
rate_matrix=self._rate_matrix,
relaxation_hamiltonian=self._relaxation_hamiltonian,
raw=raw,
)["abs"]
else:
raise QuantarheiError("System to calculate spectrum for not defined")
return spect
def _spect_from_dyn(
time: Any, HH: Any, DD: Any, prop: Any, rhoeq: Any, secular: bool = False
) -> numpy.ndarray:
"""Calculation of the first order signal field.
Parameters
----------
time : TimeAxis
Time axis on which the spectrum is calculated
HH : Hamiltonian
Hamiltonian of the system
DD : TransitionDipoleMoment
Transition dipole moment operator
prop : ReducedDensityMatrixPropagator
Propagator of the density matrix
rhoeq : ReducedDensityMatrix
Equilibrium reduced density matrix
secular : bool, optional
Should secular approximation be used? The default is False.
Returns
-------
at : numpy.array, COMPLEX
Time dependent signal field.
"""
# we will loop over transitions
rhoi = ReducedDensityMatrix(dim=HH.dim)
#
# Time dependent data
#
at = numpy.zeros(time.length, dtype=COMPLEX)
# transitions to loop over
# ig represents all states in the ground state block
for ig in range(HH.rwa_indices[1]):
if len(HH.rwa_indices) <= 2:
Nfin = HH.dim
else:
Nfin = HH.rwa_indices[2]
for ie in range(HH.rwa_indices[1], Nfin):
rhoi.data[:, :] = 0.0
rhoi.data[ie, ig] = rhoeq.data[ig, ig]
dvec_1 = DD.data[ie, ig, :]
dd = numpy.dot(dvec_1, dvec_1)
rhot = prop.propagate(rhoi)
if secular:
at += (dd / 3.0) * rhot.data[:, ie, ig]
else:
#
# non-secular loops over all coherences
#
for jg in range(HH.rwa_indices[1]):
for je in range(HH.rwa_indices[1], Nfin):
dvec_2 = DD.data[je, jg, :]
d12 = numpy.dot(dvec_2, dvec_1)
at += (d12 / 3.0) * rhot.data[:, je, jg]
return at
def _spect_from_dyn_single(
time: Any, HH: Any, DD: Any, prop: Any, rhoeq: Any, secular: bool = False
) -> numpy.ndarray:
"""Calculation of the first order signal field.
Calculation in a single propagation
Parameters
----------
time : TimeAxis
Time axis on which the spectrum is calculated
HH : Hamiltonian
Hamiltonian of the system
DD : TransitionDipoleMoment
Transition dipole moment operator
prop : ReducedDensityMatrixPropagator
Propagator of the density matrix
rhoeq : ReducedDensityMatrix
Equilibrium reduced density matrix
secular : bool, optional
Should secular approximation be used? The default is False.
Returns
-------
at : numpy.array, COMPLEX
Time dependent signal field.
"""
rhoi = ReducedDensityMatrix(dim=HH.dim)
#
# Time dependent data
#
at = numpy.zeros(time.length, dtype=COMPLEX)
# deff = numpy.sqrt(numpy.einsum("ijn,ijn->ij", DD.data,DD.data,
# dtype=REAL))
_show = False
if _show:
import matplotlib.pyplot as plt
for kk in range(3):
deff = DD.data[:, :, kk]
# excitation by an effective dipole
rhoi.data = (numpy.dot(deff, rhoeq.data) + numpy.dot(rhoeq.data, deff)) / 3.0
rhot = prop.propagate(rhoi)
if _show:
for ig in range(HH.rwa_indices[1]):
print("ig=", ig)
for ll in range(rhot.data.shape[2]):
plt.plot(time.data, rhot.data[:, ig, ll])
plt.title("kk=" + str(kk))
plt.show()
for ig in range(HH.rwa_indices[1]):
at += numpy.einsum("j,kj->k", deff[ig, :], rhot.data[:, :, ig])
return at
def _c2g(timeaxis: Any, 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 : TimeAxis
TimeAxis of the correlation function
coft : complex numpy array
Values of correlation function given at points specified
in the TimeAxis object
"""
ta = timeaxis
rr = numpy.real(coft)
ri = numpy.imag(coft)
sr = scipy.interpolate.UnivariateSpline(ta.data, rr, s=0).antiderivative()(ta.data)
sr = scipy.interpolate.UnivariateSpline(ta.data, sr, s=0).antiderivative()(ta.data)
si = scipy.interpolate.UnivariateSpline(ta.data, ri, s=0).antiderivative()(ta.data)
si = scipy.interpolate.UnivariateSpline(ta.data, si, s=0).antiderivative()(ta.data)
gt = sr + 1j * si
return gt