# 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(-My*np.pi/LyT,My*np.pi/LyT-2*np.pi/LyT,My)

k0z = np.linspace(-Mz*np.pi/LzT,Mz*np.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,My*My),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.j*kIz*z0),(My,1))

temp_y = np.tile(np.exp(-((y0-y0c)/σ_y)**2)*np.exp(1.j*kIy*y0),(Mz,1)).T

temp_1 = temp_y*temp

ψ_0[iS,:] = temp_1.reshape(1,My*Mz)##Gaussian wavepacket

##
ψ_0[iS,:] = ψ_0[iS,:]/np.sqrt(ψ_0[iS,:]@ψ_0[iS,:].T)##Normalized wavefunction

#Now we need to calculate the propagators

##Kinetic T = exp(i*hbark^2.2mdt/hbar)

TP = np.exp(-1.j*(np.tile((k0y**2)/2,(Mz,1)).T+np.tile((k0z**2)/2,(My,1)))*dt)

T =  np.tile((k0y*k0y)/2,(Mz,1)).T + np.tile((k0z*k0z)/2,(My,1))

ψ_0

##Potential Energy propagator V

ψ = ψ_0 #Initialization of wavepacket

tR = np.linspace(0,dt*Tsteps,Tsteps)
tR
ek = np.zeros((Tsteps, N), dtype = 'complex')

e = np.zeros((Tsteps, 1), dtype = 'complex')

ev = np.zeros((Tsteps,1), dtype = 'complex')

norm = np.zeros((Tsteps,1), dtype = 'complex')

meabs = np.zeros((Tsteps,N), dtype ='complex').T
#print(psi) Uncheck to make sure everything is proper
for t in range(Tsteps):
ψ
##T propgators

for n in range(N):
temp2 = ψ[n,:]
temp3 = temp2.reshape(My,Mz)
temp4 = fft.fftshift(fft.fft2(temp3))
ek[t,n] = np.real(np.sum(np.sum(np.conj(temp4)*T*temp4))/np.sum(np.sum(np.conj(temp4)*temp4)))
temp5 = temp4*TP#Kinetic Propagator
temp6 = fft.ifft2(fft.fftshift(temp5))
ψ[n,:] = temp6.reshape(1,My*Mz)

####Now we check the energy conservation####
###ek(t)+ev(t) = constant

meabs[:,t] = np.sum(np.conj(ψ)*ψ,1) ## population on each dibat
e[t] =  ek[t,0]*meabs[0,t]
norm[t] = np.real(np.sum(np.sum(np.conj(ψ)*ψ_0)))
ψ_0 = ψ

fig = plt.figure(figsize = (5,5))
plt.plot(np.real(meabs[0,:]))
plt.xlabel('Time')
plt.ylabel(' $$S_{0}$$ Population')
plt.show()

plt.xlabel('Time')
plt.ylabel(' $$S_{0_{adi}}$$ Population')
plt.show()

plt.imshow(np.abs(ψ[0,:].reshape(Mz,Mz)))
plt.show()
plt.plot(np.abs(e))
plt.xlabel('Time')
plt.ylabel('Energy')
plt.show()

plt.plot(np.abs(norm))
plt.xlabel('Time')
plt.ylabel('Norm')

#g.s ψ

fig = plt.figure(figsize = (15,15))
X,Y = np.meshgrid(y0,z0)
ax = plt.axes(projection='3d')
surf =ax.plot_surface(X,Y,np.abs(ψ[0,:].reshape(My,Mz)),cmap = cm.inferno)
ax.view_init(25, 50)
plt.xlabel('x')
plt.ylabel('y')
fig.colorbar(surf, shrink=0.5, aspect=5)


Everything looks good.

But when I try to extend my model to 3D, I am quite confused on how I should set up my momenta/position surfaces, and as such, how does this effect my Guassian wavepackets and kinetic propagators. Here I added the transpose to the 2D $$y-component$$ and to the $$x-component$$ as well in the $$3D$$ case. but I do not have a basis for doing so. Below is my code for the 3D:

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

##Split-Operator ###
# Constants
ω = np.array([1,1,1]) #Frequency for each coordinate
gs = np.array([0,0,0]) #Inital Wavepacket shifts
iS = 0 # Intial starting state
y0c = gs[1]

z0c = gs[2]

x0c = gs[0]

##intial Momenta

kIz = 0

kIy = 0

kIx = 0

##Set up a grid

My = 128*2

Mz = 128*2

Mx = 128*2

#Number of states

N = 1

#Number of time steps

Tsteps = 200

dt = 0.005

# Grid Lengths

Ly = 10

Lz = 10

Lx = 10

LyT = Ly*2

LzT = Lz*2

LxT = Lz*2

#Grid of M points

y0 = np.linspace(-Ly,Ly,My)

z0 = np.linspace(-Lz,Lz, Mx)

x0 = np.linspace(-Lx,Lx, Mz)

#Parameters
#k0[1xM] = Grid of M momenta points from 0->L

k0y = np.linspace(-My*np.pi/LyT,My*np.pi/LyT-2*np.pi/LyT,My)

k0z = np.linspace(-Mz*np.pi/LzT,Mz*np.pi/LzT-2*np.pi/LzT,Mz)

