Source code for quantarhei.qm.liouvillespace.foerstertensor

from __future__ import annotations

import numpy

from ... import COMPLEX
from ...core.managers import energy_units
from ...exceptions import QuantarheiError

# import scipy.interpolate as interp
from ..hilbertspace.hamiltonian import Hamiltonian
from ..liouvillespace.systembathinteraction import SystemBathInteraction
from .rates.foersterrates import FoersterRateMatrix
from .relaxationtensor import RelaxationTensor




class FoersterRelaxationTensor(RelaxationTensor):
    """Weak resonance coupling relaxation tensor by Foerster theory"""

    dim: int
    data: numpy.ndarray
    _has_cutoff_time: bool
    cutoff_time: float | None
    _is_initialized: bool

    def __init__(
        self,
        ham: Hamiltonian,
        sbi: SystemBathInteraction,
        initialize: bool = True,
        cutoff_time: float | None = None,
        pure_dephasing: bool = False,
    ) -> None:

        self._initialize_basis()

        if not isinstance(ham, Hamiltonian):
            raise QuantarheiError("First argument must be a Hamiltonian")

        if not isinstance(sbi, SystemBathInteraction):
            raise QuantarheiError(
                "Second argument must be of type SystemBathInteraction"
            )

        self._is_initialized = False
        self._has_cutoff_time = False
        self.as_operators = False
        self.pure_dephasing = pure_dephasing
        self.has_Iterm = False

        if cutoff_time is not None:
            self.cutoff_time = cutoff_time
            self._has_cutoff_time = True

        self.Hamiltonian = ham
        self.dim = ham.dim
        self.SystemBathInteraction: SystemBathInteraction = sbi
        self.TimeAxis = sbi.TimeAxis

        if initialize:
            self.initialize()
            self._data_initialized = True
            self._is_initialized = True

        else:
            self._data_initialized = False



    def initialize(self) -> None:

        #
        # Tensor data
        #
        Na = self.dim
        self.data = numpy.zeros((Na, Na, Na, Na), dtype=COMPLEX)

        with energy_units("int"):
            HH = self.Hamiltonian
            sbi = self.SystemBathInteraction

            if self._has_cutoff_time:
                cft = self.cutoff_time
            else:
                cft = None

            frm = FoersterRateMatrix(HH, sbi, initialize=True, cutoff_time=cft)

            #
            # Transfer rates
            #
            for aa in range(self.dim):
                for bb in range(self.dim):
                    if aa != bb:
                        self.data[aa, aa, bb, bb] = frm.data[aa, bb]

            #
            # calculate dephasing rates and depopulation rates
            #
            self.updateStructure()

            if self.pure_dephasing:
                # additional pure dephasing
                self.add_dephasing()




    def add_dephasing(self) -> None:

        # line shape function derivatives
        sbi = self.SystemBathInteraction
        Na = self.dim
        Nt = sbi.TimeAxis.length

        ht = numpy.zeros((Na, Nt), dtype=numpy.complex64)

        sbi.CC.create_one_integral()

        for ii in range(1, Na):
            ht[ii, :] = sbi.CC.get_hoft(ii - 1, ii - 1)

        for aa in range(self.dim):
            for bb in range(self.dim):
                if aa != bb:
                    self.data[aa, bb, aa, bb] -= ht[aa, Nt - 1] + ht[bb, Nt - 1]