Split Operator 2D-->3D
I am trying to modify my simulations on population dynamics using the split-operator method from 2D two-level real Hamiltonain $$H_{2D} =T(x,y)\otimes1_2 +\begin{pmatrix} -z & y\ y & z\ \end{pmatrix},$$ to a complex 3D two level Hamiltonian $$H_{3D} =T(x,y,z)\otimes1_2 + \begin{pmatrix} -z & x+iy\ x-iy & z\ \end{pmatrix}.$$ To simplify things for myself, I've thrown away the potential contributions and I am just focusing on the kinetic terms ensuring that norm/energy conservation is good. Loosely speaking, the split operator method follows as $$e^{-iH\Delta t} \approx e^{-iT\Delta t}e^{-iV\Delta t} +\mathcal{O}(\Delta t)^2$$ and for our case $e^{-iH\Delta t} \approx e^{-iT\Delta t}$, with no error as it can be shown that $\text{Error } =[T,V] = 0$. In the 2D case: ``` import numpy as np import scipy as sp import matplotlib.pyplot as plt import scipy.special as scl import numpy.matlib as mat import scipy.fftpack as fft from mpl_toolkits.mplot3d import Axes3D from matplotlib import cm from matplotlib.widgets import Slider, Button
Split-Operator
Constants
ω = np.array([1,1]) #Frequency for each coordinate
gs = np.array([0,0]) #Inital Wavepacket shifts
g = 0.0825 # g-vector:(g_j-g_i)/2
h = 0.0430
s_y = 0# s_y vector (g_jy+g_iy)/2
s_z = 0.125 # s_y vector (g_jz+g_iz)/2
Location of conical intersection: (0,0,0)--> you can translate to coordinates in article if desired
If you want to switch what surface WP is propagating on, change here
iS = 0 # Intial starting state
y0c = gs[1]
z0c = gs[0]
intial Momenta
kIz = 0
kIy = 0
Set up a grid #Lets see how this goes
My = 128*2
Mz = 128*2
Number of states
N = 1
Number of time steps
Tsteps = 600
dt = 0.005
Grid Lengths
Ly = 10
Lz = 10
LyT = Ly*2
LzT = Lz*2
Grid of M points
y0 = np.linspace(-Ly,Ly,My)
z0 = np.linspace(-Lz,Lz, Mz)
Parameters
k0[1xM] = Grid of M momenta points from 0->L
k0y = np.linspace(-Mynp.pi/LyT,Mynp.pi/LyT-2*np.pi/LyT,My)
k0z = np.linspace(-Mznp.pi/LzT,Mznp.pi/LzT-2*np.pi/LzT,Mz)
Properties
Postion and momenta grids
y0op = (np.tile(y0,(Mz,1))).T
z0op = np.tile(z0,(My,1))
k0yop = (np.tile(k0y,(My,1))).T
k0zop = (np.tile(k0z,(Mz,1)))
Inital wavefunction: 2D gaussian
ψ_0 = np.zeros((N,MyMy),dtype = 'complex') σ_x = np.sqrt(2/ω[0]) σ_y = np.sqrt(2/ω[1]) σ_z = np.sqrt(2/ω[1]) temp = np.tile(np.exp(-((z0-z0c)/(σ_z))2)np.exp(1.jkIzz0),(My,1))
temp_y = np.tile(np. 402