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 [ ]: