We solve the Schrodinger equation numerically for a Gaussian wavepacket in Morse and harmonic potential. We define our Gaussian wavepacket at as,
.
One cannot generally obtain the eigenstates of the quantum propagator. We shall therefore numerically simulate the time-dynamics of a wave packet moving in Morse potential using the split operator method in python. In this method, one propagates a wave packet alternately in real space and Fourier space. We first propagate by time-step in position space by solving,
,
where is the potential energy surface. The solution to equation above is given by
Use module
to obtain wavefunction numerically in momentum space. To obtain the
vector, use
. If
is a grid of
samples in real space,
(
) gives you a grid of
numbers in Fourier space, such that the first element corresponds to the zero-frequency mode. It is helpful to checkout the documentation of
on Fast Fourier Transform (FFT) to learn about the ordering of the FFT output. With this knowledge, propagate the next full-time step
in Fourier space by solving,
,
where is simply the wavefunction in Fourier space. The solution to the equation above is given by,
Finally use module
to obtain your wavefunction back in real space, and repeat the propagation in this space to complete the final
propagation. Therefore, in one iteration of the algorithm, the wave packet propagates by
as follows,
.
For small enough , and large number of iterations, the approximate propagator obtained from this method is identical to the exact one. Starting from the initial Gaussian wave packet
, we implement the split operator algorithm for the following Morse potential,
.
Using program, one can combine several of the snapshots of the wavepackets into a movie such as the following,
Now we perform a wavepacket propagation using the same code in a harmonic potential defined by,
.
We initialize our wave packet at with a Gaussian,
. What does the time-dynamics of the Gaussian wavepacket in a harmonic potential tell you about your initial state
? How do you explain this interesting observation below?
The time-dynamics of the Gaussian wavepacket in a harmonic oscillator can be seen in the movie below,
Run the code below to perform simulation of a Gaussian wavepacket moving in a two-dimensional Muller-Brown potential energy surface:
</pre> <pre>""" Author: Manish J. Thapa Date: 10.04.2020 ETH Zurich, Switzerland (Please keep this if you want to modify or use this code!) """ from __future__ import division, print_function import numpy as np from matplotlib import pyplot as plt from scipy.fftpack import fft2, ifft2 class propagation(object): def __init__(self, xarr, yarr,gauss2d, Vxy, m, hbar=1, dt=0.01 ): self.psi_2d = gauss2d #initial wavepacket self.xarr = xarr self.yarr=yarr self.Vxy = Vxy #potential energy surface self.hbar = hbar self.m = m self.dt=dt self.dx = self.xarr[1] - self.xarr[0] self.dy = self.yarr[1] - self.yarr[0] def evolution(self,steps): self.kx = 2 * np.pi * np.fft.fftfreq(len(self.xarr), d=self.dx) #obtain the k vector self.ky = 2 * np.pi * np.fft.fftfreq(len(self.yarr), d=self.dy) psi_x = self.psi_2d for i in range(steps): kx2=self.kx * self.kx ky2=self.ky * self.ky kx2ky2 = np.zeros((len(self.xarr), len(self.yarr))) for ixid, ix in enumerate(kx2): for iyid, iy in enumerate(ky2): kx2ky2[ixid, iyid] = ix+iy halflin = np.exp(-0.5 * 1j * self.Vxy / self.hbar * self.dt) fullnonlin=np.exp(-0.5 * 1j * self.hbar /self.m * (kx2ky2) * self.dt - 0.5 * 1j * self.hbar /self.m * (kx2ky2) * self.dt) psi_x = psi_x * halflin #half-step linear propagation in x-space psi_k = fft2(psi_x) psi_k = psi_k*fullnonlin #full-step non-linear propagation in fourier space psi_x = ifft2(psi_k) psi_x = psi_x * halflin #half-step linear propagation in x-space self.psi_2d=psi_x #set some global variables m=100 omega=1.1 def gauss2d(x,y, x0gauss,y0gauss): #initialize a wavepacket kickx = 0 #give a momentum kick in x direction kicky = 0 #give a momentum kick in y direction return np.exp(-m*omega/2*(x-x0gauss)**2)*np.exp(-m*omega/2*(y-y0gauss)**2)*np.exp(1j*(x-x0gauss)*kickx)*np.exp(1j*(y-y0gauss)*kicky) #normalize as appropriate A = [-200, -100, -170, 15] a = [-1, -1, -6.5, 0.7] b = [0, 0, 11, 0.6] c = [-10, -10, -6.5, 0.7] x0 = [1, 0, -0.5, -1] y0 = [0, 0.5, 1.5, 1] #grid for the surface (PES) Vx=np.linspace(-1.5,1,1e2) Vy=np.linspace(-0.5,1,1e2) def potentialMB(x,y): #define PES V = 0 for i in range(4): V += A[i] * np.exp(a[i] *(x-x0[i])**2+ b[i]*(x-x0[i])*(y-y0[i]) + c[i]*(y-y0[i])**2) return V Vxy=np.zeros((len(Vx),len(Vy))) for xid,xval in enumerate(Vx): for yid,yval in enumerate(Vy): Vxy[xid,yid]=potentialMB(xval,yval) #some critical parameters to centralize your wavepacket at if 0: #TS x0gauss=-0.82200249 y0gauss=0.62431232 elif 0: #min 1 x0gauss = -0.05001082 y0gauss = 0.4666941 elif 0: #min 2 x0gauss = -0.558 y0gauss = 1.442 elif 0: x0gauss = -0.81 y0gauss = 0.95 elif 1: x0gauss = -0.327 y0gauss = 0.579 else: x0gauss = -0.6 y0gauss = 0.5 #grid for your wavepacket to spread across (keep it same as for the PES) x=Vx y=Vy initGauss=np.zeros((len(x),len(y))) for xid,xval in enumerate(x): #initialize the wavepacket at t=0 for yid,yval in enumerate(y): initGauss[yid,xid]=gauss2d(xval,yval, x0gauss,y0gauss) #check transpose here but I believe it is correct Z0=np.real(abs(initGauss)**2) #probabilities #make an object of your propagation class evolve = propagation(xarr=x,yarr=y,gauss2d=initGauss,Vxy=Vxy.T,hbar=1,m=m,dt=0.001) #check transpose here but I believe it is correct (note that norm is preserved regardless of the value of dt) from matplotlib import animation fig = plt.figure(figsize=(8,8)) ax = fig.gca() ax.set_xlim(min(Vx), max(Vx)) ax.set_ylim(min(Vy), max(Vy)) ax.plot(-0.82200249,0.62431232,'-ro') ax.plot(-0.05001082, 0.4666941, '-bo') ax.plot(-0.558, 1.442, '-bo') ax.plot(0.61, 0.03, '-bo') ax.set_xlabel('X') ax.set_ylabel('Y') list=np.linspace(-500,1000,40) #for the contour lines def movie(i): evolve.evolution(i) Z = abs(evolve.psi_2d) ** 2 ax.contourf(Vx, Vy, Vxy.T, cmap='cividis',levels=list) ax.contour(Vx, Vy, Vxy.T, cmap='cividis',levels=list) ax.contourf(x, y, Z, cmap='Greys', alpha=1.0) return ax, anim = animation.FuncAnimation(fig, movie,interval=3,frames=100) #anim.save('MB.mp4', fps=3, extra_args=['-vcodec', 'libx264']) plt.show()</pre> <pre>
Note: I prepared these simulation results for Advanced Kinetics class taught by Prof. Jeremy O. Richardson at ETH Zurich.