Robert Johansson (robert@riken.jp)
%pylab inline
Populating the interactive namespace from numpy and matplotlib
from scipy import *
from scipy import optimize
from wavefunction import *
from wavefunction.wavefunction1d import *
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
def U_morse(x, args):
"""
Morse oscillator potential
"""
D = args['D']
b = args['b']
u = D * (1 - exp(-b*x)) ** 2
return u
x = linspace(x_min, x_max, N)
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
print("Found potential minima at = %f" % x_opt_min)
Found potential minima at = 0.000000
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);
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)
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);
fig, ax = subplots(figsize=(12,8))
ax.plot(x, U, 'k')
for n in range(10):
Y = evals[n] + evecs[n]
mask = where(Y > U)
ax.plot(x[mask], evals[n] * ones(shape(x))[mask], 'k--')
mask = where(Y > U-2.0)
ax.plot(x[mask], Y[mask].real)
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);
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]
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>
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
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
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
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]
%reload_ext version_information
%version_information numpy, scipy, matplotlib, wavefunction
Software | Version |
---|---|
Python | 2.7.4 (default, Apr 19 2013, 18:28:01) [GCC 4.7.3] |
IPython | 2.0.0-dev |
OS | posix [linux2] |
numpy | 1.8.0.dev-895866d |
scipy | 0.13.0.dev-7c6d92e |
matplotlib | 1.3.x |
wavefunction | 1.0.0 |
Mon Sep 02 15:16:24 2013 JST |