Quantum Dynamics

We solve the Schrodinger equation numerically for a Gaussian wavepacket in Morse and harmonic potential. We define our Gaussian wavepacket at t=0 as,

\psi(x,t=0)=\Big(\frac{m\omega}{\pi \hbar}\Big)^{1/4}\mathrm{e}^{- \frac{m \omega}{2 \hbar}(x-x_0)^2}.

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  \Delta t/2 in position space by solving,

 \mathrm{i}\hbar\frac{\partial \psi(x,t)}{\partial t}=V(x) \psi(x,t) ,

where V(x) is the potential energy surface. The solution to equation above is given by

 \psi(x,t+\Delta t/2)=\mathrm{e}^{-\frac{\mathrm{i} V(x)}{\hbar} \frac{\Delta t}{2}} \psi(x,t)

Use \texttt{numpy} module \texttt{fft.fft} to obtain wavefunction numerically in momentum space. To obtain the p vector, use \texttt{fft.fftfreq}. If \psi is a grid of n samples in real space, \texttt{fft.fft}(\psi) gives you a grid of n numbers in Fourier space, such that the first element corresponds to the zero-frequency mode. It is helpful to checkout the documentation of \texttt{scipy} on Fast Fourier Transform (FFT) to learn about the ordering of the FFT output. With this knowledge, propagate the next full-time step \Delta t in Fourier space by solving,

\mathrm{i}\hbar\frac{\partial \tilde{\psi}(p,t)}{\partial t} =\frac{p^2}{2 m } \tilde{\psi}(p,t),

where \tilde{\psi}(p,t) is simply the wavefunction in Fourier space. The solution to the equation above is given by,

\tilde{\psi}(p,t+\Delta t)=\mathrm{e}^{-\frac{\mathrm{i} p^2 }{2 m \hbar} \Delta t} \tilde{\psi}(p,t)

Finally use \texttt{numpy} module \texttt{fft.ifft} to obtain your wavefunction back in real space, and repeat the propagation in this space to complete the final \Delta t/2 propagation. Therefore, in one iteration of the algorithm, the wave packet propagates by \Delta t as follows,

\psi(x,t+ \Delta t) =\mathrm{e}^{-\frac{\mathrm{i} V(x)}{\hbar} \frac{\Delta t}{2} }\mathcal{F}^{-1}\Big(\mathrm{e}^{-\frac{\mathrm{i} p^2 }{2 m \hbar} \Delta t}\mathcal{F}\big(\mathrm{e}^{-\frac{\mathrm{i} V(x)}{\hbar} \frac{\Delta t}{2} }\psi(x,t)\big)\Big).

For small enough \Delta t, 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 \psi(x,t=0), we implement the split operator algorithm for the following Morse potential,

V(x)=D\big(1-\mathrm{e}^{-\alpha (x-x_e)}\big)^2.

morsenew

Using \texttt{ffmpeg} 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,

V(x)=\frac{1}{2} m \omega^2 x^2.

We initialize our wave packet at t=0 with a Gaussian, \psi(x,t=0). What does the time-dynamics of the Gaussian wavepacket in a harmonic potential tell you about your initial state \psi(x,t=0)? How do you explain this interesting observation below?

harmonicnew

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.