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.

Share this:

Like this: