# Wavefunctions and eigenenergies for a quantum Morse oscillator¶

Robert Johansson ([email protected])

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib

In [2]:
from scipy import *
from scipy import optimize

In [3]:
from wavefunction import *
from wavefunction.wavefunction1d import *


## Problem parameters¶

In [4]:
h = 6.626e-34
h = 1
h_ = h/(2*pi)
e = 1.602e-19
cf = 1          # if cf = h_, use units where h_ = 1

mm = 1          # oscillator mass
omega = 2 * pi  # oscillator frequency in GHz
x0 = 0          # shift in oscillator potiential minimum

D = 5.0
b = 0.5

args = {'D': D, 'b': b}

k = -h_ ** 2 / (2 * mm)

x_min = -pi
x_max =  pi
N = 750


## Potential¶

In [5]:
def U_morse(x, args):
"""
Morse oscillator potential
"""

D = args['D']
b = args['b']

u = D * (1 - exp(-b*x)) ** 2

return u


### Find the minimum point for the potential¶

In [7]:
x = linspace(x_min, x_max, N)

In [15]:
U = U_morse(x, args);

x_opt_min = optimize.fmin(U_morse, [0.0], (args,))

Optimization terminated successfully.
Current function value: 0.000000
Iterations: 3
Function evaluations: 6

In [16]:
print("Found potential minima at = %f" % x_opt_min)

Found potential minima at = 0.000000


### Plot the potential¶

In [17]:
fig, ax = subplots()

ax.plot(x, U)
ax.plot(x_opt_min, U_morse(x_opt_min, args), 'o')

#ax.set_ylim(-10, 80)
ax.set_xlabel(r'$x$', fontsize=18)
ax.set_ylabel(r'$U(x)$', fontsize=18);


## Eigenstates and eigenvalues¶

In [19]:
u = assemble_u_potential(N, U_morse, x, args)

K = assemble_K(N, k, x_min, x_max)

V = assemble_V(N, u, x_min, x_max)

H = K + V

evals, evecs = solve_eigenproblem(H)

In [21]:
NN = 10

fig, axes = subplots(NN, 1, figsize=(10,NN*1), sharex=True, sharey=True)

for n in range(NN):
Y = evecs[NN-n-1]
axes[n].plot(x, Y.real)

axes[n].set_xlabel(r'$x$', fontsize=18);

In [33]:
fig, ax = subplots(figsize=(12,8))

ax.plot(x, U, 'k')
for n in range(10):
Y = evals[n] + evecs[n]

ax.set_xlim(-2, 4)
ax.set_ylim(0, 4)
ax.set_xlabel(r'$x$', fontsize=18)
ax.set_ylabel(r'$U(x)$', fontsize=18);


## Identify bound states¶

In [35]:
nbound = 0
boundidx = []

for i in range(0,N):
if nbound < 10:
if inner(evecs[i], evecs[i]) > 0.85:
nbound = nbound + 1
boundidx.append(i)

print "Found bound states: ", boundidx

Found bound states:  [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]


## Evaluate transition matrix elements¶

In [36]:
inner_prod = zeros((len(boundidx), len(boundidx))).astype(np.complex)
expect_pos = zeros((len(boundidx), len(boundidx))).astype(np.complex)
expect_kin = zeros((len(boundidx), len(boundidx))).astype(np.complex)

for i, l in enumerate(boundidx):
for j, k in enumerate(boundidx):

psi_l = wavefunction_normalize(x, evecs[l])
psi_k = wavefunction_normalize(x, evecs[k])

inner_prod[i,j] = inner_product(x, psi_l, psi_k)                 # <psi_l|psi_k>
expect_pos[i,j] = inner_product(x, psi_l, x * psi_k)             # <psi_l|x|psi_k>
expect_kin[i,j] = inner_product(x, psi_l, derivative(x, psi_k))  # <psi_l|p|psi_k>

In [37]:
print "inner_prod:\n"
print_matrix(inner_prod)

inner_prod:

