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 = 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
args = {'w': omega, 'm': mm, 'x0': x0}
k = -h_ ** 2 / (2 * mm)
x_min = -pi
x_max = pi
N = 750
def U_ho(x, args):
"""
Harmonic oscillator potential
"""
omega = args['w']
m = args['m']
x0 = args['x0']
u = 1/2.0 * m * (omega ** 2) * ((x - x0) ** 2)
return u
x = linspace(x_min, x_max, N)
U = U_ho(x, args);
x_opt_min = optimize.fmin(U_ho, [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_ho(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_ho, 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_xlim(-1, 1)
axes[n].set_xlabel(r'$x$', fontsize=18)
#axes[n].set_ylabel(r'$\Psi_n(x)$', fontsize=18);
fig.tight_layout();
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].real * ones(shape(x))[mask], 'k--')
mask = where(Y > U-2.0)
ax.plot(x[mask], Y[mask].real)
ax.set_xlim(-1, 1)
ax.set_ylim(0, 10)
ax.set_xlabel(r'$x$', fontsize=18)
ax.set_ylabel(r'$U(x)$', fontsize=18);
boundidx = arange(10)
nbound = len(boundidx)
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 "eigenenergies = "
evals[0:nbound].real
eigenenergies =
array([ 0.50058072, 1.50156845, 2.50220863, 3.50250109, 4.50244565, 5.50204213, 6.50129033, 7.50019009, 8.49874122, 9.49694352])
print("Harmonicity:")
((evals[boundidx] - evals[boundidx[0]]) / (evals[boundidx[1]]-evals[boundidx[0]])).real
Harmonicity:
array([ 0. , 1. , 1.9996528 , 2.99895823, 3.99791609, 4.99652621, 5.99478841, 6.99270251, 7.99026831, 8.98748564])
print "orthogonality:\n"
print_matrix(inner_prod)
orthogonality: 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:\n"
print_matrix(expect_pos)
expect_pos: -0.000 -0.113 -0.000 -0.000 -0.000 -0.000 0.000 -0.000 -0.000 -0.000 -0.113 -0.000 -0.159 -0.000 -0.000 0.000 -0.000 -0.000 -0.000 0.000 -0.000 -0.159 0.000 -0.195 -0.000 -0.000 -0.000 -0.000 0.000 -0.000 -0.000 -0.000 -0.195 -0.000 -0.225 0.000 -0.000 0.000 -0.000 -0.000 -0.000 -0.000 -0.000 -0.225 -0.000 -0.252 -0.000 -0.000 -0.000 -0.000 -0.000 0.000 -0.000 0.000 -0.252 0.000 -0.276 -0.000 -0.000 0.000 0.000 -0.000 -0.000 -0.000 -0.000 -0.276 0.000 -0.298 0.000 -0.000 -0.000 -0.000 -0.000 0.000 -0.000 -0.000 -0.298 -0.000 -0.318 -0.000 -0.000 -0.000 0.000 -0.000 -0.000 -0.000 0.000 -0.318 0.000 -0.337 -0.000 0.000 -0.000 -0.000 -0.000 0.000 -0.000 -0.000 -0.337 -0.000
print "expect_kin:\n"
print_matrix(expect_kin)
expect_kin: 0.000 -4.438 -0.000 -0.003 -0.000 -0.000 0.000 -0.000 -0.000 -0.000 4.438 -0.000 -6.272 -0.000 -0.006 0.000 -0.000 -0.000 -0.000 0.000 0.000 6.272 0.000 -7.678 0.000 -0.009 -0.000 -0.000 0.000 -0.000 0.003 0.000 7.678 -0.000 -8.861 -0.000 -0.013 0.000 -0.000 -0.000 0.000 0.006 -0.000 8.861 0.000 -9.902 -0.000 -0.017 -0.000 -0.000 0.000 -0.000 0.009 0.000 9.902 -0.000 -10.842 -0.000 -0.021 0.000 -0.000 0.000 0.000 0.013 0.000 10.842 0.000 -11.704 0.000 -0.026 0.000 0.000 0.000 -0.000 0.017 0.000 11.704 -0.000 -12.506 -0.000 0.000 0.000 -0.000 0.000 0.000 0.021 -0.000 12.506 0.000 -13.257 0.000 -0.000 0.000 0.000 0.000 -0.000 0.026 0.000 13.257 0.000