"""Class representing aggregates of molecules.
The class enables building of complicated objects from objects of the Molecule
type, their mutual interactions and system-bath interactions. It also provides
an interface to various methods of open quantum systems theory.
Issues:
- Only ground to excited state transition widths and dephasing ( no 1->2 band transition widths)
- No energy gap correlation function for the 1->2 band transitions (only ground to excited ones)
* Coft is for the state and not transition for the 1->2 transition on single molecule the cofts
has to be properly subtracted
- Transformation of the transition width for multilevel molecules (check), transform the same
for the vibronic agregate for the multilevel molecule
"""
from __future__ import annotations
import warnings
from collections.abc import Generator, Iterator
from typing import Any
import numpy
from .. import COMPLEX, REAL
from ..core.managers import Manager, UnitsManaged, eigenbasis_of
from ..core.saveable import Saveable
from ..exceptions import BuildError, QuantarheiError
from ..qm.corfunctions import CorrelationFunctionMatrix
from ..qm.hilbertspace.dmoment import TransitionDipoleMoment
from ..qm.hilbertspace.hamiltonian import Hamiltonian
from ..qm.hilbertspace.operators import DensityMatrix, Operator, ReducedDensityMatrix
from ..qm.hilbertspace.statevector import StateVector
from ..qm.liouvillespace.systembathinteraction import SystemBathInteraction
from ..qm.oscillators.ho import fcstorage, operator_factory
from ..qm.propagators.dmevolution import (
DensityMatrixEvolution,
ReducedDensityMatrixEvolution,
)
from ..qm.propagators.statevectorevolution import StateVectorEvolution
from .aggregate_states import ElectronicState, VibronicState
from .interactions import dipole_dipole_interaction
from .opensystem import OpenSystem
[docs]
class AggregateBase(UnitsManaged, Saveable, OpenSystem):
"""Molecular aggregate
Parameters
----------
name : str
Specifies the name of the aggregate
molecules : list or tuple
List of molecules out of which the aggregate is built
"""
def __init__(self, molecules: list | None = None, name: str = "") -> None:
# OpenSystem.__init__(self)
self.mnames: dict[str, int] = {} #
self.monomers: list[Any] = []
self.nmono = 0 #
self.name = name #
self.mult = 0 #
self.lab: Any = None #
self._has_egcf_matrix = False #
self.egcf_matrix: Any = None
self._has_system_bath_interaction = False #
self._has_lindich_axes = False
self.coupling_initiated = False #
self.resonance_coupling: Any = None
self._has_velocity_dipoles = False
if molecules is not None:
for m in molecules:
self.add_Molecule(m)
self._init_me()
def _init_me(self) -> None:
"""Initializes all built attributes of the aggregate
This should put the object into a pre-build state
"""
self.FC = fcstorage()
self.ops = operator_factory()
self._has_relaxation_tensor = False #
self._relaxation_theory = "" #
# self._built = False #
self._diagonalized = False
self.mult = 0 #
self.sbi_mult = 0 #
self.Nel = 0 #
self.Ntot = 0 #
self.Nb: Any = 0 #
self.vibindices: list[Any] = []
self.which_band: Any = None
self.elsigs: Any = None
self.HH: Any = None
self.HamOp: Any = None
self.DD: Any = None
self.MM: Any = None
self.RR: Any = None
self.Wd: Any = None
self.Dr: Any = None
self.D2: Any = None
self.D2_max: float = 0.0
self.sbi: Any = None
self._has_wpm = False
self.WPM: Any = None
self.sv0: Any = None
self.vibsigs: Any = None
[docs]
def clean(self) -> None:
"""Cleans the aggregate object of anything built
This operation leaves the molecules of the aggregate intact and keeps
few more pieces of information it it. E. g. coupling matrix is not
deleted. You call build again after this.
"""
self._init_me()
[docs]
def wipe_out(self) -> None:
"""Removes everything except of name attribute
You have to set molecules and recalculate interactions before you can
build
"""
self.mnames = {} #
self.monomers = []
self.nmono = 0 #
self.mult = 0 #
self.sbi_mult = 0 #
self._has_egcf_matrix = False #
self.egcf_matrix = None
self._has_system_bath_interaction = False #
self.coupling_initiated = False #
self.resonance_coupling = None #
self._init_me()
########################################################################
#
# BUILDING METHODS
#
########################################################################
[docs]
def init_coupling_indexes(self) -> None:
"""Set indexes for elements of the coupling matrix
index for coupling between mon1 init1 -> final1 and mon2 init2 -> final2
is self._mol2coupling[mon1][init1][final1-init1-1][mon2-mon1-1][init2][final2-init2-1]
"""
# FIXME: Find better way how to store couplings for multilevel molecules!
self._coupling2mol: list[Any] = []
self._mol2coupling: list[Any] = []
# index for coupling between mon1 init1 -> final1 and mon2 init2 -> final2
# is self._mol2coupling[mon1][init1][final1-init1-1][mon2-mon1-1][init2][final2-init2-1]
count = 0
for mon1 in range(self.nmono - 1):
monomer1 = self.monomers[mon1]
self._mol2coupling.append([])
for init1 in range(monomer1.nel - 1):
self._mol2coupling[mon1].append([])
for fin1 in range(init1 + 1, monomer1.nel):
self._mol2coupling[mon1][init1].append([])
for mon2 in range(mon1 + 1, self.nmono):
self._mol2coupling[mon1][init1][fin1 - init1 - 1].append([])
monomer2 = self.monomers[mon2]
for init2 in range(monomer2.nel - 1):
self._mol2coupling[mon1][init1][fin1 - init1 - 1][
mon2 - mon1 - 1
].append([])
for fin2 in range(init2 + 1, monomer2.nel):
self._mol2coupling[mon1][init1][fin1 - init1 - 1][
mon2 - mon1 - 1
][init2].append(count)
self._coupling2mol.append(
((mon1, init1, fin1), (mon2, init2, fin2))
)
count += 1
[docs]
def init_coupling_vector(self) -> None:
"""initialize coupling vector"""
self.init_coupling_indexes()
Ncoupl = len(self._coupling2mol)
self.resonance_coupling_vec = numpy.zeros(Ncoupl, dtype=numpy.float64)
# TODO: add initialization flag
self.coupling_initiated = True
[docs]
def get_resonance_coupling_vec(
self, mon1: int, init1: int, fin1: int, mon2: int, init2: int, fin2: int
) -> float | numpy.ndarray:
"""returns coupling between mon1 transition init1->fin1 and
mon2 transition init2->fin2
Parameters
-----------
mon1,mon2 : integer
indexes of two monomers betweem which interaction energy is calculated
init1,fin1 : integer
index of initial and final state, respectively, of monomer 1 for
electronic transition init1->fin1. Ground state is 0 first excited
state 1, and so on
init2,fin2 : integer
index of initial and final state, respectively, of monomer 2
Return
---------
coupling : float
coupling between mon1 transition init1->fin1 and
mon2 transition init2->fin2 in current energy units
"""
if mon1 == mon2:
return 0.0
if mon1 > mon2:
# exchange excitations between monomers
mon2, mon1 = mon1, mon2
init2, init1 = init1, init2
fin2, fin1 = fin1, fin2
indx1 = mon1
indx4 = mon2 - mon1 - 1
if init1 < fin1:
indx2 = init1
indx3 = fin1 - init1 - 1
else:
indx2 = fin1
indx3 = init1 - fin1 - 1
if init2 < fin2:
indx5 = init2
indx6 = fin2 - init2 - 1
else:
indx5 = fin2
indx6 = init2 - fin2 - 1
coupling_indx = self._mol2coupling[indx1][indx2][indx3][indx4][indx5][indx6]
coupling = self.resonance_coupling_vec[coupling_indx]
return self.convert_energy_2_current_u(coupling)
[docs]
def set_coupling_by_dipole_dipole_vec(self, epsr: float = 1.0) -> None:
"""Sets resonance coupling by dipole-dipole interaction for multilevel
molecules
"""
if not self.coupling_initiated:
self.init_coupling_vector()
self.init_coupling_matrix()
for kk in range(self.resonance_coupling_vec.size):
monomer1, monomer2 = self._coupling2mol[kk]
cc = self.dipole_dipole_coupling_multilevel(
monomer1[0],
monomer1[1],
monomer1[2],
monomer2[0],
monomer2[1],
monomer2[2],
epsr=epsr,
)
c1 = self.convert_energy_2_internal_u(cc)
self.resonance_coupling_vec[kk] = c1
if ((monomer1[1] == 0) and (monomer1[2] == 1)) and (
(monomer2[1] == 0) and (monomer2[2] == 1)
):
self.resonance_coupling[monomer1[0], monomer2[0]] = c1
self.resonance_coupling[monomer2[0], monomer1[0]] = c1
[docs]
def init_coupling_matrix(self) -> None:
"""Nullifies coupling matrix"""
nstates = 0
# count the number of electronic states
for monomer in self.monomers:
nstates += monomer.nel - 1
self.resonance_coupling = numpy.zeros((nstates, nstates), dtype=numpy.float64)
self.coupling_initiated = True
#
# TESTED
[docs]
def set_resonance_coupling(
self,
i: int,
j: int,
coupling: float,
mode_linear: list | None = None,
mode_shift: list | None = None,
) -> None:
"""Sets resonance coupling value between two sites"""
if not self.coupling_initiated:
######################## My new part beginning #################################### @Vladislav Slama
self.init_coupling_matrix()
self.init_coupling_vector()
self._set_coupling_vec(i, 0, 1, j, 0, 1, coupling)
#
# This has to stay for compatibility reasons
#
coup = self.convert_energy_2_internal_u(coupling)
self.resonance_coupling[i, j] = coup
self.resonance_coupling[j, i] = coup
#
# TESTED
#
# Coordinate dependence of the resonance coupling
#
# this is a temporary solution
self.mode_linear = mode_linear
if mode_linear is not None:
self.mode_shift = len(mode_linear) * [0.0]
self.has_mode_dependence = True
if mode_shift is not None:
self.mode_shift = mode_shift
def _set_coupling_vec(
self,
mon1: int,
init1: int,
fin1: int,
mon2: int,
init2: int,
fin2: int,
coupling: float,
) -> None:
"""Sets resonance coupling value between two transitions. Works
for multilevel molecules
Parameters
-----------
mon1,mon2 : integer
indexes of two monomers betweem which interaction energy is calculated
init1,fin1 : integer
index of initial and final state, respectively, of monomer 1 for
electronic transition init1->fin1. Ground state is 0 first excited
state 1, and so on
init2,fin2 : integer
index of initial and final state, respectively, of monomer 2
"""
if not self.coupling_initiated:
self.init_coupling_vector()
if mon1 == mon2:
print(
"Trying to set couling between states within the same "
"molecule. This coupling is set to zero!"
)
return
if mon1 > mon2:
# exchange excitations between monomers
mon2, mon1 = mon1, mon2
init2, init1 = init1, init2
fin2, fin1 = fin1, fin2
indx1 = mon1
indx4 = mon2 - mon1 - 1
if init1 < fin1:
indx2 = init1
indx3 = fin1 - init1 - 1
else:
indx2 = fin1
indx3 = init1 - fin1 - 1
if init2 < fin2:
indx5 = init2
indx6 = fin2 - init2 - 1
else:
indx5 = fin2
indx6 = init2 - fin2 - 1
coupling_indx = self._mol2coupling[indx1][indx2][indx3][indx4][indx5][indx6]
self.resonance_coupling_vec[coupling_indx] = self.convert_energy_2_internal_u(
coupling
)
return
[docs]
def set_coupling_by_Hamiltonian(self, HH: numpy.ndarray) -> None:
"""set the resonance coupling values according to given electronic
Hamiltonian (single exciton band is expected - without the ground state).
Only excitations from the ground state are allowed.
Parameters
-----------
HH : numpy.array of float
First exciton band hamiltonian. Only excitations from the ground
state are supported so far. The ordering of the hamiltonian must
be the same and the nomomers and its transitions.
"""
# FIXME: Allow higher excited states for the monomer
count = [0, 0]
Nmon = len(self.monomers)
for mon1indx in range(Nmon):
mon1 = self.monomers[mon1indx]
for elst1 in range(1, mon1.nel): # without the ground state
count[1] = 0
for mon2indx in range(Nmon):
mon2 = self.monomers[mon2indx]
for elst2 in range(1, mon2.nel):
if mon2indx > mon1indx:
coupling = HH[count[0], count[1]]
self._set_coupling_vec(
mon1indx, 0, elst1, mon2indx, 0, elst2, coupling
)
count[1] += 1
count[0] += 1
[docs]
def set_resonance_coupling_vec(self, i: int, j: int, coupling: float) -> None:
"""Sets resonance coupling value between two sites (between !electronic
states!)
Parameters
----------
i,j : integer
indexes of electronic levels (in aggregate) for which the coupling
is set
coupling : float
coupling between the electronic states in current energy units
"""
if not self.coupling_initiated:
self.init_coupling_vector()
# TODD: Delete unnesesary
self.init_coupling_matrix()
# TODO: Delete unnesesary
coup = self.convert_energy_2_internal_u(coupling)
self.resonance_coupling[i, j] = coup
self.resonance_coupling[j, i] = coup
text_warning = "Trying to define coupling between states within single molecule"
if self.nmono > 1:
elst1 = self.elsigs[i]
elst2 = self.elsigs[j]
if self.which_band[i] == self.which_band[j]:
if self.which_band[i] == 1:
mon1 = numpy.nonzero(elst1)[0][0]
mon2 = numpy.nonzero(elst2)[0][0]
if mon1 == mon2:
print(text_warning)
else:
fin1 = elst1[mon1]
fin2 = elst2[mon2]
self._set_coupling_vec(mon1, 0, fin1, mon2, 0, fin2, coupling)
else:
Ns = len(elst1)
sites = [0, 0]
n_diff = 0
# count differences
for i in range(Ns):
if elst1[i] != elst2[i]:
if (n_diff == 0) or (n_diff == 1):
sites[n_diff] = i
n_diff += 1
# if there are exactly 2 differences, the differing
# two molecules are those coupled; sites[n_diff] contains
# indiced those coupled molecules
if n_diff == 2:
mon1 = sites[0]
mon2 = sites[1]
init1 = elst1[mon1]
fin1 = elst2[mon1]
init2 = elst1[mon2]
fin2 = elst2[mon2]
self._set_coupling_vec(
mon1, init1, fin1, mon2, init2, fin2, coupling
)
else:
print(text_warning)
else:
print(text_warning)
else:
print(text_warning)
[docs]
def get_resonance_coupling(self, i: int, j: int) -> float | numpy.ndarray:
"""Returns resonance coupling value between two sites"""
return self.get_resonance_coupling_vec(i, 0, 1, j, 1, 0)
[docs]
def set_resonance_coupling_matrix(
self, coupmat: numpy.ndarray | list | tuple
) -> None:
"""Sets resonance coupling values from a matrix"""
if type(coupmat) in (list, tuple):
coupmat = numpy.array(coupmat)
coup = self.convert_energy_2_internal_u(numpy.asarray(coupmat))
self.resonance_coupling = coup
if not self.coupling_initiated:
self.coupling_initiated = True
#
# TESTED
[docs]
def dipole_dipole_coupling(
self, kk: int, ll: int, epsr: float = 1.0, delta: float = 1.0e-5
) -> float | numpy.ndarray:
"""Calculates dipole-dipole coupling"""
if kk == ll:
raise QuantarheiError(
"Only coupling between different molecules \
can be calculated"
)
# FIXME: this works only for first excited states of two-level molecules
d1 = self.monomers[kk].dmoments[0, 1, :]
r1 = self.monomers[kk].position
d2 = self.monomers[ll].dmoments[0, 1, :]
r2 = self.monomers[ll].position
if numpy.sqrt(numpy.dot(r1 - r2, r1 - r2)) < delta:
raise QuantarheiError()
val = dipole_dipole_interaction(r1, r2, d1, d2, epsr)
return self.convert_energy_2_current_u(val)
#
# TESTED
[docs]
def dipole_dipole_coupling_multilevel(
self,
mon1: int,
in1: int,
fin1: int,
mon2: int,
in2: int,
fin2: int,
epsr: float = 1.0,
) -> float | numpy.ndarray:
"""Calculates dipole-dipole coupling for multilevel monomers
Parameters
-----------
mon1,mon2 : integer
indexes of two monomers betweem which dipole-dipole interaction
energy is calculated
init1,fin1 : integer
index of initial and final state, respectively, of monomer 1 for
electronic transition init1->fin1. Ground state is 0 first excited
state 1, and so on
init2,fin2 : integer
index of initial and final state, respectively, of monomer 2
epsr : float
relative permitivity of the environment
Returns
----------
val : float
dipole-dipole interaction energy in current energy units
"""
if mon1 == mon2:
raise QuantarheiError(
"Only coupling between different molecules \
can be calculated"
)
# FIXME: this works only for first excited states of two-level molecules
d1 = self.monomers[mon1].dmoments[in1, fin1, :]
r1 = self.monomers[mon1].position
d2 = self.monomers[mon2].dmoments[in2, fin2, :]
r2 = self.monomers[mon2].position
val = dipole_dipole_interaction(r1, r2, d1, d2, epsr)
return self.convert_energy_2_current_u(val)
[docs]
def set_coupling_by_dipole_dipole(
self, epsr: float = 1.0, delta: float = 1.0e-5
) -> None:
"""Sets resonance coupling by dipole-dipole interaction"""
self.set_coupling_by_dipole_dipole_vec(epsr=epsr)
[docs]
def calculate_resonance_coupling(
self, method: str = "dipole-dipole", params: dict | None = None
) -> None:
"""Sets resonance coupling calculated by a given method
Parameters
----------
method: string
Method to be used for calculation of resonance coupling
"""
if params is None:
params = {"epsr": 1.0}
if method == "dipole-dipole":
epsr = params["epsr"]
self.set_coupling_by_dipole_dipole(epsr=epsr)
else:
raise QuantarheiError(
"Unknown method for calculation of resonance coupling"
)
#
# TESTED
# FIXME: This must be set in coordination with objects describing laboratory
[docs]
def set_lindich_axes(self, axis_orthog_membrane: numpy.ndarray) -> None:
"""Creates a coordinate system with one axis supplied by the user
(typically an axis orthogonal to the membrane), and two other axes, all
of which are orthonormal.
"""
qr = numpy.vstack(
(axis_orthog_membrane, numpy.array([1, 0, 0]), numpy.array([0, 1, 0]))
).T
self.q, r = numpy.linalg.qr(qr)
self._has_lindich_axes = True
# FIXME: This must be set in coordination with objects describing laboratory
[docs]
def get_lindich_axes(self) -> numpy.ndarray:
if self._has_lindich_axes:
return self.q
raise QuantarheiError("No linear dichroism coordinate system supplied")
# FIXME: This should be delegated SystemBathInteraction
[docs]
def set_egcf_matrix(self, cm: object) -> None:
"""Sets a matrix describing system bath interaction"""
self.egcf_matrix = cm
self._has_egcf_matrix = True
#
# Molecues
#
[docs]
def add_Molecule(self, mono: Any) -> None:
"""Adds monomer to the aggregate"""
# If at least one monomer has energy gap correlation function
# we will try to build system bath interaction for a the aggregate.
# Exception will be thrown if not all monomers have the same egcf
if mono._has_egcf:
self._has_system_bath_interaction = True
self.monomers.append(mono)
self.mnames[mono.name] = len(self.monomers) - 1
self.nmono += 1
#
# TESTED
[docs]
def get_Molecule_by_name(self, name: str) -> object:
try:
im = self.mnames[name]
return self.monomers[im]
except (KeyError, IndexError):
raise QuantarheiError()
[docs]
def get_Molecule_index(self, name: str) -> int:
try:
im = self.mnames[name]
return im
except KeyError:
raise QuantarheiError()
[docs]
def remove_Molecule(self, mono: object) -> None:
self.monomers.remove(mono)
self.nmono -= 1
[docs]
def get_nearest_Molecule(self, molecule: Any) -> tuple:
"""Returns a molecule nearest in the aggregate to a given molecule
Parameters
----------
molecule : Molecule
Molecule whose neighbor we look for
"""
tol = 1.0e-3
rmin = 1.0e20
r1 = molecule.position
mmin = None
for m in self.monomers:
r2 = m.position
r_vec = r1 - r2
dist = numpy.sqrt(numpy.dot(r_vec, r_vec))
if (dist > tol) and (dist < rmin):
mmin = m
rmin = dist
return mmin, rmin
#
# Vibrational modes
#
#
# Transition dipole
#
[docs]
def get_dipole(
self, n: Any, N: Any = None, M: Any = None, **kwargs: Any
) -> numpy.ndarray:
nm = self.monomers[n]
return nm.get_dipole(N, M)
[docs]
def get_velocity_dipole(self, n: int, N: int, M: int) -> numpy.ndarray:
nm = self.monomers[n]
return nm.get_velocity_dipole(N, M)
[docs]
def get_magnetic_dipole(self, n: int, N: int, M: int) -> numpy.ndarray:
nm = self.monomers[n]
return nm.get_magnetic_dipole(N, M)
#
# Various info
#
[docs]
def get_width(self, n: int, N: int, M: int) -> float:
nm = self.monomers[n]
return nm.get_transition_width((N, M))
[docs]
def get_max_excitations(self) -> list:
"""Returns a list of maximum number of excitations on each monomer"""
omax = []
for nm in self.monomers:
omax.append(nm.nel - 1)
return omax
[docs]
def set_dipole(self, N: Any, M: Any = None, vec: Any = None) -> None:
"""Sets the dipole moment of a given transition"""
raise QuantarheiError(
"Transition dipole moments cannot be set directly. Use aggregate components."
)
[docs]
def fc_factor(self, state1: Any, state2: Any) -> float:
"""Franck-Condon factors between two vibrational states
Calculates Franck-Condon factor between two aggregate_states
regardless of their electronic parts
"""
inx1 = state1.vsig
inx2 = state2.vsig
sta1 = state1.elstate.vibmodes
sta2 = state2.elstate.vibmodes
if not (len(sta1) == len(sta2)):
raise QuantarheiError("Incompatible states")
res = 1.0
for kk in range(len(sta1)):
smod1 = sta1[kk]
smod2 = sta2[kk]
# difference in shifts
shft = smod1.shift - smod2.shift
# quantum numbers
qn1 = inx1[kk]
qn2 = inx2[kk]
# calculate FC factors
#
# Best implementation would be a table look-up. First we calculate
# a table of FC factors from known omegas and shifts and here we
# just consult the table.
if not self.FC.lookup(shft):
fc = self.ops.shift_operator(shft)[:20, :20]
self.FC.add(shft, fc)
ii = self.FC.index(shft)
rs = self.FC.get(ii)[qn1, qn2]
res = res * rs
return numpy.real(res)
[docs]
def map_egcf_to_states(self, mpx: numpy.ndarray) -> numpy.ndarray:
"""Maps the participation matrix on g(t) functions storage
This should work for a simple mapping matrix and first excited band.
Specialized mapping such as those for 2 exciton states is implemented
in classes that inherite from here.
"""
ss = self.SS
Ng = self.Nb[0]
Ne1 = self.Nb[1] + Ng
def _nonzero(vec: numpy.ndarray) -> int | None:
"""Returns a possition in the vector on which value == 1 is found"""
for kk, vl in enumerate(vec):
if vl == 1:
return kk
return None
if self.mult > 1:
Ne2 = self.Nb[2] + Ne1
if self.mult == 1:
#
# This works only for aggregate without vibrations
# FIXME: generalize to vibrations
WPM = numpy.einsum(
"na,nb,ni->abi", ss[Ng:Ne1, Ng:Ne1] ** 2, ss[Ng:Ne1, Ng:Ne1] ** 2, mpx
)
elif self.mult == 2:
# version including higher excited state band ( excitons etc.)
Del = numpy.zeros((Ne2 - Ng, Ne2 - Ng, mpx.shape[1]), dtype=REAL)
# looping over the states
for ii, state1 in self.allstates(mult=self.mult):
for jj, state2 in self.allstates(mult=self.mult):
sts = state1.get_shared_sites(state2)
Del[ii - Ng, jj - Ng, :] = numpy.einsum("ij->j", mpx[sts[:], :])
# if len(sts) == 1:
# # one site is shared
# Del[ii-Ng,jj-Ng,:] = mpx[sts[0],:]
#
# elif len(sts) == 2:
# # two sites shared (for two-excitons states are equal)
# Del[ii-Ng,jj-Ng,:] = mpx[sts[0],:] + mpx[sts[1],:]
WPM = numpy.einsum(
"nmk,an,bm->abk", Del, ss[Ng:Ne2, Ng:Ne2] ** 2, ss[Ng:Ne2, Ng:Ne2] ** 2
)
"""
# we assume that the mapping is onto sites and we construct
# two-exiciton mapping here. In this mapping we count
# how many times a given correlation function ocurs in a give
# double-exciton.
mpx2 = numpy.zeros((Ne2-Ng,mpx.shape[1]),dtype=numpy.int32)
mpx2[:(Ne1-Ng),:mpx.shape[1]] = mpx
for kk, state in self.allstates(mult=2):
#estate = state.get_ElectronicState()
if kk >= Ne1:
sts = state.get_excited_sites()
f_index_st0 = _nonzero(mpx[sts[0],:])
f_index_st1 = _nonzero(mpx[sts[1],:])
# we add one to the correlation function whose
# index we found
mpx2[kk-Ng, f_index_st0] += 1
mpx2[kk-Ng, f_index_st1] += 1
print("mpx2:")
print(mpx2)
#WPM = numpy.einsum("na,nb,ni->abi",ss[Ng:Ne2,Ng:Ne2]**2,
# ss[Ng:Ne2,Ng:Ne2]**2,mpx2)
WPM = numpy.einsum("na,nb->ab",ss[Ng:Ne2,Ng:Ne2]**2,
ss[Ng:Ne2,Ng:Ne2]**2)
print(WPM)
print(WPM[2,0,0])
print(WPM[2,0,1])
print(WPM[2,1,0])
print(WPM[2,1,1])
"""
else:
raise QuantarheiError(
"Participation matrix not implemented for"
"multiplicity higher than mult=2."
)
return WPM
[docs]
def get_transition_width(self, state1: Any, state2: Any | None = None) -> float:
"""Returns phenomenological width of a given transition
Parameters
----------
state1 : {ElectroniState/VibronicState, tuple}
If both state1 and state2 are specified, it is assumed they are
of the type of Electronic of Vibronic state. Otherwise, if state2
is None, it is assumed that it is a tuple representing
a transition
state2 : {ElectroniState/VibronicState, None}
If not None it is of the type of Electronic of Vibronic state
"""
if state2 is not None:
b1 = state1.elstate.band
b2 = state2.elstate.band
if abs(b2 - b1) == 1:
# index of a monomer on which the transition occurs
exindx = self._get_exindx(state1, state2)
# TODO: Here it should not be 01 transition but for multilevel molecule appropriate transition
exct1 = state1.elstate.elsignature[exindx]
exct2 = state2.elstate.elsignature[exindx]
width = self.monomers[exindx].get_transition_width((exct2, exct1))
return width
if abs(b2 - b1) == 2:
s1_signature = state1.elstate.elsignature
if numpy.nonzero(s1_signature)[0].size == 1:
# Todo: repare _get_exindx for state from double excited block but on single molecule
exindx = self._get_exindx(state1, state2)
exct = state1.elstate.elsignature[exindx]
width = self.monomers[exindx].get_transition_width((0, exct))
# --------------------- My version is more general and takes the width defined externaly
# sig1 = state1.elstate.get_signature()
# sig2 = state2.elstate.get_signature()
# if (2 in sig1) or (2 in sig2):
# exindx = self._get_exindx(state1, state2)
# width = \
# 2.0*self.monomers[exindx].get_transition_width((0,1))
# -------------------------------------------------------------------------------
return width
twoex = self._get_twoexindx(state1, state2)
if isinstance(twoex, int):
return -1.0
(indx1, indx2) = twoex
# state2.elstate.elsignature, indx1, indx2)
exct1 = state1.elstate.elsignature[indx1]
exct2 = state1.elstate.elsignature[indx2]
width = self.monomers[indx1].get_transition_width((0, exct1))
width += self.monomers[indx2].get_transition_width((0, exct2))
return width
return -1.0
transition = state1
Nf = transition[0]
Ni = transition[1]
eli = self.elinds[Ni]
elf = self.elinds[Nf]
# g -> 1 exciton band transitions
if (self.which_band[eli] == 0) and (self.which_band[elf] == 1):
# this simulates bath correlation function
return self.Wd[Nf, Nf] ** 2
if (self.which_band[eli] == 1) and (self.which_band[elf] == 0):
# this simulates bath correlation function
return self.Wd[Ni, Ni] ** 2
# 1 exciton -> 2 exciton transitions
if (self.which_band[eli] == 1) and (self.which_band[elf] == 2):
# this simulates the term g_ff + g_ee - 2Re g_fe
ret = (
self.Wd[Ni, Ni] ** 2
+ self.Wd[Nf, Nf] ** 2
- 2.0 * (self.Wd[Nf, Ni] ** 2)
)
return ret
if (self.which_band[eli] == 2) and (self.which_band[elf] == 1):
# this simulates the term g_ff + g_ee - 2Re g_fe
ret = (
self.Wd[Ni, Ni] ** 2
+ self.Wd[Nf, Nf] ** 2
- 2.0 * (self.Wd[Nf, Ni] ** 2)
)
return ret
return 0.0
[docs]
def get_transition_dephasing(self, state1: Any, state2: Any | None = None) -> float:
"""Returns phenomenological dephasing of a given transition
Parameters
----------
state1 : {ElectroniState/VibronicState, tuple}
If both state1 and state2 are specified, it is assumed they are
of the type of Electronic of Vibronic state. Otherwise, if state2
is None, it is assumed that it is a tuple representing
a transition
state2 : {ElectroniState/VibronicState, None}
If not None it is of the type of Electronic of Vibronic state
"""
if state2 is not None:
# index of a monomer on which the transition occurs
exindx = self._get_exindx(state1, state2)
if exindx < 0:
return 0.0
deph = self.monomers[exindx].get_transition_dephasing((0, 1))
return deph
transition = state1
Nf = transition[0]
Ni = transition[1]
eli = self.elinds[Ni]
elf = self.elinds[Nf]
# g -> 1 exciton band transitions
if (self.which_band[eli] == 0) and (self.which_band[elf] == 1):
return self.Dr[Nf, Nf] ** 2
if (self.which_band[eli] == 1) and (self.which_band[elf] == 0):
return self.Dr[Ni, Ni] ** 2
# 1 exciton -> 2 exciton band transitions
if (self.which_band[eli] == 1) and (self.which_band[elf] == 2):
# this simulates the term g_ff + g_ee - 2Re g_fe
return self.Dr[Ni, Ni] ** 2 + self.Dr[Nf, Nf] - 2.0 * self.Dr[Nf, Ni]
if (self.which_band[eli] == 2) and (self.which_band[elf] == 1):
# this simulates the term g_ff + g_ee - 2Re g_fe
return self.Dr[Ni, Ni] ** 2 + self.Dr[Nf, Nf] - 2.0 * self.Dr[Nf, Ni]
return -1.0
[docs]
def transition_dipole(self, state1: Any, state2: Any) -> numpy.ndarray | float:
"""Transition dipole moment between two states
Parameters
----------
state1 : class VibronicState
state 1
state2 : class VibronicState
state 2
"""
exindx = self._get_exindx(state1, state2)
if exindx < 0:
return 0.0
######################## My new part beginning #################################### @Vladislav Slama
# get excitation indexes
st1 = state1.elstate.elsignature[exindx]
st2 = state2.elstate.elsignature[exindx]
if st1 < st2:
eldip = self.get_dipole(exindx, st1, st2)
else:
eldip = self.get_dipole(exindx, st2, st1)
# Franck-Condon factor between the two states
fcfac = self.fc_factor(state1, state2)
return eldip * fcfac
[docs]
def transition_velocity_dipole(
self, state1: Any, state2: Any
) -> numpy.ndarray | float:
"""Transition dipole moment between two states
Parameters
----------
state1 : class VibronicState
state 1
state2 : class VibronicState
state 2
"""
exindx = self._get_exindx(state1, state2)
if exindx < 0:
return 0.0
# get excitation indexes
st1 = state1.elstate.elsignature[exindx]
st2 = state2.elstate.elsignature[exindx]
if st1 < st2:
elvdip = self.get_velocity_dipole(exindx, st1, st2)
else:
elvdip = self.get_velocity_dipole(exindx, st2, st1)
# Franck-Condon factor between the two states
fcfac = self.fc_factor(state1, state2)
return elvdip * fcfac
[docs]
def transition_magnetic(self, state1: Any, state2: Any) -> numpy.ndarray | float:
"""Transition magnetic dipole moment between two states
Parameters
----------
state1 : class VibronicState
state 1
state2 : class VibronicState
state 2
"""
exindx = self._get_exindx(state1, state2)
if exindx < 0:
return 0.0
# get excitation indexes
st1 = state1.elstate.elsignature[exindx]
st2 = state2.elstate.elsignature[exindx]
magdip = self.get_magnetic_dipole(exindx, st1, st2)
# Franck-Condon factor between the two states
fcfac = self.fc_factor(state1, state2)
return magdip * fcfac
def _get_twoexindx(self, state1: Any, state2: Any) -> tuple | int:
"""Indices of two molecule with transitions or negative number
if not found
Parameters
----------
state1 : class VibronicState
state 1
state2 : class VibronicState
state 2
"""
# get electronic signatures
els1 = state1.elstate.elsignature
els2 = state2.elstate.elsignature
# only states in neighboring bands can be connected by dipole moment
b1 = state1.elstate.band
b2 = state2.elstate.band
if abs(b1 - b2) != 2:
return -1
# count the number of differences
sig_idx = 0
count = 0
for kk in els1:
if kk != els2[sig_idx]:
count += 1
sig_idx += 1
if count != 2:
return -1
# now that we know that the states differ by two excitations, let
# us find on which molecule they are
exstates = []
exindxs = []
sig_idx = -1
for kk in els1: # signature is just a tuple; iterate over it
sig_idx += 1
if kk != els2[sig_idx]: # this is the index where they differ
# which of them is excited
if kk > els2[sig_idx]:
exstates.append(els1)
else:
exstates.append(els2)
exindxs.append(sig_idx)
if len(exstates) == 0:
raise QuantarheiError()
return exindxs[0], exindxs[1]
def _get_exindx(self, state1: Any, state2: Any) -> int:
"""Index of molecule with transition or negative number if not found
Parameters
----------
state1 : class VibronicState
state 1
state2 : class VibronicState
state 2
"""
# get electronic signatures
els1 = state1.elstate.elsignature
els2 = state2.elstate.elsignature
# count the number of differences
sig_idx = 0
count = 0
for kk in els1:
if kk != els2[sig_idx]:
count += 1
sig_idx += 1
if count != 1:
return -1
# now that we know that the states differ by one excitation, let
# us find on which molecule it is
exstate = None
sig_idx = -1
for kk in els1: # signature is just a tuple; iterate over it
sig_idx += 1
if kk != els2[sig_idx]: # this is the index where they differ
# which of them is excited
if kk > els2[sig_idx]:
exstate = els1
else:
exstate = els2
exindx = sig_idx
if exstate is None:
raise QuantarheiError()
return exindx
[docs]
def total_number_of_states(
self,
mult: int = 1,
vibgen_approx: str | None = None,
Nvib: int | None = None,
vibenergy_cutoff: float | None = None,
save_indices: bool = False,
band_external: list | None = None,
) -> int:
"""Total number of states in the aggregate
Counts all states of the aggregate by iterating through them. States
are generated with a set of constraints.
"""
nret = 0
for state in self.allstates(
mult=mult,
save_indices=save_indices,
vibgen_approx=vibgen_approx,
Nvib=Nvib,
vibenergy_cutoff=vibenergy_cutoff,
band_external=band_external,
):
nret += 1
return nret
[docs]
def total_number_of_electronic_states(self, mult: int = 1) -> int:
"""Total number of electronic states in the aggregate"""
nret = 0
for elsig in self.elsignatures(mult=mult):
nret += 1
return nret
[docs]
def number_of_states_in_band(
self,
band: int = 1,
vibgen_approx: str | None = None,
Nvib: int | None = None,
vibenergy_cutoff: float | None = None,
) -> int:
"""Number of states in a given excitonic band"""
nret = 0
for state in self.allstates(
mult=band,
mode="EQ",
save_indices=False,
vibgen_approx=vibgen_approx,
Nvib=Nvib,
vibenergy_cutoff=vibenergy_cutoff,
):
nret += 1
return nret
[docs]
def number_of_electronic_states_in_band(self, band: int = 1) -> int:
"""Number of states in a given excitonic band"""
nret = 0
for elsig in self.elsignatures(mult=band, mode="EQ"):
nv = 1
nret += nv
return nret
[docs]
def get_ElectronicState(
self, sig: tuple, index: int | None = None
) -> ElectronicState:
"""Returns electronic state corresponding to this aggregate
Parameters
----------
sig : tuple
Tuple defining the electronic state of the aggregate
index : integer or None
If integer is specified, this number is recorded as an index
of this state in the aggregate. It is recorded only internally
in the state object. Aggregate keeps its own record which is
created during the build.
"""
return ElectronicState(self, sig, index)
[docs]
def get_StateBand(self, state: Any) -> int:
"""Returns band of the corresponding electronic state
This effectively counts the number of excitations in the state
by counting the band numbers in each molecules.
Parameters
----------
state : ElectronicState
Aggregate electronic state.
"""
band = 0
elsig = state.elsignature
for n in range(len(elsig)):
mon = self.monomers[n]
band += mon.which_band[elsig[n]]
return band
[docs]
def get_VibronicState(self, esig: tuple, vsig: tuple) -> VibronicState:
"""Returns vibronic state corresponding to the two specified signatures"""
elstate = self.get_ElectronicState(sig=esig)
return VibronicState(elstate, vsig)
[docs]
def coupling_vec(self, state1: Any, state2: Any) -> float | numpy.ndarray:
"""Coupling between two aggregate states
Parameters
----------
state1,state2 : {ElectronicState, VibronicState}
States for which coupling should be calculated
Returns
---------
coup : float
Resonance coupling in current units
"""
#
# Coupling between two purely electronic states
#
coup: float | numpy.ndarray = 0.0
if isinstance(state1, ElectronicState) and isinstance(state2, ElectronicState):
if self.nmono > 1:
# coupling within the bands
if state1.band == state2.band:
els1 = state1.elsignature
els2 = state2.elsignature
if state1.band == 1:
mon1 = numpy.nonzero(els1)[0][0]
mon2 = numpy.nonzero(els2)[0][0]
if mon1 == mon2:
coup = 0.0
else:
fin1 = els1[mon1]
fin2 = els2[mon2]
coup = self.get_resonance_coupling_vec(
mon1, 0, fin1, mon2, 0, fin2
)
else:
Ns = len(els1)
sites = [0, 0]
n_diff = 0
# count differences
for i in range(Ns):
if els1[i] != els2[i]:
if (n_diff == 0) or (n_diff == 1):
sites[n_diff] = i
n_diff += 1
# if there are exactly 2 differences, the differing
# two molecules are those coupled; sites[n_diff] contains
# indiced those coupled molecules
if n_diff == 2:
mon1 = sites[0]
mon2 = sites[1]
init1 = els1[mon1]
fin1 = els2[mon1]
init2 = els1[mon2]
fin2 = els2[mon2]
coup = self.get_resonance_coupling_vec(
mon1, init1, fin1, mon2, init2, fin2
)
#
# kk = sites[0]
# ll = sites[1]
# coup = self.resonance_coupling[kk,ll]
else:
coup = 0.0
else:
coup = 0.0
else:
coup = 0.0
#
# Coupling between two general states
#
elif isinstance(state1, VibronicState) and isinstance(state2, VibronicState):
es1 = state1.elstate
es2 = state2.elstate
fc = self.fc_factor(state1, state2)
# it make sense to calculate coupling only when the number
# of molecules is larger than 1
if self.nmono > 1:
# coupling within the bands
if es1.band == es2.band:
els1 = es1.elsignature
els2 = es2.elsignature
# single exciton band
if es1.band == 1:
mon1 = numpy.nonzero(els1)[0][0]
mon2 = numpy.nonzero(els2)[0][0]
if mon1 == mon2:
coup = 0.0
else:
fin1 = els1[mon1]
fin2 = els2[mon2]
coup = self.get_resonance_coupling_vec(
mon1, 0, fin1, mon2, 0, fin2
)
coup *= fc
else:
Ns = len(els1)
sites = [0, 0]
n_diff = 0
# count differences
for i in range(Ns):
if els1[i] != els2[i]:
if (n_diff == 0) or (n_diff == 1):
sites[n_diff] = i
n_diff += 1
# if there are exactly 2 differences, the differing
# two molecules are those coupled; sites[n_diff] contains
# indiced those coupled molecules
if n_diff == 2:
mon1 = sites[0]
mon2 = sites[1]
init1 = els1[mon1]
fin1 = els2[mon1]
init2 = els1[mon2]
fin2 = els2[mon2]
coup = self.get_resonance_coupling_vec(
mon1, init1, fin1, mon2, init2, fin2
)
coup *= fc
else:
coup = 0.0
else:
coup = 0.0
else:
coup = 0.0
return self.convert_energy_2_current_u(coup)
[docs]
def coupling(
self, state1: Any, state2: Any, full: bool = False
) -> float | numpy.ndarray:
"""Coupling between two aggregate states
Parameters
----------
state1 : {ElectronicState, VibronicState}
States for which coupling should be calculated
"""
#
# Coupling between two purely electronic states
#
if isinstance(state1, ElectronicState) and isinstance(state2, ElectronicState):
if self.nmono > 1:
# coupling within the bands
if state1.band == state2.band:
if state1.band == 1:
kk = state1.index - 1
ll = state2.index - 1
if (kk >= 0) and (ll >= 0):
coup = self.resonance_coupling[kk, ll]
else:
coup = 0.0
else:
els1 = state1.elsignature
els2 = state2.elsignature
Ns = len(els1)
sites = [0, 0]
n_diff = 0
# count differences
for i in range(Ns):
if els1[i] != els2[i]:
if (n_diff == 0) or (n_diff == 1):
sites[n_diff] = i
n_diff += 1
# if there are exactly 2 differences, the differing
# two molecules are those coupled; sites[n_diff] contains
# indiced those coupled molecules
if n_diff == 2:
kk = sites[0]
ll = sites[1]
coup = self.resonance_coupling[kk, ll]
else:
coup = 0.0
elif state1.band + 2 == state2.band:
coup = 0.0
else:
coup = 0.0
else:
coup = 0.0
#
# Coupling between two general states
#
elif isinstance(state1, VibronicState) and isinstance(state2, VibronicState):
es1 = state1.elstate
es2 = state2.elstate
fc = self.fc_factor(state1, state2)
# it make sense to calculate coupling only when the number
# of molecules is larger than 1
if self.nmono > 1:
# coupling within the bands
if es1.band == es2.band:
# single exciton band
if es1.band == 1:
kk = es1.index - 1
ll = es2.index - 1
if (kk >= 0) and (ll >= 0):
coup = self.resonance_coupling[kk, ll] * fc
else:
coup = 0.0
else:
els1 = es1.elsignature
els2 = es2.elsignature
Ns = len(els1)
sites = [0, 0]
n_diff = 0
# count differences
for i in range(Ns):
if els1[i] != els2[i]:
if (n_diff == 0) or (n_diff == 1):
sites[n_diff] = i
n_diff += 1
# if there are exactly 2 differences, the differing
# two molecules are those coupled; sites[n_diff] contains
# indiced those coupled molecules
if n_diff == 2:
kk = sites[0]
ll = sites[1]
ar1 = numpy.array(els1)
ar2 = numpy.array(els2)
df = numpy.abs(ar1 - ar2)
sdf = numpy.sum(df)
if sdf == 2:
mx1 = numpy.max([ar1[kk], ar2[kk]])
mx2 = numpy.max([ar1[ll], ar2[ll]])
harm_fc = numpy.sqrt(numpy.real(mx1))
harm_fc = harm_fc * numpy.sqrt(numpy.real(mx2))
fc = fc * harm_fc
coup = self.resonance_coupling[kk, ll] * fc
else:
coup = 0.0
else:
coup = 0.0
elif (numpy.abs(es1.band - es2.band) == 2) and full:
els1 = es1.elsignature
els2 = es2.elsignature
Ns = len(els1)
sites = [0, 0]
n_diff = 0
# count differences
for i in range(Ns):
if els1[i] != els2[i]:
if (n_diff == 0) or (n_diff == 1):
sites[n_diff] = i
n_diff += 1
# if there are exactly 2 differences, the differing
# two molecules are those coupled; sites[n_diff] contains
# indiced those coupled molecules
if n_diff == 2:
kk = sites[0]
ll = sites[1]
coup = self.resonance_coupling[kk, ll] * fc
else:
coup = 0.0
else:
coup = 0.0
else:
coup = 0.0
return self.convert_energy_2_current_u(coup)
#######################################################################
#
# Generators
#
#######################################################################
[docs]
def elsignatures(
self, mult: int = 1, mode: str = "LQ", emax: list | None = None
) -> Iterator[tuple[Any, ...]]:
"""Generator of electronic signatures
Here we create signature tuples of electronic states. The signature
is a tuple with as many integer numbers as the members of
the aggregate. Each integer represents the state in which the
member of the aggregate is, e.g. 0 for ground state, 1 for the first
excited state etc.
Provide correct ordering for the multilevel molecues
Parameters
----------
mult : int
multiplicity of excitons
mode : str {"LQ", "EQ"}
mode of the functions.
mode="LQ" returns all signatures of states with
multiplicity less than or equal to the `mult`
mode="EQ" returns signatures of states with a multiplicity
given by `mult`
"""
if mode not in ["LQ", "EQ"]:
raise QuantarheiError("Unknown mode")
n_sites = len(self.monomers)
# list of maximum numbers of excitations on each sites
if emax is None:
omax = self.get_max_excitations()
else:
omax = emax
if mult < 0:
raise QuantarheiError("mult must be larger than or equal to zero")
# define internal multiplicity which is dependent on the monomers
for band_n in range(mult + 1):
mult_int = band_n
for mon in self.monomers:
Nexct = (
numpy.sum(mon.Nb[: band_n + 1]) - (mon.Nb[: band_n + 1]).size
) # Count extra states in the bands
mult_int += Nexct
mlt = 0
# iterate over all excition multiplicities
while mlt <= mult_int:
# no excitations (ground state)
out = [0 for k in range(n_sites)]
# if this is the multiplicity 0, yield the ground state
# if (((mlt == 0) and (mode == "LQ")) or (mult==0)):
# if band_n==0:
if ((band_n == 0) and (mode == "LQ")) or (mult == 0):
yield tuple(out)
else:
excit_lvl = 1
# first we have only ground state signature
ins = [out]
strt = [0]
while excit_lvl <= mlt:
nins = []
nstr = []
# take all signatures in "ins" and add one excitation
for out_added, last in self._add_excitation(ins, strt, omax):
# if mlt excitation was added yield
if ((excit_lvl == mlt) and (mode == "LQ")) or (
(mult_int == excit_lvl) and (mult_int == mlt)
):
band = 0
for ii in range(len(out_added)):
band += self.monomers[ii].which_band[out_added[ii]]
if band == band_n:
yield tuple(out_added)
else:
# make a list of all new signatures
nins.append(out_added)
# for each signature, save the index
# on which an excitation was added last
nstr.append(last)
# set the new signatures for processing in the iteration
ins = nins
strt = nstr
excit_lvl += 1
mlt += 1
def _add_excitation(
self, inlists: list, strt: list, omax: list
) -> Iterator[tuple[Any, Any]]:
"""Adds one excitation to all submitted electronic signatures"""
list_idx = 0
# go through all signatures
for inlist in inlists:
n_pos = len(inlist)
if len(omax) != n_pos:
raise QuantarheiError(
"arg omax has to be a list of the same \
length as arg inlist"
)
# go through all positions from the last index on (in order
# to create unique signatures)
for i in range(strt[list_idx], n_pos):
# if it is possible to add an excitation
# make a new list and add
if inlist[i] < omax[i]:
out = inlist.copy()
out[i] += 1
# yield the list and the index of the last added exitation
yield out, i
list_idx += 1
[docs]
def vibsignatures(self, elsignature: tuple, approx: str | None = None) -> Any:
"""Generator of vibrational signatures
Parameters
----------
approx : None or str
Approximation used in generation of vibrational states
Allowed values are None or "SPA"
"""
cs = ElectronicState(self, elsignature)
return cs.vsignatures(approx=approx)
[docs]
def allstates(
self,
mult: int = 1,
mode: str = "LQ",
all_vibronic: bool = True,
save_indices: bool = False,
vibgen_approx: str | None = None,
Nvib: int | None = None,
vibenergy_cutoff: float | None = None,
band_external: list | None = None,
) -> Generator[tuple[int, Any], None, None]:
"""Generator of all aggregate states
Iterator generating all states of the aggregate given a set
of constraints.
Parameters
----------
mult : integer {0, 1, 2}
Exciton multiplicity (ground state, single and double excitons).
All excitons with the multiplicity smaller or equal to ``mult``
are generated by default
mode : str {"LQ", "EQ"}
If set to "LQ" generates excitons with smaller or equal
multiplicity than specified. If "EQ" is specified, generates only
excitons with given multiplicity
save_indices : bool
If True, saves indices of all generated states, so that they can
be later used.
all_vibronic : bool
If True, all generated states are of the type ``VibronicState``,
even if no vibrational modes are specified. If False,
``ElectronicState`` is returned for pure electronic states
vibgen_approx : str {"ZPA", "SPA", "TPA", "NPA", "SPPMA", "TPPMA", "NPPMA"}
Type of approximation in generating vibrational states
Nvib : integer
Number of vibrational states that goes into "NPA" and "NPPMA"
approximations
vibenergy_cutoff: float
Maximum vibrational energy allowed in generation of vibrational
states
"""
ast = 0 # index counting all states
ist = 0 # index counting electronic states
######################## My new part beginning #################################### @Vladislav Slama
# create list of electronic signatures
es_list = []
bands = []
for ess1 in self.elsignatures(mult=mult, mode=mode):
es1 = self.get_ElectronicState(ess1, ist)
es_list.append(es1)
if band_external is None:
es1.band = self.get_StateBand(
es1
) # this should produce right bands without external definition
bands.append(es1.band)
ist += 1
# If needed reorder according to the external band definition
if band_external is not None:
# change band for electronic system
for band in band_external:
bands[band[0]] = band[1]
es_list[band[0]].band = band[1]
new_indx = numpy.argsort(bands, kind="mergesort")
# reorder according to bands
es_list = [es_list[ii] for ii in new_indx]
else:
new_indx = numpy.argsort(bands, kind="mergesort")
# reorder according to bands
es_list = [es_list[ii] for ii in new_indx]
bands = [bands[ii] for ii in new_indx]
#
# ist = 0 # index counting electronic states
# run over all electronic states
# for ess1 in self.elsignatures(mult=mult, mode=mode):
for nn, es1 in enumerate(es_list):
#
# # generate electronic state
# es1 = self.get_ElectronicState(ess1, ist)
ess1 = es1.elsignature
ist = es1.index
if bands[nn] <= mult:
# loop over all vibrational signatures in electronic states
nsig = 0
for vsig1 in es1.vsignatures(
approx=vibgen_approx, N=Nvib, vibenergy_cutoff=vibenergy_cutoff
):
# create vibronic state with a given signature
s1: VibronicState | ElectronicState = VibronicState(es1, vsig1)
s1.band = bands[nn]
# s1.elstate.band = bands[nn]
if save_indices:
# save indices corresponding to vibrational sublevels
# of a given electronic state
self.vibindices[ist].append(ast)
self.vibsigs[ast] = (ess1, vsig1)
self.elinds[ast] = ist
yield ast, s1
ast += 1 # count all states
nsig += 1 # count the number of vibrational signatures
# if no vibrational signatures
if nsig == 0:
# if True return vibronic states even
# for purely electronic state
if all_vibronic:
s1 = VibronicState(es1, None)
else:
s1 = es1
if save_indices:
# save electronic signatures to be searchable later
self.elsigs[ist] = ess1
# save the band to which this electronic index corresponds
if band_external is None:
self.which_band[ist] = bands[nn]
else:
self.which_band[ist] = numpy.sum(ess1)
[docs]
def elstates(
self, mult: int = 1, mode: str = "LQ", save_indices: bool = False
) -> Generator[tuple[int, Any], None, None]:
"""Generator of electronic states"""
a = 0
for ess1 in self.elsignatures(mult=mult, mode=mode):
es1 = self.get_ElectronicState(ess1, a)
yield a, es1
a += 1
def __repr__(self) -> str:
return f"Aggregate(name={self.name!r}, nmono={self.nmono})"
def __str__(self) -> str:
out = "\nquantarhei.Aggregate object"
out += "\n==========================="
out += f"\nname = {self.name}"
out += f"\nnumber of molecules = {self.nmono:d} "
count = 0
for nm in self.monomers:
out += f"\n\nMonomer {count:d}"
out += str(nm)
count += 1
out += "\n\nResonance coupling matrix: "
out += "\n-------------------------- "
out += "\n" + str(self.resonance_coupling)
out += "\n\nAggregate built = " + str(self._built)
out += "\n\nSelected attributes"
out += "\n--------------------"
out += "\nmult = " + str(self.mult)
out += "\nNel = " + str(self.Nel)
out += "\nNtot = " + str(self.Ntot)
return out
###########################################################################
#
# BUILDING
#
###########################################################################
[docs]
def build(
self,
mult: int = 1,
sbi_for_higher_ex: bool = False,
vibgen_approx: str | None = None,
Nvib: int | None = None,
vibenergy_cutoff: float | None = None,
band_external: list | None = None,
el_blocks: bool = False,
) -> None:
# fem_full=False):
# TODO: Add Full Frenkel exciton model
"""Builds aggregate properties
Calculates Hamiltonian and transition dipole moment matrices and
sets up system-bath interaction
Parameters
----------
mult : int
exciton multiplicity
sbi_for_higher_ex: bool
If set True, system-bath information is explicitely created for
higher exciton states (consistent with the specified parameters
`mult`). If set False, it is expected that if system-bath
interaction for higher excitons is needed, it will be reconstructed
from the single exciton part of this object
vibge_approx:
Approximation used in the generation of vibrational state.
band_external : list of integers (dimension Nx2)
redefinition of correspondence of individual electronic states to
excitonic bands, e.g. [[4,1],[5,2]] means that 4th electronic state
is moved to single exciton band and 5th electronic state
to the second exciton band (zeroth electronic state is the ground
state). The mult must be high enough to allow the building the
original states, e.g for 2 two level carotenoids where S2-S2
interaction should mult=2 Because it is interaction between second
excited state and secodn excited state
NOTES:
--------------
good to check the bands before building the molecule with:
for a, s1 in aggreg.allstates(mult=mult,
vibgen_approx=vibgen_approx, Nvib=Nvib,
vibenergy_cutoff=vibenergy_cutoff):
print(a,"elsign:",s1.elstate.elsignature,"vibsign:",s1.vsig,"elindex:",s1.elstate.index,"band:",s1.elstate.band)
"""
manager = Manager()
manager.set_current_units("energy", "int")
# maximum multiplicity of excitons handled by this aggregate
self.mult = mult
if sbi_for_higher_ex:
self.sbi_mult = mult
else:
self.sbi_mult = 1
#######################################################################
#
# Electronic and vibrational states
#
#######################################################################
# total number of electronic states
self.Nel = self.total_number_of_electronic_states(mult=mult)
# storage for indices of vibrational states
self.vibindices = []
# there are as many lists of indices as there are electronic states
for i in range(self.Nel):
self.vibindices.append([])
# number of states in the aggregate (taking into account
# approximations in generation of vibrational states)
Ntot = self.total_number_of_states(
mult=mult,
vibgen_approx=vibgen_approx,
Nvib=Nvib,
save_indices=False,
vibenergy_cutoff=vibenergy_cutoff,
)
# save total number of states (including vibrational)
self.Ntot = Ntot
# information about the band to which a state belongs
# self.which_band = numpy.zeros(self.Ntot, dtype=int)
self.which_band = numpy.zeros(self.Nel, dtype=int)
# electronic signature for every state
self.elsigs = [None] * self.Nel
# vibrational signature for each state
self.vibsigs = [None] * self.Ntot
# FIXME: what is this???
self.elinds = numpy.zeros(self.Ntot, dtype=int)
# Hamiltonian matrix
HH = numpy.zeros((Ntot, Ntot), dtype=numpy.float64)
# Transition dipole moment matrix
DD = numpy.zeros((Ntot, Ntot, 3), dtype=numpy.float64)
# Magnetic dipole moment matrix (in coordinate system centered on the molecule)
MM = numpy.zeros((Ntot, Ntot, 3), dtype=numpy.complex128)
# Rotatory strength matrix
RR = numpy.zeros((Ntot, Ntot), dtype=numpy.float64)
RRv = numpy.zeros((Ntot, Ntot), dtype=numpy.float64)
RRm = numpy.zeros((Ntot, Ntot), dtype=numpy.float64)
# Matrix of Franck-Condon factors
FC = numpy.zeros((Ntot, Ntot), dtype=numpy.float64)
# Matrix of the transition widths (their square roots)
Wd: numpy.ndarray = numpy.zeros((Ntot, Ntot), dtype=REAL)
# Matrix of dephasing rates
Dr: numpy.ndarray = numpy.zeros((Ntot, Ntot), dtype=REAL)
# Matrix of dephasing transformation coefficients
self.Xi: numpy.ndarray = numpy.zeros((Ntot, self.Nel), dtype=REAL)
# electronic indices of twice excited state (zero for all other states)
twoex_indx = numpy.zeros((Ntot, 2), dtype=int)
# two-exciton state index by pair of single excitations
# - originally this was only alocated with the number of molecules/monomers
twoex_state = numpy.zeros((self.Nel, self.Nel), dtype=int)
# Initialization of the matrix of couplings between states
if not self.coupling_initiated:
# FIXME: Duplicate. coupling_vector is more general
self.init_coupling_matrix()
self.init_coupling_vector()
Ntot = self.total_number_of_states(
mult=mult,
vibgen_approx=vibgen_approx,
Nvib=Nvib,
save_indices=True,
vibenergy_cutoff=vibenergy_cutoff,
band_external=band_external,
)
# repair band assignments for electronic states
if band_external is not None:
for ii in range(self.Ntot):
elindx = self.elinds[ii]
for newband in band_external:
if elindx == newband[0]:
self.which_band[elindx] = newband[1]
self.all_states = []
for a, s1 in self.allstates(
mult=self.mult,
vibgen_approx=vibgen_approx,
Nvib=Nvib,
vibenergy_cutoff=vibenergy_cutoff,
band_external=band_external,
):
self.all_states.append((a, s1))
# if el_blocks:
# self.electronic_blocks = dict()
# for sig in self.elsigs:
# self.electronic_blocks[sig] = "ciao"
#
#
# Set up Hamiltonian and Transition dipole moment matrices
for a, s1 in self.all_states:
if a == 0:
s0 = s1
# diagonal Hamiltonian elements
######################## My new part beginning #################################### @Vladislav Slama
HH[a, a] = s1._energy() # for energy in current units use s1.energy()
# get dephasing and width from the ground-state
# for each excited state
elind = self.elinds[a]
if (self.which_band[elind] == 1) or (self.which_band[elind] == 2):
Wd[a, a] = numpy.sqrt(self.get_transition_width(s1, s0))
Dr[a, a] = numpy.sqrt(self.get_transition_dephasing(s1, s0))
#
# save composition of twice excited states
# (This only stores info about aggregate states composed of excitations on different molecules)
#
if self.which_band[elind] == 2:
# k_s counts excited molecules in the doubly exc. state
# there are molecules 0 and 1 in diad (n,m)
k_s = 0
# counts positons in the electronic signature
# i.e. it counts molecular index
sig_position = 0
# we loop through the states in the signature
for i_s in s1.elstate.elsignature:
# record excitations
# if i_s == 1:
if i_s != 0: # allow the doubly excited states on single molecule
# we save indices of electronic states and
# 0 is taken by the ground state
state_is_list = [
0,
] * self.nmono
state_is_list[sig_position] = i_s
state_is = tuple(state_is_list)
els = self.get_ElectronicState(state_is)
# two alternatives (older vs. newer)
# twoex_indx[a, k_s] = sig_position + 1
assert els.index is not None
twoex_indx[a, k_s] = els.index
k_s += 1
sig_position += 1
#
# If this is a two-exciton state, then save the info here
#
if (twoex_indx[a, 0] - 1 >= 0) and (twoex_indx[a, 1] - 1 >= 0):
twoex_state[twoex_indx[a, 0] - 1, twoex_indx[a, 1] - 1] = a
twoex_state[twoex_indx[a, 1] - 1, twoex_indx[a, 0] - 1] = a
# This works for single member bands
# self.twoex_state = twoex_state
for b, s2 in self.all_states: # self.allstates(mult=self.mult,
# vibgen_approx=vibgen_approx, Nvib=Nvib,
# vibenergy_cutoff=vibenergy_cutoff):
DD[a, b, :] = numpy.real(self.transition_dipole(s1, s2))
MM[a, b, :] = self.transition_magnetic(s1, s2)
FC[a, b] = numpy.real(self.fc_factor(s1, s2))
# FIXME: Here we assume only excitation from the lowest state (lowest vibrational state)
mon1 = s1.get_monomer()
mon2 = s2.get_monomer()
try:
if mon1 != -1 and mon2 != -1:
da = self.transition_dipole(s0, s1)
db = self.transition_dipole(s0, s2)
mb = self.transition_magnetic(s0, s2)
Ra = numpy.array(self.monomers[mon1].position, "f8")
# alternative definition of rotatory strength
# Rb = numpy.array(self.monomers[mon2].position,"f8")
# RR[a,b] = numpy.dot( (Ra - Rb), numpy.cross(da, db))
Ea = s1._energy() - s0._energy()
# for energy in current units use s1.energy()
try:
dav = self.transition_velocity_dipole(s0, s1)
self._has_velocity_dipoles = True
except (AttributeError, TypeError):
dav = -1j * Ea * da
RRv[a, b] = numpy.real(1j * numpy.dot(Ra, numpy.cross(dav, db)))
RR[a, b] = numpy.dot(Ra, numpy.cross(da, db))
RRm[a, b] = numpy.real(numpy.dot(dav, mb))
except Exception:
pass
if a != b:
HH[a, b] = self.coupling_vec(s1, s2)
# TODO: ADD FULL Frenkel exciton model
# HH[a,b] = numpy.real(self.coupling(s1, s2, full=fem_full))
####### New !!!!CHECK!!!!!
# nz1 = numpy.nonzero(s1.elstate.elsignature)[0]
# nz2 = numpy.nonzero(s2.elstate.elsignature)[0]
# if (nz1.size == 1) and (nz2.size == 1) and (nz1 == nz2) and (a!=b):
# if s1.elstate.elsignature[nz1[0]] < s2.elstate.elsignature[nz2[0]]:
# trwidth = self.get_transition_width(s2, s1)
# else:
# trwidth = 0.0
#
# if trwidth >= 0:
# Wd[b,a] = numpy.sqrt(trwidth)
# else:
# Wd[b,a] = 0.0
self.twoex_state = twoex_state
####### End new
# Storing Hamiltonian and dipole moment matrices
self.HH = HH
# Hamiltonian operator
self.HamOp = Hamiltonian(data=HH)
# dipole moments
self.DD = DD
# dipole moments
self.MM = MM
# Oscilatory strength
self.RR = RR
self.RRv = RRv
self.RRm = RRm
# FIXME: make this on-demand (if poissible)
trdata: numpy.ndarray = numpy.zeros(
(DD.shape[0], DD.shape[1], DD.shape[2]), dtype=REAL
)
trdata[:, :, :] = DD[:, :, :]
self.TrDMOp = TransitionDipoleMoment(data=trdata)
self.TrMMOp = TransitionDipoleMoment(data=MM)
# Franck-Condon factors
self.FCf = FC
# widths
self.Wd = Wd
# dephasings
self.Dr = Dr
# composition of two-ex states
# first index of state a is twoex_indx[0, a]
self.twoex_indx = twoex_indx
# squares of transition dipoles
dd2 = numpy.sum(self.DD**2, axis=2)
self.D2 = dd2
# FIXME: do I need this??? Is it even corrrect??? (maybe amax?)
# maximum of transition dipole moment elements
self.D2_max = numpy.max(dd2)
# Number of states in individual bands
self.Nb = numpy.zeros(self.mult + 1, dtype=int)
# for ii in range(self.mult+1):
# self.Nb[ii] = self.number_of_states_in_band(band=ii,
# vibgen_approx=vibgen_approx, Nvib=Nvib,
# vibenergy_cutoff=vibenergy_cutoff)
for ii in range(self.Ntot):
elindex = self.elinds[ii]
band = self.which_band[elindex]
self.Nb[band] += 1
# Number of electronic states in individual bands
self.Nbe = numpy.zeros(self.mult + 1, dtype=int)
# for ii in range(self.mult+1):
# self.Nbe[ii] = self.number_of_electronic_states_in_band(band=ii)
for ii in range(self.Nel):
band = self.which_band[ii]
self.Nbe[band] += 1
# prepare RWA indices and set info for Rotating Wave Approximation
rwa_indices = numpy.zeros(self.mult + 1, int)
for ii in range(self.mult):
rwa_indices[ii + 1] = rwa_indices[ii] + self.Nb[ii]
self.HamOp.set_rwa(rwa_indices)
engtmp = numpy.diag(self.HH)
self.vibenergy = engtmp.copy() # to keep site basis info also after diag.
self.vibdipoles = self.DD.copy() # to keep site basis info also after diag.
indxsorted = numpy.argsort(engtmp)
self.vibsigs_engsort = [None] * self.Ntot
for ii in range(self.Ntot):
self.vibsigs_engsort[ii] = self.vibsigs[indxsorted[ii]]
#######################################################################
#
# System-bath interaction
#
#######################################################################
#
# There are two methods to set system-bath interaction
# 1) Each molecule gets its bath correlation function
# 2) Global energy gap correlation function matrix is set
#
# is energy gap correlation function matrix present?
if self._has_egcf_matrix:
nmonst = 0
for monomer in self.monomers:
nmonst += monomer.nel - 1
# Check the consistency of the energy gap correlation matrix
if self.egcf_matrix.nob != nmonst:
raise QuantarheiError(
"Correlation matrix has a size different"
" from the number of monomeric states"
)
# FIXME The aggregate having a egcf matrix does not mean the monomers
# have egcf matrices. They could just have correlation funtions.
for i in range(self.nmono):
if (
self.monomers[i]._is_mapped_on_egcf_matrix
and self.monomers[i].egcf_matrix is not self.egcf_matrix
):
# TODO: Ask about this - it would mean that every monomer has the same energy gap correlation function
raise QuantarheiError(
"Correlation matrix in the monomer"
" has to be the same as the one of"
" the aggregate."
)
# seems like everything is consistent -> we can calculate system-
# -bath interaction
self._has_system_bath_interaction = True
# if not, try to get one from monomers later
else:
self._has_system_bath_interaction = False
# try to set energy gap correlation matrix from monomers
if not self._has_system_bath_interaction:
# let's assume we can calculate EGCF matrix from monomers
egcf_ok = True
# get correlation function from a monomer
try:
egcf1 = self.monomers[0].get_transition_environment((0, 1))
except Exception:
# we cannot calculate EGCF matrix, there is no system-bath
# interaction, or it is not based on correlation functions
egcf_ok = False
# if we have correlation functions for nonomers, let's construct
# the EGCF matrix
if egcf_ok:
# time axis of the first monomer
time = egcf1.axis
# Number of correlation functions is the number of electronic
# states minus ground state (this assumes that only electronic
# states are coupled to the bath)
Nelg = 1 # ASSUMPTION: here we assume a single electronic
# ground state
if sbi_for_higher_ex:
# except for ground state, all electronic states have EGCF
Ncf = self.Nel - Nelg
else:
# in single exciton, two-level molecule picture, there is
# a single correlation function per monomer
# ASSUMPTION: Two-level molecules
# Ncf = self.nmono
# Ncf = self.Nel - Nelg # We construct bigger correlation
# matrix, but fill only part
Ncf = self.Nbe[1]
# instantiate the EGCF matrix object
self.egcf_matrix = CorrelationFunctionMatrix(time, Ncf)
# run over all electronic states
for i in range(self.Nel):
elsig = self.elsigs[i] # eletronic signature of state i
nzr = numpy.nonzero(elsig)[0] # all nonzero elements of the
# signature
nnzr = nzr.size # number of nonzero elements in the
# signature
# in single exciton band
if self.which_band[i] == 1 or nnzr == 1:
j = i - Nelg
# mon = self.monomers[j]
mon = self.monomers[nzr[0]]
exct = elsig[nzr[0]]
# get correlation for a monomer
# ASSUMPTION: Two-level molecule
cfce = mon.get_transition_environment((0, exct))
# set correlation function into the diagonal of the
# EGCF matrix. Index corresponds to the monomer
mapi = self.egcf_matrix.set_correlation_function(cfce, [(j, j)])
# FIXME: what is returned?
if mapi <= 0:
raise QuantarheiError("Something's wrong")
# in two-exciton band
elif (self.which_band[i] == 2) and sbi_for_higher_ex and nnzr != 1:
egcf_idx = i - Nelg
# monomers of a two-exciton state are obtaines
# FIXME: is this correct???
# j = self.elsigs[i][0]
# k = self.elsigs[i][1]
# mon1 = self.monomers[j]
# mon2 = self.monomers[k]
mon1 = self.monomers[nzr[0]]
mon2 = self.monomers[nzr[1]]
exct1 = elsig[nzr[0]]
exct2 = elsig[nzr[1]]
# we get correlation functions of the two monomers
# ASSUMPTION: Two-level molecules
cfce1 = mon1.get_transition_environment((0, exct1))
cfce2 = mon2.get_transition_environment((0, exct2))
# correlation functions are added to form a two-exciton
# correlation function
cfce = cfce1 + cfce2
# Two-exciton correlation function is set into
# EGCF matrix
mapi = self.egcf_matrix.set_correlation_function(
cfce, [(egcf_idx, egcf_idx)]
)
# FIXME: cross-correlation between double excitons
# needs to be handled.
for k2 in range(self.Nel):
# Index for the correlation functions (without the ground states)
egcf_idx2 = k2 - Nelg
# find which monomer(s) is(are) excited
elsig2 = self.elsigs[k2] # eletronic signature of state i
nzr2 = numpy.nonzero(elsig2)[0]
nnzr2 = nzr2.size
if self.which_band[k2] == 1:
# The cross correlation function will be the one with common index
# which monomer excited
mon2 = self.monomers[nzr2[0]]
# to which state it is excited
exct2 = elsig2[nzr2[0]]
if k2 in self.twoex_indx[i]:
cfce = mon2.get_transition_environment((0, exct2))
mapi = self.egcf_matrix.set_correlation_function(
cfce, [(egcf_idx, egcf_idx2)]
)
mapi = self.egcf_matrix.set_correlation_function(
cfce, [(egcf_idx2, egcf_idx)]
)
if (
(self.which_band[k2] == 2)
and (nnzr2 != 1)
and (i != k2)
):
mon21 = self.monomers[nzr2[0]]
mon22 = self.monomers[nzr2[1]]
exct21 = elsig2[nzr2[0]]
exct22 = elsig2[nzr2[1]]
elexcit1 = self.twoex_indx[i]
elexcit2 = self.twoex_indx[k2]
if (
elexcit1[0] == elexcit2[0]
or elexcit1[1] == elexcit2[0]
):
cfce = mon21.get_transition_environment((0, exct21))
mapi = self.egcf_matrix.set_correlation_function(
cfce, [(egcf_idx, egcf_idx2)]
)
elif (
elexcit1[0] == elexcit2[1]
or elexcit1[1] == elexcit2[1]
):
cfce = mon22.get_transition_environment((0, exct22))
mapi = self.egcf_matrix.set_correlation_function(
cfce, [(egcf_idx, egcf_idx2)]
)
# TODO:
# Add else:
# cfce = 0
# mapi = self.egcf_matrix.set_correlation_function(cfce,[(egcf_idx,egcf_idx2)])
if mapi <= 0:
raise QuantarheiError("Something's wrong")
# some effective theory here
# FIXME: make sure we know what the ``sbi_for_higher_ex``
# switch actually means
elif (self.which_band[i] == 2) and (not sbi_for_higher_ex):
# this should be handled by
# a map between double excitons and site cor. functions
pass
# no theory for higher bands so-far
elif (self.which_band[i] > 2) and sbi_for_higher_ex:
pass
self._has_system_bath_interaction = True
self._has_egcf_matrix = True
# TODO: Repair the following for multilevel systems
# if all needed for system-bath interaction is present
# we can construct the SystemBathInteraction object
if self._has_system_bath_interaction:
# system interaction operators
iops = []
# how many operators should be created
if sbi_for_higher_ex:
Nop = self.Nel - 1 # all electronic states
else:
Nop = self.Nbe[1] # we count only single excited states
# count electronic states of monomeric basis
nmonst = 0
for monomer in self.monomers:
nmonst += monomer.nel - 1
# if there are more states in the single exciton block
# than the number of sites, it means we have vibrational states
if nmonst != self.Nb[1]:
# create a projection operator for each monomer
# a monomer corresponds to one single excited state starting
# with electronic index 1 (0 is the ground state)
# ASSUMPTION: Two-level molecules
for i in range(1, Nop + 1):
op1 = Operator(dim=self.HH.shape[0], real=True)
# here we make a projector on a given electronic state |i>
# ASSUMPTION: Oscillator is represented by its eigenstates
# FIXME: if it should be projector on electronic states there should be elinds
# instead of vibindexes
for j in self.vibindices[i]:
op1.data[j, j] = 1.0
iops.append(op1)
# standard case with only electronic states
else:
# create a projection operator for each monomer
# a monomer corresponds to one single excited state starting
# with electronic index 1 (0 is the ground state)
# ASSUMPTION: Two-level molecules
for i in range(1, Nop + 1):
op1 = Operator(dim=self.HH.shape[0], real=True)
op1.data[i, i] = 1.0
iops.append(op1)
# we create SystemBathInteraction object
self.sbi = SystemBathInteraction(iops, self.egcf_matrix, system=self)
else:
# system-bath interaction is not present
pass
self._built = True
manager.unset_current_units("energy")
[docs]
def rebuild(
self,
mult: int = 1,
sbi_for_higher_ex: bool = False,
vibgen_approx: str | None = None,
Nvib: int | None = None,
vibenergy_cutoff: float | None = None,
) -> None:
"""Cleans the object and rebuilds it"""
self.clean()
self.build(
mult=mult,
sbi_for_higher_ex=sbi_for_higher_ex,
vibgen_approx=vibgen_approx,
Nvib=Nvib,
vibenergy_cutoff=vibenergy_cutoff,
)
###########################################################################
#
# POST BUILDING METHODS
#
###########################################################################
def _allocate_converted_operator(
self,
operator: (
ReducedDensityMatrix
| DensityMatrix
| Hamiltonian
| TransitionDipoleMoment
| ReducedDensityMatrixEvolution
| DensityMatrixEvolution
),
dim: int,
Nt: int | None = None,
allow_tdm: bool = False,
) -> tuple[
ReducedDensityMatrix
| Hamiltonian
| TransitionDipoleMoment
| ReducedDensityMatrixEvolution,
int,
bool,
bool,
]:
n_indices = 2
evolution = False
whole = False
nop: (
ReducedDensityMatrix
| Hamiltonian
| TransitionDipoleMoment
| ReducedDensityMatrixEvolution
| None
) = None
if isinstance(operator, (ReducedDensityMatrix, DensityMatrix)):
nop = ReducedDensityMatrix(dim=dim)
elif isinstance(operator, Hamiltonian):
nop = Hamiltonian(dim=dim)
elif allow_tdm and isinstance(operator, TransitionDipoleMoment):
nop = TransitionDipoleMoment(dim=dim)
n_indices = 3
elif isinstance(
operator, (ReducedDensityMatrixEvolution, DensityMatrixEvolution)
):
if Nt is not None:
nop = ReducedDensityMatrix(dim=dim)
evolution = True
else:
nop = ReducedDensityMatrixEvolution(operator.TimeAxis)
rhoi = ReducedDensityMatrix(dim=dim)
nop.set_initial_condition(rhoi)
evolution = True
whole = True
else:
raise Exception(
"Operation not implemented for this type: "
+ operator.__class__.__name__
)
return nop, n_indices, evolution, whole
def _fc_basis_transform(
self,
nop: (
ReducedDensityMatrix
| Hamiltonian
| TransitionDipoleMoment
| ReducedDensityMatrixEvolution
),
operator: (
ReducedDensityMatrix
| DensityMatrix
| Hamiltonian
| TransitionDipoleMoment
| ReducedDensityMatrixEvolution
| DensityMatrixEvolution
),
n_components: int,
) -> None:
for a in range(n_components):
for n in range(self.Nel):
i_ng = -1
for i_n in self.vibindices[n]:
i_ng += 1
for m in range(self.Nel):
i_mg = -1
for i_m in self.vibindices[m]:
i_mg += 1
if n_components > 1:
nop._data[i_n, i_m, a] = 0.0
for j_n in self.vibindices[n]:
for j_m in self.vibindices[m]:
nop._data[i_n, i_m, a] += (
self.FCf[i_ng, j_n]
* operator._data[j_n, j_m, a]
* self.FCf[j_m, i_mg]
)
else:
nop._data[i_n, i_m] = 0.0
for j_n in self.vibindices[n]:
for j_m in self.vibindices[m]:
nop._data[i_n, i_m] += (
self.FCf[i_ng, j_n]
* operator._data[j_n, j_m]
* self.FCf[j_m, i_mg]
)
[docs]
def convert_to_ground_vibbasis(
self,
operator: (
ReducedDensityMatrix
| DensityMatrix
| Hamiltonian
| TransitionDipoleMoment
| ReducedDensityMatrixEvolution
| DensityMatrixEvolution
),
Nt: int | None = None,
) -> (
ReducedDensityMatrix
| Hamiltonian
| TransitionDipoleMoment
| ReducedDensityMatrixEvolution
):
"""Converts an operator to a ground state vibrational basis repre
Default representation in Quantarhei is that with a specific shifted
vibrational basis in each electronic state. Here we convert to a
representation where there is a single basis used for all vibrational
states regardless of elecrtronic state
The conversion MUST be done in site basis. Only in site basis
we can distinguish the vibrational states properly
"""
if operator.dim == self.Ntot:
nop, n_indices, evolution, whole = self._allocate_converted_operator(
operator, dim=self.Ntot, Nt=Nt, allow_tdm=True
)
# FIXME: This limitation might not be necessary
# in the ground states of all monomers, there must be the same
# or greater number of levels than in the excited state
# over all monomers
for k in range(self.nmono):
mono = self.monomers[k]
# over all modes
n_mod = mono.get_number_of_modes()
for i in range(n_mod):
mod = mono.get_Mode(i)
n_g = mod.get_nmax(0)
# over all excited states
# FIXME: this should be mono.Nel as in Aggregate
for j in range(mono.nel):
if j > 0:
n_e = mod.get_nmax(j)
if n_e > n_g:
raise QuantarheiError(
"Number of levels"
" in the excited state of a molecule has to be \n"
"the same or smaller than in the ground state"
)
#
# ground state vibrational states
#
stgs = []
for i_g in self.vibindices[0]:
vs_g = self.vibsigs[i_g]
stg = self.get_VibronicState(vs_g[0], vs_g[1])
stgs.append(stg)
if n_indices in (2, 3):
n_components = 3 if n_indices == 3 else 1
self._fc_basis_transform(nop, operator, n_components)
return nop
raise QuantarheiError("Incompatible operator")
raise Exception("Incompatible operator dimension")
[docs]
def trace_converted(
self,
operator: (
ReducedDensityMatrix
| DensityMatrix
| Hamiltonian
| ReducedDensityMatrixEvolution
| DensityMatrixEvolution
),
Nt: int | None = None,
) -> (
ReducedDensityMatrix
| Hamiltonian
| TransitionDipoleMoment
| ReducedDensityMatrixEvolution
):
if operator.dim == self.Ntot:
nop, n_indices, evolution, whole = self._allocate_converted_operator(
operator, dim=self.Nel, Nt=Nt
)
if n_indices == 2:
for n in range(self.Nel):
i_ng = -1
for i_n in self.vibindices[n]:
i_ng += 1
for m in range(self.Nel):
i_mg = -1
for i_m in self.vibindices[m]:
i_mg += 1
if i_ng == i_mg:
nop._data[n, m] += operator._data[i_n, i_m]
return nop
raise QuantarheiError("Incompatible operator")
raise Exception("Incompatible operator dimension")
def _accumulate_fc_trace(
self,
nop: ReducedDensityMatrix
| Hamiltonian
| TransitionDipoleMoment
| ReducedDensityMatrixEvolution,
op_slice: numpy.ndarray,
FcProd: numpy.ndarray,
whole: bool,
) -> None:
for n in range(self.Nel):
for i_n in self.vibindices[n]:
for m in range(self.Nel):
for i_m in self.vibindices[m]:
fc = FcProd[i_n, i_m]
if whole:
nop._data[:, n, m] += op_slice[:, i_n, i_m] * fc
else:
nop._data[n, m] += op_slice[0, i_n, i_m] * fc
[docs]
def trace_over_vibrations(
self,
operator: (
ReducedDensityMatrix
| DensityMatrix
| Hamiltonian
| ReducedDensityMatrixEvolution
| DensityMatrixEvolution
),
Nt: int | None = None,
) -> (
ReducedDensityMatrix
| Hamiltonian
| TransitionDipoleMoment
| ReducedDensityMatrixEvolution
):
"""Average an operator over vibrational degrees of freedom
Average MUST be done in site basis. Only in site basis
we can distinguish the vibrational states properly
"""
if operator.dim == self.Ntot:
nop, n_indices, evolution, whole = self._allocate_converted_operator(
operator, dim=self.Nel, Nt=Nt
)
if n_indices == 2:
# convert to representation by ground-state oscillator
# FIXME: This limitation might not be necessary
# in the ground states of all monomers, there must be the same
# or greater number of levels than in the excited state
# over all monomers
for k in range(self.nmono):
mono = self.monomers[k]
# over all modes
n_mod = mono.get_number_of_modes()
for i in range(n_mod):
mod = mono.get_Mode(i)
n_g = mod.get_nmax(0)
# over all excited states
# FIXME: this should be mono.Nel as in Aggregate
for j in range(mono.nel):
if j > 0:
n_e = mod.get_nmax(j)
if n_e > n_g:
raise QuantarheiError(
"Number of levels"
" in the excited state of a molecule has to be \n"
"the same or smaller than in the ground state"
)
# do the conversion
#
# ground state vibrational states
#
stgs = []
for i_g in self.vibindices[0]:
vs_g = self.vibsigs[i_g]
stg = self.get_VibronicState(vs_g[0], vs_g[1])
stgs.append(stg)
Ng = self.Nb[0]
FcProd = numpy.einsum("gi,jg->ij", self.FCf[:Ng, :], self.FCf[:, :Ng])
if evolution and whole:
op_slice = operator._data
elif evolution:
op_slice = operator._data[Nt, :, :][numpy.newaxis, :, :]
else:
op_slice = operator._data[numpy.newaxis, :, :]
self._accumulate_fc_trace(nop, op_slice, FcProd, evolution and whole)
else:
raise QuantarheiError(
"Cannot trace over this object: wrong number of indices"
)
return nop
raise QuantarheiError("Incompatible operator")
[docs]
def cast_to_vibronic(self, KK: numpy.ndarray) -> numpy.ndarray:
"""Casts an electronic operator to a vibronic basis"""
agg = self
newkk = numpy.zeros((agg.Ntot, agg.Ntot), dtype=numpy.float64)
# populate the operator
for i_el in range(agg.Nel):
for i_vib in agg.vibindices[i_el]:
vs_i = agg.vibsigs[i_vib]
st_i = agg.get_VibronicState(vs_i[0], vs_i[1])
for j_el in range(agg.Nel):
for j_vib in agg.vibindices[j_el]:
vs_j = agg.vibsigs[j_vib]
st_j = agg.get_VibronicState(vs_j[0], vs_j[1])
# electronic transition operator
# dressed in Franck-Condon factors
newkk[i_vib, j_vib] = (
numpy.real(agg.fc_factor(st_i, st_j)) * KK[i_el, j_el]
)
return newkk
[docs]
def convert_to_DensityMatrix(
self, psi: object, trace_over_vibrations: bool = True
) -> object:
"""Converts StateVector into DensityMatrix (possibly reduced one)"""
if trace_over_vibrations:
if isinstance(psi, StateVector):
rho = psi.get_DensityMatrix()
rho = self.trace_over_vibrations(rho)
elif isinstance(psi, StateVectorEvolution):
# FIXME: Implement direct conversion
rho = psi.get_DensityMatrixEvolution()
rho = self.trace_over_vibrations(rho)
else:
if isinstance(psi, StateVector):
rho = psi.get_DensityMatrix()
elif isinstance(psi, StateVectorEvolution):
rho = psi.get_DensityMatrixEvolution()
return rho
[docs]
def get_RWA_suggestion(self) -> float:
"""Returns average transition energy
Average transition energy of the monomer as a suggestion for
RWA frequency
"""
# Nn = self.Nb[1] # number of monomers
esum = 0.0
count = 0
for mn in self.monomers:
for i in range(mn.Nb[1]):
omeg = mn.get_energy(i + 1) - mn.get_energy(0)
esum += omeg
count += 1
return esum / count
def _bath_reorg(self, cfm: Any, indx: int) -> float:
coft = cfm.cfuncs[cfm.get_index_by_where((indx, indx))]
reorg_bath = 0.0
for parm in coft.params:
if parm["ftype"] == "OverdampedBrownian":
reorg_bath += parm["reorg"]
return reorg_bath
def _site_reorg_diag(self, subtract_bath: bool = True) -> numpy.ndarray:
"""Returns the reorganisation energy of an exciton state"""
# SystemBathInteraction
sbi = self.get_SystemBathInteraction()
# CorrelationFunctionMatrix
cfm = sbi.CC
reorg_site = numpy.zeros(self.Ntot)
# electronic states corresponding to single excited states
elst = numpy.where(self.which_band == 1)[0]
for el1 in elst:
if subtract_bath:
reorgB = self._bath_reorg(cfm, el1 - 1)
else:
reorgB = 0.0
reorg = cfm.get_reorganization_energy(el1 - 1, el1 - 1) - reorgB
for kk in self.vibindices[el1]:
reorg_site[kk] += reorg
elst_dbl1 = numpy.where(self.which_band == 2)[0]
for el1 in elst_dbl1:
if subtract_bath:
reorgB = self._bath_reorg(cfm, el1 - 1)
else:
reorgB = 0.0
reorg = cfm.get_reorganization_energy(el1 - 1, el1 - 1) - reorgB
for kk in self.vibindices[el1]:
reorg_site[kk] += reorg
return reorg_site
def _excitonic_reorg_diag(
self, SS: numpy.ndarray, subtract_bath: bool = True
) -> numpy.ndarray:
"""Returns the reorganisation energy of an exciton state"""
# SystemBathInteraction
sbi = self.get_SystemBathInteraction()
# CorrelationFunctionMatrix
cfm = sbi.CC
reorg_exct = numpy.zeros(self.Ntot)
# electronic states corresponding to single excited states
elst = numpy.where(self.which_band == 1)[0]
for n in range(1, self.Nb[1] + 1):
for el1 in elst:
if subtract_bath:
reorgB = self._bath_reorg(cfm, el1 - 1)
else:
reorgB = 0.0
reorg = cfm.get_reorganization_energy(el1 - 1, el1 - 1) - reorgB
for kk in self.vibindices[el1]:
reorg_exct[n] += (SS[kk, n] ** 2) * (SS[kk, n] ** 2) * reorg
elst_dbl1 = numpy.where(self.which_band == 2)[0]
elst_dbl2 = numpy.where(self.which_band == 2)[0]
for n in range(self.Nb[0] + self.Nb[1], self.Ntot):
for el1 in elst_dbl1:
for el2 in elst_dbl2:
if subtract_bath:
reorgB = self._bath_reorg(cfm, el1 - 1)
else:
reorgB = 0.0
reorg = cfm.get_reorganization_energy(el1 - 1, el2 - 1) - reorgB
for kk in self.vibindices[el1]:
for ll in self.vibindices[el2]:
reorg_exct[n] += (SS[kk, n] ** 2) * (SS[ll, n] ** 2) * reorg
return reorg_exct
# FIXME: All these properties have to be meaningfully defined for all open systems
def _get_exciton_prop(self, adiabatic: Any = None, HH_in: Any = None) -> tuple:
is_adiabatic = False
adiabatic_noBath = False
if HH_in is None:
HH = self.HH.copy()
else:
HH = HH_in.copy()
if self._diagonalized:
raise OSError(
"Not possible to obtain the exciton properties for diagonalized aggregate"
)
if adiabatic is not None:
if adiabatic != False:
is_adiabatic = True
else:
is_adiabatic = False
if (adiabatic == "SubtractBath") or (adiabatic == "NoBath"):
adiabatic_noBath = True
if is_adiabatic:
reorg_site = self._site_reorg_diag(subtract_bath=adiabatic_noBath)
for kk in range(self.Ntot):
HH[kk, kk] -= reorg_site[kk]
val, SS = numpy.linalg.eigh(HH.data)
reorg_excit = self._excitonic_reorg_diag(SS, subtract_bath=adiabatic_noBath)
val += reorg_excit
else:
val, SS = numpy.linalg.eigh(HH.data)
return val, SS
def _extract_diagonal_widths(
self, Nel: int, Wd_ini: numpy.ndarray, Dr_ini: numpy.ndarray
) -> tuple[numpy.ndarray, numpy.ndarray]:
Wd_in = numpy.zeros(Nel, dtype=REAL)
Dr_in = numpy.zeros(Nel, dtype=REAL)
for ii in range(Nel):
for k in self.vibindices[ii]:
Wd_in[ii] = Wd_ini[k, k]
Dr_in[ii] = Dr_ini[k, k]
return Wd_in, Dr_in
def _compute_coupling_coefficients(
self, N_states: int, Nel: int, SS: numpy.ndarray, band_filter: int | None = None
) -> numpy.ndarray:
kap = numpy.zeros((N_states, Nel), dtype=REAL)
for aa in range(N_states):
if band_filter is not None:
ela = self.elinds[aa]
if self.which_band[ela] != band_filter:
continue
st = 0
for ii in range(Nel):
for ialph in self.vibindices[ii]:
kap[aa, ii] += numpy.abs(SS[st, aa]) ** 2
st += 1
return kap
def _accumulate_band_linewidths(
self,
N_output: int,
Nel: int,
kap: numpy.ndarray,
Wd_in: numpy.ndarray,
Dr_in: numpy.ndarray,
) -> tuple[numpy.ndarray, numpy.ndarray]:
Wd = numpy.zeros(N_output, dtype=REAL)
Dr = numpy.zeros(N_output, dtype=REAL)
for aa in range(N_output):
for nn in range(Nel):
Wd[aa] += (kap[aa, nn] ** 2) * (Wd_in[nn] ** 2)
Dr[aa] += (kap[aa, nn] ** 2) * (Dr_in[nn] ** 2)
return numpy.sqrt(Wd), numpy.sqrt(Dr)
[docs]
def diagonalize(self) -> None:
"""Transforms some internal quantities into diagonal basis"""
if self._diagonalized:
return
ee, SS = numpy.linalg.eigh(self.HH)
self.Hs = self.HH.copy()
self.HD = ee
self.SS = SS
self.S1 = numpy.linalg.inv(SS)
self.HH = numpy.dot(self.S1, numpy.dot(self.HH, self.SS))
have_vibs = False
if len(self.vibindices[0]) > 1:
have_vibs = True
# Kronecker delta over all states
delta = operator_factory(self.Ntot).unity_operator()
if not have_vibs:
######################################################################
# CASE OF NO VIBRATIONAL MODES
######################################################################
#
# some quantities to be precalculated for two-ex lineshape
# 1->2 has to be trasformed first because we need untransformed 0->1
# for such a transformation
#
N1b = self.Nb[0] + self.Nb[1]
# \kappa_{nA} =
# \sum_{K}(\delta_{nk}+\delta_{nl})*|\langle A | K\rangle|^2
#
# where K is a two-exc. state K = (k,l), A is a two-ex. state
# and n is a single exciton state
#
# Below aa1 = A, aa2 = n, aa3 = K, st_a = k and st_b = l
#
kappa: numpy.ndarray = numpy.zeros((self.Ntot, self.Ntot), dtype=REAL)
if self.mult >= 2:
N2b = self.Nb[0] + self.Nb[1] + self.Nb[2]
# all states (and 2-ex band selected)
for el1 in range(self.Nel):
if self.which_band[el1] == 2:
# all states corresponding to electronic two-exc. state kk
for aa1 in self.vibindices[el1]:
# all states and (1-ex band selected)
for el2 in range(self.Nel):
if self.which_band[el2] == 1:
for aa2 in self.vibindices[el2]:
# all states and (2-ex band selected)
for el3 in range(self.Nel):
if self.which_band[el3] == 2:
for aa3 in self.vibindices[el3]:
st_a = self.twoex_indx[aa3, 0]
st_b = self.twoex_indx[aa3, 1]
if st_b != 0:
kappa[aa2, aa1] += (
delta[aa2, st_a]
+ delta[aa2, st_b]
) * (SS[aa3, aa1] ** 2)
else:
# (transitions on the single molecule (1,0)->(2,0)
e1_trans = (
numpy.nonzero(
self.vibsigs[aa2][0]
)[0]
== numpy.nonzero(
self.vibsigs[aa3][0]
)[0]
).all()
if e1_trans:
kappa[aa2, aa1] += (
SS[aa3, aa1] ** 2
)
#
# Cross terms
#
Wd_tmp = self.Wd.copy()
for aa_2x in range(N1b, N2b):
for alpha in range(N1b):
self.Wd[aa_2x, alpha] = 0.0
for nn_2x in range(N1b, N2b):
for k_1x in range(N1b):
st_c = self.twoex_indx[nn_2x, 0]
st_d = self.twoex_indx[nn_2x, 1]
if st_d != 0:
# FIXME: I think here should be the indexes of the opposite state if st_c == k_1x then the width should be self.Wd[st_d, st_d]
self.Wd[aa_2x, alpha] += (
(
(self.Wd[st_c, st_c] ** 2)
* delta[st_c, k_1x]
+ (self.Wd[st_d, st_d] ** 2)
* delta[st_d, k_1x]
)
* (SS[nn_2x, aa_2x] ** 2)
* (SS[k_1x, alpha] ** 2)
)
else:
e1_trans = (
numpy.nonzero(self.vibsigs[nn_2x][0])[0]
== numpy.nonzero(self.vibsigs[k_1x][0])[0]
).all()
if e1_trans:
self.Wd[aa_2x, alpha] += (
(Wd_tmp[nn_2x, k_1x] ** 2)
* (SS[nn_2x, aa_2x] ** 2)
* (SS[k_1x, alpha] ** 2)
)
self.Wd[N1b:N2b, 0:N1b] = numpy.sqrt(self.Wd[N1b:N2b, 0:N1b])
self.Wd[0:N1b, N1b:N2b] = numpy.transpose(self.Wd[N1b:N2b, 0:N1b])
#
# Transform line shapes for 1->2 transitions
#
Wd_a = numpy.zeros(N2b, dtype=REAL)
Dr_a = numpy.zeros(N2b, dtype=REAL)
for aa in range(N1b, N2b):
for nn in range(N1b, N2b):
st_c = self.twoex_indx[nn, 0]
st_d = self.twoex_indx[nn, 1]
if st_d != 0:
Wd_a[aa] += (SS[nn, aa] ** 2) * (
(self.Wd[st_c, st_c] ** 2) * kappa[st_c, aa]
+ (self.Wd[st_d, st_d] ** 2) * kappa[st_d, aa]
)
else:
Wd_a[aa] += (
(SS[nn, aa] ** 2)
* (Wd_tmp[aa, st_c] ** 2)
* kappa[st_c, aa]
)
W_aux = numpy.diag(numpy.sqrt(Wd_a))
self.Wd[N1b:N2b, N1b:N2b] = W_aux[N1b:N2b, N1b:N2b]
#
# Transform line shapes for 0->1 transitions
#
SS4 = numpy.abs(SS[:N1b, :N1b]) ** 4
Wd_a = numpy.sqrt(
numpy.einsum("in,n->i", SS4, numpy.diag(self.Wd[:N1b, :N1b]) ** 2)
)
Dr_a = numpy.sqrt(
numpy.einsum("in,n->i", SS4, numpy.diag(self.Dr[:N1b, :N1b]) ** 2)
)
self.Wd[0:N1b, 0:N1b] = numpy.diag(Wd_a)
self.Dr[0:N1b, 0:N1b] = numpy.diag(Dr_a)
else:
N1b = self.Nb[0] + self.Nb[1]
# The rest should be uncommented and reviewed for the multilevel systems
# raise IOError("Diagonalization of the transition width with vibronic states and multilevel molecules not yet implemented")
warn_me = False
if warn_me:
# Check for the double excited state on single molecule
for n in range(N1b, self.Ntot):
if self.twoex_indx[n, 1] == 0:
text = (
f"State {n} is double excited state on single molecule."
" For this state transition width transformation is not defined."
" Correct the diagonalization rutine in aggregate_base!!! "
)
warnings.warn(text)
######################################################################
# CASE OF VIBRATIONAL MODES
######################################################################
#
# Transform line shapes for 0->1 transitions
#
Nel = self.Nbe[0] + self.Nbe[1] # GS + singly excited states
Wd_ini = self.Wd.copy()
Dr_ini = self.Dr.copy()
Wd_in, Dr_in = self._extract_diagonal_widths(Nel, Wd_ini, Dr_ini)
kap = self._compute_coupling_coefficients(N1b, Nel, SS, band_filter=1)
Wd_a, Dr_a = self._accumulate_band_linewidths(N1b, Nel, kap, Wd_in, Dr_in)
self.Wd[0:N1b, 0:N1b] = numpy.diag(Wd_a)
self.Dr[0:N1b, 0:N1b] = numpy.diag(Dr_a)
if self.mult >= 2:
Nel = self.Nel
N2b = self.Nb[0] + self.Nb[1] + self.Nb[2]
Wd_in, Dr_in = self._extract_diagonal_widths(Nel, Wd_ini, Dr_ini)
kap2 = self._compute_coupling_coefficients(N2b, Nel, SS)
Wd_b, Dr_b = self._accumulate_band_linewidths(
N1b, Nel, kap2, Wd_in, Dr_in
)
#
# Single exciton band
#
self.Wd[0:N1b, 0:N1b] = numpy.diag(Wd_b)
self.Dr[0:N1b, 0:N1b] = numpy.diag(Dr_b)
Wd_c: numpy.ndarray = numpy.zeros((self.Ntot, self.Ntot), dtype=REAL)
for aa in range(N1b, N2b):
for bb in range(N1b, N2b):
for nn in range(Nel):
vind = self.vibindices[nn]
nni = vind[0]
n = self.twoex_indx[nni, 0]
m = self.twoex_indx[nni, 1]
for mm in range(Nel):
vind = self.vibindices[mm]
mmi = vind[0]
st_a = self.twoex_indx[mmi, 0]
st_b = self.twoex_indx[mmi, 1]
Wd_c[aa, bb] += (
(
(Wd_in[n] ** 2)
* (delta[n, st_a] + delta[n, st_b])
+ (Wd_in[m] ** 2)
* (delta[m, st_a] + delta[m, st_b])
)
* kap2[aa, nn]
* kap2[bb, mm]
)
W_cc: numpy.ndarray = numpy.zeros(Wd_c.shape[0], dtype=REAL)
for diag_idx in range(N2b):
W_cc[diag_idx] = Wd_c[diag_idx, diag_idx]
W_aux = numpy.diag(numpy.sqrt(W_cc))
#
# Two-exciton band
#
self.Wd[N1b:N2b, N1b:N2b] = W_aux[N1b:N2b, N1b:N2b]
Wd_c = numpy.zeros((self.Ntot, self.Ntot), dtype=REAL)
for aa in range(N1b, N2b):
for bb in range(N1b):
for nn in range(Nel):
vind = self.vibindices[nn]
nni = vind[0]
n = self.twoex_indx[nni, 0]
m = self.twoex_indx[nni, 1]
for k in range(1 + self.nmono):
Wd_c[aa, bb] += (
(
(Wd_in[n] ** 2) * delta[n, k]
+ (Wd_in[m] ** 2) * delta[m, k]
)
* kap2[aa, nn]
* kap2[bb, k]
)
W_aux = numpy.sqrt(Wd_c)
self.Wd[N1b:N2b, 0:N1b] = W_aux[N1b:N2b, 0:N1b]
self.Wd[0:N1b, N1b:N2b] = numpy.transpose(self.Wd[N1b:N2b, 0:N1b])
#
#
# Coefficients xi_{ai} to transform pure dephasing of electronic coherence
#
#
# all states and (1-ex band selected)
for aa in range(N1b):
for ii in range(self.Nel):
# el2 = self.elinds[ii]
if self.which_band[ii] == 1:
for ialph in self.vibindices[ii]:
self.Xi[aa, ii] += SS[ialph, aa] ** 2
#
# Transform transition dipole moments
#
for n in range(3):
self.DD[:, :, n] = numpy.dot(self.S1, numpy.dot(self.DD[:, :, n], self.SS))
Ntot = self.HH.shape[0]
dd2 = numpy.sum(self.DD**2, axis=2)
self.D2 = dd2
self.D2_max = numpy.max(dd2)
self.rho0 = numpy.zeros(self.HH.shape, dtype=COMPLEX)
self.rho0[0, 0] = 1.0
self._diagonalized = True
def _thermal_population(
self,
temp: float = 0.0,
subtract: numpy.ndarray | None = None,
relaxation_hamiltonian: numpy.ndarray | None = None,
start: int = 0,
) -> numpy.ndarray:
"""Thermal populations at temperature temp
Thermal populations calculated from the diagonal elements
of the Hamiltonian.
Parameters
----------
temp : float
Temperature in Kelvins
subtract : list like
Reoreganization energies to subtract from the Hamiltonian
relaxation_hamiltonian: array
Hamiltonian according to which we form thermal equilibrium
"""
from ..core.units import kB_intK
kBT = kB_intK * temp
# if not relaxation_hamiltonian:
# HH = self.get_Hamiltonian()
# else:
# HH = relaxation_hamiltonian
HH = relaxation_hamiltonian
# This is all done with arrays, not with Qrhei objects
# HH = HH.data
dim = HH.shape[0]
if subtract is None:
subtract = numpy.zeros(dim, dtype=numpy.float64)
rho0 = numpy.zeros((dim, dim), dtype=numpy.complex128)
if temp == 0.0:
rho0[start, start] = 1.0
else:
# FIXME: we assume only single exciton band
ens = numpy.zeros(dim - start, dtype=numpy.float64)
# we specify the basis from outside. This allows to choose
# canonical equilibrium in arbitrary basis
for i in range(start, dim):
ens[i - start] = numpy.real(HH[i, i] - subtract[i - start])
ne = numpy.exp(-ens / kBT)
sne = numpy.sum(ne)
rho0_diag = ne / sne
rho0[start:, start:] = numpy.diag(rho0_diag)
return rho0
def _thermal_band_population(
self,
temperature: float,
relaxation_theory_limit: str,
relaxation_hamiltonian: Hamiltonian | None,
start: int,
n_states: int,
) -> numpy.ndarray:
if relaxation_theory_limit == "strong_coupling":
if not relaxation_hamiltonian:
HH = self.get_Hamiltonian()
else:
HH = relaxation_hamiltonian
Ndim = HH.dim
re = numpy.zeros(Ndim - start, dtype=numpy.float64)
if not relaxation_hamiltonian:
for i in range(n_states):
re[i] = self.sbi.get_reorganization_energy(i)
ham = HH.data
return self._thermal_population(
temperature, subtract=re, relaxation_hamiltonian=ham, start=start
)
if relaxation_theory_limit == "weak_coupling":
if not relaxation_hamiltonian:
Ham = self.get_Hamiltonian()
else:
Ham = relaxation_hamiltonian
with eigenbasis_of(Ham):
H = Ham.data
subt = numpy.zeros(H.shape[0])
subtfil = numpy.amin(
numpy.array([H[ii, ii] for ii in range(start, H.shape[0])])
)
subt.fill(subtfil)
return self._thermal_population(
temperature, subtract=subt, relaxation_hamiltonian=H, start=start
)
raise Exception("Unknown relaxation_theory_limit")
def _impulsive_population(
self,
relaxation_theory_limit: str = "weak_coupling",
temperature: float = 0.0,
DD: numpy.ndarray | None = None,
) -> numpy.ndarray:
"""Impulsive excitation of the density matrix from ground state"""
rho = self.get_DensityMatrix(
condition_type="thermal",
relaxation_theory_limit=relaxation_theory_limit,
temperature=temperature,
)
rho0 = rho.data
if DD is None:
DD = self.TrDMOp.data
# abs value of the transition dipole moment
dabs = numpy.sqrt(DD[:, :, 0] ** 2 + DD[:, :, 1] ** 2 + DD[:, :, 2] ** 2)
# excitation from bra and ket
rho0 = numpy.dot(dabs, numpy.dot(rho0, dabs))
return rho0
[docs]
def get_StateVector(
self, condition_type: str | None = None, DD: numpy.ndarray | None = None
) -> Any:
"""Returns state vector accordoing to specified conditions
Parameters
----------
condition_type : str
Type of the initial condition. If None, the property sv0, which
was presumably calculated in the past, is returned.
Condition types
---------------
impulsive_excitation
Excitation by ultrabroad laser pulse
"""
# aggregate must be built before we call this method
if not self._built:
raise QuantarheiError(
"Aggregate must be built before get_StateVector can be invoked."
)
# if no condition is specified, it is understood that we return
# internal sv0, which was calculated sometime in the past
if condition_type is None:
return StateVector(data=self.sv0)
if condition_type == "impulsive_excitation":
if DD is None:
DD = self.TrDMOp.data
# abs value of the transition dipole moment
dabs = numpy.sqrt(DD[:, :, 0] ** 2 + DD[:, :, 1] ** 2 + DD[:, :, 2] ** 2)
sv0: numpy.ndarray = numpy.zeros(self.Ntot, dtype=COMPLEX)
# zero temperature
sv0[0] = 1.0
self.sv0 = numpy.dot(dabs, sv0)
return StateVector(data=self.sv0)
raise QuantarheiError("Unknown condition type")
[docs]
def get_DensityMatrix(
self,
condition_type: str | None = None,
relaxation_theory_limit: str = "weak_coupling",
temperature: float | None = None,
relaxation_hamiltonian: Hamiltonian | None = None,
DD: numpy.ndarray | None = None,
) -> DensityMatrix:
"""Returns density matrix according to specified condition
Returs density matrix to be used e.g. as initial condition for
propagation.
Parameters
----------
condition_type : str
Type of the initial condition. If None, the property rho0, which
was presumably calculated in the past, is returned.
relaxation_theory_limits : str {weak_coupling, strong_coupling}
Type of the relaxation theory limits;
We mean the system bath coupling. When `weak_coupling` is chosen,
the density matrix is returned in form of a canonical equilibrium
in terms of the exciton basis. For `strong_coupling`,
the canonical equilibrium is calculated in site basis with site
energies striped of reorganization energies.
temperature : float
Temperature in Kelvin
relaxation_hamiltonian :
Hamiltonian according to which we form thermal equilibrium. In case
of `strong_coupling`, no reorganization energies are subtracted -
we assume that the supplied energies are already void of them.
DD : real array
Submit your own transition dipole moment matrix (for x, y and z)
Condition types
---------------
thermal
Thermally equilibriated population of the whole density matrix
thermal_excited_state
Thermally equilibriuated excited state
impulsive_excitation
Excitation by ultrabroad laser pulse
"""
# aggregate must be built before we call this method
if not self._built:
raise QuantarheiError(
"Aggregate must be built before get_DensityMatrix can be invoked."
)
# if Aggregate has interaction with the bath, temperature
# is already defined
if temperature is None:
if self.sbi is None:
temperature = 0.0
elif self.sbi.has_temperature():
temperature = self.sbi.get_temperature()
else:
temperature = 0.0
# if no condition is specified, it is understood that we return
# internal rho0, which was calculated sometime in the past
if condition_type is None:
return DensityMatrix(data=self.rho0)
# impulsive excitation from a thermal ground state
if condition_type == "impulsive_excitation":
rho0 = self._impulsive_population(
relaxation_theory_limit=relaxation_theory_limit,
temperature=temperature,
DD=DD,
)
self.rho0 = rho0
return DensityMatrix(data=self.rho0)
# thermal population based on the total Hamiltonian
if condition_type == "thermal":
if not relaxation_hamiltonian:
Ham = self.get_Hamiltonian()
else:
Ham = relaxation_hamiltonian
# FIXME: weak and strong limits not distinguished
rho0 = self._thermal_population(
temperature, relaxation_hamiltonian=Ham.data
)
self.rho0 = rho0
return DensityMatrix(data=self.rho0)
if condition_type == "thermal_excited_state":
rho0 = self._thermal_band_population(
temperature,
relaxation_theory_limit,
relaxation_hamiltonian,
start=self.Nb[0],
n_states=self.Nb[1],
)
self.rho0 = rho0
return DensityMatrix(data=self.rho0)
if condition_type == "thermal_twoexciton":
rho0 = self._thermal_band_population(
temperature,
relaxation_theory_limit,
relaxation_hamiltonian,
start=self.Nb[1] + self.Nb[0],
n_states=self.Nb[2],
)
self.rho0 = rho0
return DensityMatrix(data=self.rho0)
raise QuantarheiError("Unknown condition type")
#
# TESTED
[docs]
def get_temperature(self) -> float:
"""Returns temperature associated with this aggregate
The temperature originates from the system-bath interaction
"""
# aggregate must be built before we call this method
if not self._built:
raise QuantarheiError()
if self.sbi is None:
return 0.0
return self.sbi.CC.get_temperature()
[docs]
def get_electronic_groundstate(self) -> tuple:
"""Indices of states in electronic ground state
Returns indices of all states in the electronic
ground state of the system.
"""
return tuple(range(self.Nb[0]))
[docs]
def get_excitonic_band(self, band: int = 1) -> tuple:
"""Indices of states in a given excitonic band.
Returns indices of all states in the excitonic band
with number of excitons equal to `band`
Parameters
----------
band : int
Specifies which band should be returned.
"""
return self.get_band(band=band)
[docs]
def get_transition(self, Nf: Any, Ni: Any) -> tuple:
"""Returns relevant info about the energetic transition
Parameters
----------
Nf : {int, ElectronicState, VibronicState}
Final state of the transition
Ni : {int, ElectronicState VibronicState}
Initial state of the transition
"""
if isinstance(Nf, ElectronicState) and isinstance(Ni, ElectronicState):
if self.Ntot == self.Nel:
iNf = Nf.index
iNi = Ni.index
else:
raise QuantarheiError("The Hamiltonian is not pure electronic")
elif isinstance(Nf, VibronicState) and isinstance(Ni, VibronicState):
vsig = Nf.get_vibsignature()
esig = Nf.get_ElectronicState().get_signature()
iNf = self.vibsigs.index((esig, vsig))
vsig = Ni.get_vibsignature()
esig = Ni.get_ElectronicState().get_signature()
iNi = self.vibsigs.index((esig, vsig))
else:
iNf = Nf
iNi = Ni
#
# if Nf and Ni are not of the same type, it will lead to Exception
#
energy = self.convert_energy_2_current_u(self.HH[iNf, iNf] - self.HH[iNi, iNi])
trdipm = self.DD[iNf, iNi, :]
return (energy, trdipm)
[docs]
def has_SystemBathInteraction(self) -> bool:
"""Returns True if the Aggregate is embedded in a defined environment"""
# aggregate must be built before we call this method
if not self._built:
raise QuantarheiError("Aggregate must be built before this call.")
if (self.sbi is not None) and self._has_system_bath_interaction:
return True
return False
[docs]
def get_SystemBathInteraction(self) -> Any:
"""Returns the aggregate SystemBathInteraction object"""
if self._built:
return self.sbi
raise BuildError("Aggregate object not built.")
[docs]
def set_SystemBathInteraction(self, sbi: Any) -> None:
"""Sets the SystemBathInteraction operator for this aggregate"""
# FIXME: check its compatibility
self.sbi = sbi
self.sbi.set_system(self)
[docs]
def get_Hamiltonian(self) -> Any:
"""Returns the aggregate Hamiltonian"""
if self._built:
return self.HamOp # Hamiltonian(data=self.HH)
raise BuildError("Aggregate object not built")
[docs]
def get_electronic_Hamiltonian(self, full: bool = False) -> Any:
"""Returns the aggregate electronic Hamiltonian
In case this is a purely electronic aggregate, the output
is identical to get_Hamiltonian()
"""
HH: numpy.ndarray = numpy.zeros((self.Nel, self.Nel), dtype=REAL)
for a, sta in self.elstates(mult=self.mult):
HH[a, a] = sta.energy()
for b, stb in self.elstates(mult=self.mult):
if a != b:
HH[a, b] = self.coupling(sta, stb, full=full)
HHel = Hamiltonian(data=HH)
return HHel
[docs]
def get_TransitionDipoleMoment(self) -> Any:
"""Returns the aggregate transition dipole moment operator"""
if self._built:
return self.TrDMOp # TransitionDipoleMoment(data=self.DD)
raise BuildError("Aggregate object not built")
[docs]
def get_TransitionMagneticDipoleMoment(self) -> Any:
"""Returns the aggregate transition dipole moment operator"""
if self._built:
return TransitionDipoleMoment(data=self.MM)
raise BuildError("Aggregate object not built")