Quantum Mechanics with Python
Numerical Methods for Physicists, Lecture 4
Mat Medo, Yi-Cheng Zhang
Physics Department, Fribourg University, Switzerland
March 11, 2013
Medo & Zhang (UNIFR)
Quantum Python
1 / 19
Schrdinger equation (here in 1D)
i~
(x, t)
~2 2 (x, t)
+ V (x)(x, t)
=
t
2m
x 2
PDE with complex numbers
scipy.integrate.complex_ode re-maps a given
complex-valued differential equation to a set of real-valued ones
and uses scipy.integrate.ode to integrate them
But we have partial derivatives. . .
Medo & Zhang (UNIFR)
Quantum Python
2 / 19
Schrdinger equation (here in 1D)
i~
(x, t)
~2 2 (x, t)
+ V (x)(x, t)
=
t
2m
x 2
PDE with complex numbers
scipy.integrate.complex_ode re-maps a given
complex-valued differential equation to a set of real-valued ones
and uses scipy.integrate.ode to integrate them
But we have partial derivatives. . .
1
Do re-mapping to real values and use Finite-Difference
Time-Domain (FDTD) as we did for the heat equation
Use the split-step Fourier method
QuTiP: a quantum toolbox in Python
Re-mapping & Escript (fast finite elements for PDEs)
Medo & Zhang (UNIFR)
Quantum Python
2 / 19
FDTD: re-mapping
Re-mapping to real-valued variables:
(x, t) = (x, t) + i(x, t)
Schrdinger equation decouples:
(x, t)
~2 2 (x, t)
+ V (x)(x, t)
=
t
2m x 2
~2 2 (x, t)
(x, t)
=
~
V (x)(x, t)
t
2m x 2
Discretization in time and space: xl = l x, tn = n t
(xl , tn ) := n (l),
Medo & Zhang (UNIFR)
(xt , tn ) := n (l)
Quantum Python
3 / 19
FDTD: discretization
Real and complex part of shifted by half-interval from each
otherthis allows us to use the second-order approximation for
time derivatives which is more precise (O(t) vs O(t 2 ))
(xl , tn+1/2 )
(xl , tn+1 ) (xl , tn )
n+1 (l) n (l)
=
t
t
t
(x
,
t
)
(x
,
t
)
(xl , tn )
n+1/2 (l) n1/2 (l)
l n+1/2
l n1/2
=
t
t
t
Spatial derivatives exactly as before
2 n (l)
n (l + 1) 2n (l) + n (l 1)
x 2
x 2
Medo & Zhang (UNIFR)
Quantum Python
4 / 19
FDTD: equations
Compute the future state from the current state
{n , n1/2 } {n+1 , n+1/2 }
Update equations are
n+1/2 (l) = n1/2 (l) c2 V (l)n (l) +
+ c1 n (l 1) 2n (l) + n (l + 1)
n+1 (l) = n (l) + c2 V (l) n+1/2 (l)
c1 n+1/2 (l 1) 2 n+1/2 (l) + n+1/2 (l + 1)
Where c1 := ~t/(2mx 2 ) and c2 := t/~
Medo & Zhang (UNIFR)
Quantum Python
5 / 19
FDTD: practical issues 1
It is easier to work in units where ~ = m = 1
2
The scheme is stable only if t ~ m ~x 2 +
maxx,t V (x,t) 1
2
See http://www.scipy.org/Cookbook/SchrodingerFDTD
(there is the rest of theory as well)
This is only orientational (e.g., V (x) somewhere does not
necessarily mean t 0)
Finer grid = smaller time steps needed
To reach the same physical time, simulation time grows as 1/x 3
Check whether we are really in the stable (convergent) regime:
halve t and see whether the results change substantially (they
should not)
We are typically interested in the probability density , not in
and themselves
Medo & Zhang (UNIFR)
Quantum Python
6 / 19
FDTD: practical issues 2
Its nice to have a reversible algorithm (as in the leap-frog kind of
integration schemes)
Evolve your system to the future and back and get the initial state
Evolved wave-function needs to stay normalized
Force normalization by hand after each step (possible but
non-systematic)
Find a unitary algorithm which preserves the norm
See Numerical recipes in C (Section 19.2) by Press et al
We approximated the time evolution operator U(t) = exp[iHt] as
U(t) 1 iHt whereas it is better to use
U(t) = exp[iHt]
1 21 iHt
1 + 12 iHt
Then naturally U(t)U(t) = 1
Medo & Zhang (UNIFR)
Quantum Python
7 / 19
Detour: animations with Python
Sometimes better than static figures
Especially valuable if you want to impress. . .
First approach:
Write a script that saves many figures, say frame-001.gif,
frame-002.gif, etc.
Use avconv (formerly ffmpeg) to create an animation
avconv -an -r 4 -i frame-%03d.png animation.avi
Or use gifsicle to create an animated gif
gifsicle -delay=100 -loop -no-transparent -O2
list_of_files > output.gif
Second approach: animate with Matplotlib (v1.1 and higher)
http://jakevdp.github.com/blog/2012/08/18/
matplotlib-animation-tutorial/
Medo & Zhang (UNIFR)
Quantum Python
8 / 19
Split-step Fourier method: the idea
The potential part of the Schrdinger equation is simple to solve
i~
pot
= V (x) = (x, t + t) = (x, t)eiV (x)t/~
t
Inverse Fourier transform of
1
(x, t) =
2
, t)eikx dk
(k
By substituting this into the Schrdinger equation
i~
~2 k 2
= (k
, t + t) kin
, t)ei~k 2 t/2m
=
+ iV
= (k
t
2m
k
Fast Fourier transform (FFT) allows us to switch between and
Medo & Zhang (UNIFR)
Quantum Python
9 / 19
Split-step Fourier method: the algorithm
Discretize the x and k -space by choosing a, b, N, k0 :
x = (b a)/N,
k = 2/(b a),
xn = a + n x,
k0 = /x,
km = k0 + m k
This should be sufficient to represent states of (initial and later)
The choice can be tricky, experiment a bit to see if it works fine
2
Discretize the wave-functions:
n (t) := (xn , t),
Vn := V (xn ),
Evolve the system by one time step
Repeat 3 until the end time is reached
Medo & Zhang (UNIFR)
Quantum Python
m (t) = (k
m , t)
10 / 19
Split-step Fourier method: point 3
A half step in x: n n exp[iVn t/2~]
m from n
FFT: computes
m
m exp[i~k 2 t/2m]
A full step in k :
m
Inverse FFT: computes n from
A second half step in x: n n exp[iVn t/2~]
Splitting the x-step in two parts produces a numerically more
stable algorithm (as in leapfrog integration in mechanics which is
reversible and conserves energy)
See http://jakevdp.github.com/blog/2012/09/05/
quantum-python/ for more details
Medo & Zhang (UNIFR)
Quantum Python
11 / 19
Split-step Fourier method: tunelling
Output of ex4-1.py (see also the animation at the url above):
1.0
norm = 1.000
1.0
t=4
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0 150
1.0
100
50
50
norm = 1.000
100
150
0.0 150
1.0
t=8
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0 150
100
Medo & Zhang (UNIFR)
50
50
100
150
norm = 1.000
100
50
t=6
50
norm = 1.000
0.0 150
Quantum Python
100
50
100
150
t = 10
50
100
150
12 / 19
Fast Fourier transform
Discrete Fourier transform is rather slow to compute: O(N 2 )
where N is the number of points where the transform is computed
The beauty of FFT: O(N 2 ) is reduced to O(N log N) (similar as
quicksort improving over bubble sort and alike)
The idea behind: discrete Fourier transform with N points can be
written as a sum of two transforms with N/2 points
Using this recursively, we move to N/4 points, etc.
Thats why its best to use FFT for N which is a power of two
fft and ifft from scipy.fftpack implement it
FFT works best when is a power of two (because of how it is
constructed)
N 6= 2x is possible but usually leads to slower computation
Medo & Zhang (UNIFR)
Quantum Python
13 / 19
QuTiP & Escript
QuTiP:
QM is not as difficult as one might think: Its only linear algebra!
States, operators, time evolution, gates, visualization
http://qutip.blogspot.ch and
http://code.google.com/p/qutip
Escript:
For non-linear, time-dependent partial differential equations
https://launchpad.net/escript-finley
Medo & Zhang (UNIFR)
Quantum Python
14 / 19
What about the stationary states?
Until now we addressed only time evolution in QM
Stationary state n satisfies
~2 2
Hn = En n
+ V (r) n = En n
2m
Boundary conditions in a box: (a) = (b) = 0
We fix (a) = 0 and find E so that (b) is zero too
In open space: limx (x) = limx (x) = 0
For a symmetric potential, we simplify by realizing that eigenstates
are also symmetric
Even eigenstates: 0 (0) = 0; odd eigenstates: (0) = 0
We need E for which (x) does go to zero as x grows
Medo & Zhang (UNIFR)
Quantum Python
15 / 19
Stationary states: integration
Again (xn ) := n , V (xn ) = Vn
In addition, ~ = m = 1 and h = x
The approximation for f 00 (x) plus the Schrdinger equation
1
00n + Vn n = En = 00n = 2(Vn E)n
2
n1 2n + n+1
00n =
= n+1 = 2n n1 + 2(Vn E)h2
h2
Once 0 and 1 are given, the rest can be computed
Medo & Zhang (UNIFR)
Quantum Python
16 / 19
Stationary states: integration
Again (xn ) := n , V (xn ) = Vn
In addition, ~ = m = 1 and h = x
The approximation for f 00 (x) plus the Schrdinger equation
1
00n + Vn n = En = 00n = 2(Vn E)n
2
n1 2n + n+1
00n =
= n+1 = 2n n1 + 2(Vn E)h2
h2
Once 0 and 1 are given, the rest can be computed
Much better (order of 4 instead of 2 above) is Numerovs method
2
h2
(2 5h6 fn )yn (1 + 12
fn1 )yn1
d2
+ f (x) y (x) = 0 = yn+1 =
2
2
h
dx
1 + 12 fn+1
Medo & Zhang (UNIFR)
Quantum Python
16 / 19
Stationary states: continuation
How to find the right E: the shooting method
Medo & Zhang (UNIFR)
Quantum Python
17 / 19
Stationary states: continuation
How to find the right E: the shooting method
It is best to write a function, say psi_end(E), which takes E as
an argument and returns (x) at some (not too) distant point x
Finding E efficiently: e.g., brentq(f,a,b) from
scipy.optimize finds a root of f between a and b
Medo & Zhang (UNIFR)
Quantum Python
17 / 19
Stationary states: LHO
We first plot (xm ) for a range of energies
40
even
odd
(4)
20
0
20
400
Medo & Zhang (UNIFR)
Quantum Python
18 / 19
Stationary states: LHO
We first plot (xm ) for a range of energies
40
even
odd
(4)
20
0
20
400
Then individual roots can be found
ex4-1.py gives 0.501 and 1.500 for the first two eigenstates
Medo & Zhang (UNIFR)
Quantum Python
18 / 19
Try this at home
1: Reversibility and unitarity
Implement the Finite-Difference Time-Domain algorithm. Is it reversible? Is it
unitary?
2: Ground state
p
Find the ground state energy of potential V (x) = exp( |x|) in the units
where ~ = m = 1. Can you find also the first excited state?
Medo & Zhang (UNIFR)
Quantum Python
19 / 19
