Test: Marcus Model

For this test, we use the class HarmonicOscillator from section 6.2.1.

mass, omega = 1.0, 1.0
c = 1.0
d = 1.0/(2.0*6.0**2)

ho = HarmonicOscillator(mass, omega)
xmin = -20.0
xmax = -xmin
ngrid = 256
displacement = 3.0
e_ge = 0.0
x = np.linspace(xmin, xmax, ngrid)
nstate = 2

v11 = ho.get_potential(x - displacement)
v22 = e_ge + ho.get_potential(x + displacement)
v12 = c * np.exp(-d * x**2)

vgrid = np.zeros((x.shape[0], nstate, nstate)) # Dimension: ngrid * nstate * nstate
psi0 = np.zeros((x.shape[0], nstate))          # Dimension: ngrid * nstate
vgrid[:, 0, 0] = v11
vgrid[:, 1, 1] = v22
vgrid[:, 0, 1] = v12
vgrid[:, 1, 0] = v12
psi00 = ho.get_eigenfunction(x + displacement, 0)
psi0[:, 1] = psi00
psi0[:, 0] = np.zeros(x.shape[0])

mu = np.zeros((x.shape[0], nstate, nstate))
mu_eg = 1.0
mu[:, 0, 1] = mu_eg * np.ones(x.shape[0])
mu[:, 1, 0] = mu_eg * np.ones(x.shape[0])

wp_multin = QuantumDynamics1D(mass)
wp_multin.setup_grid(xmin, xmax, ngrid)
wp_multin.set_potential(vgrid)
wp_multin.set_dipole_coupling(mu)

n = 200
nstep = 2000
tstep = (2*np.pi/omega) / n

def electric_field(t, e0, omega, sigma, td):
    return e0 * np.cos(omega*t) * np.exp(-(t - td)**2/(2.0*sigma**2))

e0 = 5.0
omega_eg = e_ge + ho.get_potential(-2 * displacement)
sigma_t = 20 * tstep
td = 300 * tstep
wp_multin.set_efield(lambda t: electric_field(t, e0, omega_eg, sigma_t, td))

wp_multin.propagate(tstep, nstep, psi0, method="TD Split Operator")
anim_multin = wp_multin.animate(
        nstep=nstep, delay=10, scalewf=50.0, plot_potential=True,
)