2

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.exp(-((y0-y0c)/σ_y)*2)np.exp(1.jkIyy0),(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((k0y2)/2,(Mz,1)).T+np.tile((k0z2)/2,(My,1)))dt)


T =  np.tile((k0yk0y)/2,(Mz,1)).T + np.tile((k0zk0z)/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(-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)


k0x = np.linspace(-Mxnp.pi/LxT,Mxnp.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.jkIzz0),(Mx,1))


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


temp_1 = temp_ytempnp.tile(np.exp(-((x0-x0c)/σ_x)*2)np.exp(1.jkIxx0),(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((k0y2)/2,(Mz,1)).T+np.tile((k0x2)/2,(My,1)).T+np.tile((k0z2)/2,(Mx,1)))dt)


T =  np.tile((k0yk0y)/2,(Mz,1)).T +np.tile((k0xk0x)/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.

New2Python
  • 21
  • 3
  • 1
    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 – nicoguaro Nov 11 '21 at 02: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. – New2Python Nov 11 '21 at 14:16
  • 1
    I'd suggest to improve the description of the problem and add the math if possible. – nicoguaro Nov 11 '21 at 16:13
  • 1
    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. – Bill Barth Nov 11 '21 at 18:32
  • Thanks @BillBarth and nicogruaro. I'll edit it up and come back. – New2Python Nov 11 '21 at 20:53

0 Answers0