拆分运算符 2D-->3D

计算科学 Python
2021-12-10 09:49:55

我正在尝试使用来自 2D 两级真实哈密顿的分裂算子方法修改我对种群动力学的模拟

H2D=T(x,y)⊗12+(−zyyz),
到复杂的 3D 两级哈密顿量
H3D=T(x,y,z)⊗12+(−zx+iyx−iyz).
为了简化我自己的事情，我放弃了潜在的贡献，我只关注动力学术语，确保规范/能量守恒是好的。松散地说，拆分运算符方法如下
e−iHΔt≈e−iTΔte−iVΔt+O(Δt)2
对于我们的情况e−iHΔt≈e−iTΔt，没有错误，因为可以证明Error =[T,V]=0. 在二维情况下：

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)

一切看起来都很好。

但是当我尝试将我的模型扩展到 3D 时，我对如何设置我的动量/位置表面非常困惑，因此，这对我的高斯波包和动力学传播器有何影响。在这里，我将转置添加到 2Dy−component 并到x−component以及在3D案子。但我没有这样做的依据。以下是我的 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)

如果有人可以帮助指导我解决这个问题，或者向我解释更正，那将有很大帮助。

谢谢，如果需要更多信息，请告诉我。

0个回答
没有发现任何回复~
其它你可能感兴趣的问题
上一篇ICCG负残留产品r吨米− 1rrTM−1r 下一篇LSMR：更新∥ x ∥‖x‖便宜