We are solving a 1D Schrödinger equation (in atomic units): $$ \left(-{1\over 2\mu} {d^2\over dx^2} + V(x)\right)\psi(x) = E\psi(x) $$ For example for double well: $$ V(x) = D_0 {(x^2 - a^2)^2\over a^4} $$ we want to get:
We follow the Appendix 2 in [1], which explains how to use sine basis to solve the problem.
[1] Horáček, J., Gemperle, F., Meyer, H.-D. (1996). Calculation of dissociative attachment of electrons to diatomic molecules by the Schwinger–Lanczos approach. The Journal of Chemical Physics, 104(21), 8433. doi:10.1063/1.471593
In [1]:
%pylab inline
from sympy.interactive import init_printing
init_printing()
from sympy import sqrt, sin, pi, var, integrate, Symbol, S, Integral
var("L x mu")
i = Symbol("i", integer=True)
j = Symbol("j", integer=True)
def phik(k, x):
return sqrt(2/L)*sin(k*pi*x/L) # domain is [0, L]
The basis functions are $\phi_i = \langle x|i\rangle$:
In [2]:
phik(i, x)
Out[2]:
In [3]:
-S(1)/(2*mu) * Integral(phik(i, x) * phik(i, x).diff(x, 2), (x, 0, L))
Out[3]:
In [4]:
_.doit()
Out[4]:
The off diagonal terms are zero:
In [5]:
-S(1)/(2*mu) * Integral(phik(i, x) * phik(j, x).diff(x, 2), (x, 0, L))
Out[5]:
In [6]:
_.doit()
Out[6]:
In [7]:
integrate(phik(1, x)*phik(2, x), (x, 0, L))
Out[7]:
In [8]:
integrate(phik(1, x)*phik(3, x), (x, 0, L))
Out[8]:
Which is equal to the solution in [1]: $$ T_{ij} = \cases{{1\over 2\mu}\left[i\pi\over L\right]^2 &if $i=j$;\cr 0,&otherwise.\cr} $$
In [9]:
Integral(phik(i, x) * x * phik(i, x), (x, 0, L))
Out[9]:
In [10]:
_.doit()
Out[10]:
But notice that:
In [11]:
Integral(phik(i, x) * x * phik(i, x), (x, -L/2, L/2))
Out[11]:
In [12]:
_.doit()
Out[12]:
The off diagonal terms are:
In [13]:
Integral(phik(i, x) * x * phik(j, x), (x, 0, L))
Out[13]:
In [14]:
_.doit()
Out[14]:
In [15]:
_.simplify()
Out[15]:
Which is equal to the solution in [1]: $$ X_{ij} = \cases{-{8ijL\over \pi^2 (i+j)^2(i-j)^2} &if $i+j$ is odd;\cr 0,&otherwise.\cr} $$
Authors of the code:
Original repository: https://github.com/OndrejMarsalek/Fourier-DVR-1D
In [16]:
import numpy as np
def _eig_sorted(M):
"""Convenience wrapper over numpy.linalg.eig.
:param M: matrix
:type M: 2D numpy array
Returns eigenvalues and eigenvectors of M sorted by eigevalue magnitude.
"""
val, vec = np.linalg.eig(M)
idx_sort = val.argsort()
return val[idx_sort], vec[:, idx_sort]
class Domain_Fourier_DVR_1D(object):
"""
Solves the Schroedinger equation on a finite one-dimensional interval using
the discrete variable representation method with the Fourier sine basis:
phi_k = sqrt(2 / (b-a)) * sin(k * pi * (x-a) / (b-a)), k=1..n_DVR
For details, see for example Appendix 2 of:
J. Chem. Phys. 104, 8433 (1996)
http://dx.doi.org/10.1063/1.471593
"""
def __init__(self, a, b, n_DVR):
"""Constructs the domain with given end points and basis size.
:param a: lower bound of domain
:param b: upper bound of domain
:param n_DVR: number of basis functions
"""
# store domain parameters
self._a = a
self._b = b
self._n_DVR = n_DVR
# no previous calculation was performed
self._m = None
self._V = None
# update spectral decomposition of the position operator
self._update_X_in_Fourier_basis()
def _update_X_in_Fourier_basis(self):
"""Decompose the position operator in the Fourier sine basis.
Eigenvalues and eigenvectors of the position operator
in the Fourier sine basis are stored in `self`.
"""
a = self._a
b = self._b
n_DVR = self._n_DVR
# construct position operator in Fourier sine basis
i = np.arange(1, n_DVR+1, dtype=float)[:, np.newaxis].repeat(n_DVR, axis=1)
j = i.transpose()
ipj = i + j
ipjmod1 = (ipj) % 2
div = np.where(ipjmod1, ((i-j) * ipj)**2, 1)
fact = -8.0 * (b-a) / np.pi**2
X_four = fact * np.where(ipjmod1, (i*j) / div, 0.0)
x, phi_four = _eig_sorted(X_four)
x = x + (a+b) / 2.0
self._x = x
self._phi_four = phi_four
def _update_T_four(self, m):
"""Build the kinetic energy operator and store it in `self`.
:param m: mass
"""
if self._m is None or (m != self._m):
self._m = m
l = self._b - self._a
t = (0.5 / m) * (np.pi / l)**2 * np.arange(1, self._n_DVR + 1)**2
self._T_four = np.diagflat(t)
def _update_V_four(self, V):
"""Build the potential energy operator and store it in `self`.
:param V: potential energy - real-space vectorized function
"""
if self._V is None or (V != self._V):
self._V = V
phi_four = self._phi_four
V_x = np.diagflat(V(self._x))
self._V_four = np.dot(np.dot(phi_four, V_x), phi_four.transpose())
def solve(self, m, V):
"""Solve the Schroedinger equation on this domain.
:param m: mass
:param V: potential energy - real-space vectoried function
Returns eigenenergies and eigenstates of the Hamiltonian sorted
by eigenenergy magnitude.
"""
# update kinetic energy and potential operators, if needed
self._update_T_four(m)
self._update_V_four(V)
# construct the Hamiltonian
H_four = self._T_four + self._V_four
# solve
E, psi_four = _eig_sorted(H_four)
return E, psi_four
def grid(self, x, psi_four):
"""Evaluate states on a real-space grid.
:param x: real-space grid - 1D array
:param psi_four: states in Fourier basis
Returns states on grid x.
"""
n_DVR, n_out = psi_four.shape
if n_DVR != self._n_DVR:
data = (self._n_DVR, n_DVR)
err = 'Wrong dimension of states. Expected %d, got %d.' % data
raise ValueError(err)
# convenience
a = self._a
b = self._b
# flag points outside the [a, b] interval
outside = np.logical_and((x > a), (x < b))
# construct Fourier sine basis functions on real-space grid
norm = np.sqrt(2 / (b-a))
four_1_grid = np.exp(1.0j * np.pi * (x-a) / (b-a))
four_1_grid *= outside
k = np.arange(1, n_DVR + 1, dtype=complex)
four_grid = norm * np.imag(four_1_grid[:, np.newaxis]**k)
# project states to real-space grid
psi_grid = np.dot(four_grid, psi_four).transpose()
return psi_grid
In [17]:
# settings
m = 1.0
x_min = -5.0
x_max = 5.0
n_DVR = 200 # increase this for convergence
n_g = 601
D0 = 43.0
a = 1.5
V = lambda x: (D0 / a**4) * (x**2 - a**2)**2
n_plot = 15
scale = 4.0
# solve
domain = Domain_Fourier_DVR_1D(x_min, x_max, n_DVR)
E, E_four = domain.solve(m, V)
# evaluate eigenstates on grid
x = np.linspace(x_min, x_max, n_g)
psi_x = domain.grid(x, E_four[:,:n_plot])
# print energies
for i, e in enumerate(E[:n_plot]):
print '%3i %18.12f' % (i, e)
# plot eigenstates
for i in range(n_plot):
plot([x[0], x[-1]], [E[i], E[i]], '--', color='gray')
plot(x, V(x), 'k-', lw=2)
for i in range(n_plot):
plot(x, scale * psi_x[i] + E[i])
xlim(-3, 3)
ylim(-E[0], E[n_plot])
xlabel("$x$ [a.u.]")
ylabel("E [a.u.]")
tight_layout()
In [18]:
# settings
m = 1.0
x_min = -15.0
x_max = 15.0
n_DVR = 300 # increase this for convergence
n_g = 1001
k = 1.0
V = lambda x: 0.5 * k * x * x
n_plot = 50
scale = 1.0
# solve
domain = Domain_Fourier_DVR_1D(x_min, x_max, n_DVR)
E, E_four = domain.solve(m, V)
# evaluate eigenstates on grid
x = np.linspace(x_min, x_max, n_g)
psi_x = domain.grid(x, E_four[:,:n_plot])
# print energies
for i, e in enumerate(E[:n_plot]):
print '%3i %12.6f' % (i, e)
# plot eigenstates
subplots_adjust(left=0.05, right=0.95,
bottom=0.05, top=0.95)
plot(x, V(x), 'k-', lw=2)
for i in range(n_plot):
plot([x[0], x[-1]], [E[i], E[i]], '--', color='gray')
for i in range(n_plot):
plot(x, scale * psi_x[i] + E[i])
xlim(x_min, x_max)
ylim(-E[0], E[n_plot])
xlabel("$x$ [a.u.]")
ylabel("E [a.u.]");
In [18]: