In [1]:
%matplotlib inline
from pylab import *
from numpy import *

In [2]:
#physical constants
hbar = 1.0
m = 1.0

In [3]:
#discretization
nx = 100

minx = -1.
maxx = 1.

x = linspace(minx, maxx, nx)

In [4]:
#assemble the matrix Hamiltonian
from d2matrix import d2matrix
from scipy.sparse import dia_matrix

def make_H(V):
nx = len(V)
kinetic = -hbar**2/(2*m)*d2matrix(nx)
potential = diag(V)
H = kinetic + potential
return H

V = zeros(nx)

H = make_H(V)

/usr/local/lib/python3.4/dist-packages/numpy/core/fromnumeric.py:2507: VisibleDeprecationWarning: rank is deprecated; use the ndim attribute or function instead. To find the rank of a matrix see numpy.linalg.matrix_rank.
VisibleDeprecationWarning)

In [5]:
# calculate the eigenvectors
eigvals, eigvecs = linalg.eig(H)
eigvecs = array(eigvecs)

In [6]:
for i in arange(nx):
plot(x, eigvecs[:,i]+eigvals[i]*150)
ylim([0, 10])

Out[6]:
(0, 10)
In [7]:
def solve(H):
eigvals, eigvecs = linalg.eig(H)

I = eigvals.argsort()
eigvals = eigvals[I]
eigvecs = array(eigvecs)[:,I]
return eigvals, eigvecs

In [8]:
V = x**2
H = make_H(V)
eigvals, eigvecs = solve(H)

In [9]:
figure(figsize(5,6))
for i in arange(10):
plot(x, eigvecs[:,i]*.05+eigvals[i])
plot(x, V, lw=3)
ylim([0, 0.3])

Out[9]:
(0, 0.3)

Two atoms:

In [10]:
#discretization
nx = 100

minx = -1.
maxx = 1.

x = linspace(minx, maxx, nx)

V = zeros(nx)
a1 = -0.5
a2 = 0.5
asize = 0.2
apot = -0.07
V[(x>a1-asize) & (x<a1+asize) | (x>a2-asize) & (x<a2+asize) ] = apot

In [11]:
eigvals, eigvecs = solve(make_H(V))

In [12]:
figure(figsize(6,6))
for i in arange(6):
plot(x, eigvecs[:,i]*.05+eigvals[i])
plot(x, V, lw=2)
set_printoptions(precision=16)
print(eigvals[:6])

#ylim([-1, -0.6])

[-0.0624387366052539 -0.062438578859082  -0.0406528132169814
-0.0406493136442565 -0.0091801384671852 -0.0089414977975879]


Evolution of a localized electron? Evolution of each eigenvector follows from Schrodinger equation $$i\hbar\frac{\partial\psi_j}{\partial t} = \hat H\psi_j = E_j\psi_j$$ $$\psi_j(t) = e^{{E_j t}/{i \hbar}}\psi_j(0)$$ And for general initial state $$\psi(0) = \sum_j c_j\psi_j(0), {\rm where }\quad c_j = \langle\psi_j(0)|\psi(0)\rangle$$ is given by $$\psi(t) = \sum_j c_j\psi_j(t) = \sum_j c_j e^{{E_j t}/{i \hbar}}\psi_j(0)$$

In [15]:
coeffs = zeros(nx)
coeffs[0] = 1/sqrt(2.)
coeffs[1] = 1/sqrt(2.)
E_mean = dot(coeffs*eigvals, coeffs)

In [16]:
def psi_t(eigvecs, eigvals, coeffs, t):
phase = exp(eigvals*t/1j*hbar)
#print phase

psi = zeros(len(coeffs), dtype=complex)
for i in arange(len(coeffs)):
if coeffs[i] != 0.:
psi += eigvecs[:,i]*phase[i]
return psi

In [17]:
from matplotlib import animation
from JSAnimation.IPython_display import display_animation

In [18]:
fig = figure()
ax = axes(xlim=(minx,maxx))
line1, = ax.plot([],[],lw=1)
ax.plot(x, V)

def animate(data):
line1.set_data(x, abs(psi_t(eigvecs, eigvals, coeffs, 2000000.*data))*.1+E_mean)

return line1,

anim = animation.FuncAnimation(fig, animate, frames=arange(20), interval=100)
display_animation(anim, default_mode='loop')

Out[18]:

Once Loop Reflect
In [ ]: