TENSO

from tenso.prototypes.bath import gen_bcf

bath = gen_bcf(
    re_d    = [540],    # λ for DL component (cm⁻¹)
    width_d = [70],     # ωc for DL component (cm⁻¹)
    freq_b  = [1663],   # ω0 for BO component (cm⁻¹)
    re_b    = [330],    # λ  for BO component (cm⁻¹)
    width_b = [4],      # η  for BO component (cm⁻¹)
    temperature          = 300,    # K
    decomposition_method = 'Pade',
    n_ltc                = 1,      # low-temperature corrections
)

import numpy as np

eps   = 1500.0   # cm⁻¹  energy gap ε
delta =  300.0   # cm⁻¹  tunneling Δ

sys_ham = np.array([[ eps/2,  delta],
                    [ delta, -eps/2]], dtype=np.complex128)

sys_op  = np.array([[0.0, 1.0],
                    [1.0, 0.0]], dtype=np.complex128)

# Initial state: |↑⟩ (excited state)
wfn      = np.array([1.0, 0.0], dtype=np.complex128)
init_rdo = np.outer(wfn, wfn.conj())

from math import ceil
from tqdm import tqdm
from tenso.prototypes.heom import system_multibath

end_time  = 100.0   # fs
dt = 0.05    # fs

propagator = system_multibath(
    fname             = 'out',
    init_rdo          = init_rdo,
    sys_ham           = sys_ham,
    sys_ops           = [sys_op],
    bath_correlations = [bath],
    dim               = 20,
    end_time          = end_time,
    step_time         = dt,
)

progress_bar = tqdm(propagator, total=ceil(end_time / dt))
for _t in progress_bar:
    progress_bar.set_description(f'@{_t:.2f} fs')

import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt('out.dat.log', dtype=np.complex128)

t       = data[:, 0].real    # time (fs)
pop_exc = data[:, 1].real    # excited state population
pop_gnd = data[:, 4].real    # ground state population
coh_eg  = data[:, 2]         # coherence ρ_eg (complex)

fig, axes = plt.subplots(1, 2, figsize=(8, 3))
axes[0].plot(t, pop_exc, label=r'$[\rho_S]_{00}$')
axes[0].plot(t, pop_gnd, '--', label=r'$[\rho_S]_{11}$')
axes[0].set(xlabel='Time (fs)', ylabel='Population')
axes[0].legend()
axes[1].plot(t, np.abs(coh_eg), color='tab:purple')
axes[1].set(xlabel='Time (fs)', ylabel=r'$|\rho_{eg}|$')
plt.tight_layout()
plt.show()