k0x = np.linspace(-Mx*np.pi/LxT,Mx*np.pi/LxT-2*np.pi/LxT,Mx)

##Properties

##Postion and momenta surfaces##

y0op = (np.tile(y0,(Mz,1))).T

z0op = np.tile(z0,(Mx,1))

x0op = (np.tile(x0,(My,1))).T

k0yop = (np.tile(k0y,(My,1))).T

k0zop = (np.tile(k0z,(Mx,1)))

k0xop = (np.tile(k0x,(My,1))).T

##Inital wavefunction: 2D gaussian
ψ_0 = np.zeros((N,My*My),dtype = 'complex')
σ_y = np.sqrt(2/ω[1])
σ_z = np.sqrt(2/ω[2])
σ_x = np.sqrt(2/ω[0])

temp = np.tile(np.exp(-((z0-z0c)/(σ_z))**2)*np.exp(1.j*kIz*z0),(Mx,1))

temp_y = np.tile(np.exp(-((y0-y0c)/σ_y)**2)*np.exp(1.j*kIy*y0),(Mz,1)).T

temp_1 = temp_y*temp*np.tile(np.exp(-((x0-x0c)/σ_x)**2)*np.exp(1.j*kIx*x0),(My,1)).T

ψ_0[iS,:] = temp_1.reshape(1,My*Mz)##Gaussian wavepacket

##
ψ_0[iS,:] = ψ_0[iS,:]/np.sqrt(ψ_0[iS,:]@ψ_0[iS,:].T)##Normalized wavefunction

#Now we need to calculate the propagators

##Kinetic T = exp(i*hbark^2.2mdt/hbar)

TP = np.exp(-1.j*(np.tile((k0y**2)/2,(Mz,1)).T+np.tile((k0x**2)/2,(My,1)).T+np.tile((k0z**2)/2,(Mx,1)))*dt)

T =  np.tile((k0y*k0y)/2,(Mz,1)).T +np.tile((k0x*k0x)/2,(My,1)).T + np.tile((k0z*k0z)/2,(Mx,1))

ψ = ψ_0 #Initialization of wavepacket

tR = np.linspace(0,dt*Tsteps,Tsteps)

ek = np.zeros((Tsteps, N), dtype = 'complex')

e = np.zeros((Tsteps, 1), dtype = 'complex')

ev = np.zeros((Tsteps,1), dtype = 'complex')

norm = np.zeros((Tsteps,1), dtype = 'complex')

meabs = np.zeros((Tsteps,N), dtype ='complex').T
#print(psi) Uncheck to make sure everything is proper
for t in range(Tsteps):
ψ
##T propgators

for n in range(N):
temp2 = ψ[n,:]
temp3 = temp2.reshape(My,Mz)
temp4 = fft.fftshift(fft.fft2(temp3))
ek[t,n] = np.real(np.sum(np.sum(np.conj(temp4)*T*temp4))/np.sum(np.sum(np.conj(temp4)*temp4)))
temp5 = temp4*TP    #Kinetic Propagator
temp6 = fft.ifft2(fft.fftshift(temp5))
ψ[n,:] = temp6.reshape(1,My*Mz)

####Now we check the energy conservation####
###ek(t)+ev(t) = constant

meabs[:,t] = np.sum(np.conj(ψ)*ψ,1) ## population on each dibat
e[t] =  ek[t,0]*meabs[0,t]
norm[t] = np.real(np.sum(np.sum(np.conj(ψ)*ψ_0)))
ψ_0 = ψ

fig = plt.figure(figsize = (5,5))
plt.plot(np.real(meabs[0,:]))
plt.xlabel('Time')
plt.ylabel(' $$S_{0}$$ Population')
plt.show()

###Visualization###
plt.plot(np.abs(e))
plt.xlabel('Time')
plt.ylabel('Energy')
plt.show()

plt.plot(np.abs(norm))
plt.xlabel('Time')
plt.ylabel('Norm')

#g.s ψ

fig = plt.figure(figsize = (15,15))
X,Y = np.meshgrid(y0,z0)
ax = plt.axes(projection='3d')
surf =ax.plot_surface(X,Y,np.abs(ψ[0,:].reshape(My,Mz)),cmap = cm.inferno)
ax.view_init(25, 50)
plt.xlabel('x')
plt.ylabel('y')
fig.colorbar(surf, shrink=0.5, aspect=5)



If anyone could help guide me through this problem, or explain to me the correction it would help greatly.

Thanks, and if further information is required please let me know.

• I think that the topic of your question is related to this site. Your question... I'm not sure. Right now it looks like you want us to understand what you're trying to achieve from your code Nov 11 '21 at 2:38
• Hi @Nicoguaro thanks for your comment. I've updated the code to include the x-coordinate I am still not sure about the transpose of coordinate within the TP function. I guess I do not fully understand the split-operator method. Nov 11 '21 at 14:16
• I'd suggest to improve the description of the problem and add the math if possible. Nov 11 '21 at 16:13
• Having the continuous equations and your approximation scheme before the code would make the question clearer. Also, removing most of the blank lines from the code would greatly help with its readability in this format. Nov 11 '21 at 18:32
• Thanks @BillBarth and nicogruaro. I'll edit it up and come back. Nov 11 '21 at 20:53