1.000    0.000   -0.000    0.000   -0.000    0.000   -0.000   -0.000    0.000   -0.000
0.000    1.000   -0.000    0.000   -0.000    0.000   -0.000   -0.000    0.000   -0.000
-0.000   -0.000    1.000   -0.000    0.000    0.000    0.000   -0.000    0.000   -0.000
0.000    0.000   -0.000    1.000   -0.000    0.000    0.000   -0.000    0.000   -0.000
-0.000   -0.000    0.000   -0.000    1.000   -0.000    0.000   -0.000    0.000   -0.000
0.000    0.000    0.000    0.000   -0.000    1.000   -0.000    0.000   -0.000    0.000
-0.000   -0.000    0.000    0.000    0.000   -0.000    1.000   -0.000    0.000   -0.000
-0.000   -0.000   -0.000   -0.000   -0.000    0.000   -0.000    1.000   -0.000   -0.000
0.000    0.000    0.000    0.000    0.000   -0.000    0.000   -0.000    1.000   -0.000
-0.000   -0.000   -0.000   -0.000   -0.000    0.000   -0.000   -0.000   -0.000    1.000

In [34]:
print "expect_pos ="
print2DMatrix(expect_pos)

expect_pos =
0.038    0.226   -0.018   -0.002    0.000    0.000   -0.000    0.000    0.000    0.000
0.226    0.116    0.321    0.032   -0.005   -0.001    0.000   -0.000   -0.000   -0.000
-0.018    0.321    0.196   -0.396    0.046    0.008   -0.002    0.000    0.000    0.000
-0.002    0.032   -0.396    0.279    0.460    0.060   -0.011    0.003    0.001    0.000
0.000   -0.005    0.046    0.460    0.364   -0.518    0.074   -0.015   -0.004   -0.001
0.000   -0.001    0.008    0.060   -0.518    0.453    0.571   -0.089   -0.020   -0.005
-0.000    0.000   -0.002   -0.011    0.074    0.571    0.544    0.621    0.104    0.025
0.000   -0.000    0.000    0.003   -0.015   -0.089    0.621    0.639   -0.668   -0.120
0.000   -0.000    0.000    0.001   -0.004   -0.020    0.104   -0.668    0.739   -0.713
0.000   -0.000    0.000    0.000   -0.001   -0.005    0.025   -0.120   -0.713    0.843

In [39]:
print "expect_kin:\n"
print_matrix(expect_kin)

expect_kin:

0.000    2.184   -0.346    0.068   -0.016    0.004   -0.001    0.000   -0.000    0.000
-2.184   -0.000    3.028   -0.592    0.135   -0.035    0.010   -0.003    0.001   -0.000
0.346   -3.028    0.000    3.634   -0.825    0.211   -0.060    0.019   -0.006    0.002
-0.068    0.592   -3.634   -0.000    4.108   -1.049    0.296   -0.091    0.030   -0.011
0.016   -0.135    0.825   -4.108   -0.000    4.494   -1.265    0.388   -0.129    0.046
-0.004    0.035   -0.211    1.049   -4.494   -0.000    4.813   -1.472    0.486   -0.172
0.001   -0.010    0.060   -0.296    1.265   -4.813   -0.000    5.079   -1.671    0.589
-0.000    0.003   -0.019    0.091   -0.388    1.472   -5.079   -0.000    5.300   -1.862
0.000   -0.001    0.006   -0.030    0.129   -0.486    1.671   -5.300   -0.000    5.482
-0.000    0.000   -0.002    0.011   -0.046    0.172   -0.589    1.862   -5.482   -0.000

In [41]:
print "eigenenergies = "
print evals[0:nbound].real

eigenenergies =
[ 0.12519193  0.37080402  0.61004865  0.84292835  1.06944556  1.28960261
1.5034017   1.71084493  1.91193429  2.10667171]


### Software¶

In [42]:
%reload_ext version_information
%version_information numpy, scipy, matplotlib, wavefunction

Out[42]:
SoftwareVersion
Python2.7.4 (default, Apr 19 2013, 18:28:01) [GCC 4.7.3]
IPython2.0.0-dev
OSposix [linux2]
numpy1.8.0.dev-895866d
scipy0.13.0.dev-7c6d92e
matplotlib1.3.x
wavefunction1.0.0
Mon Sep 02 15:16:24 2013 JST
In [ ]: