Solving Schrodinger equation numerically

Question

Here is the Schrodinger equation that is to be solved:

A 1D hard wall potential in $[0, 1]$. The potential within the potential well is given by a linear combination of Gaussian dips $$v(x) = - \sum_{i=1}^3a_ie^{-\frac{(x-b_i)^{2}}{2c_i^2}}$$

Thus the time-independent Schrodinger equation to be solved is $$(-\frac{h^2}{2m}\partial_x^{2}+v(x))\psi(x)=E\psi(x)$$

In the paper that I read, the interval $[0, 1]$ is discretized into dense grid points in the first place. How can I solve it numerically in such a way?

I welcome to scicomp.stackexchange. To encourage people to answer your question, I suggest to edit the posts with: What did you already try? What problems did you meet? Cite the paper. — Mauro Vanzetto, Jul 09 '18 at 08:58
This question has been asked previously in this site, have you checked the answers there. Do you have values for $b_i$ and $c_i$? — nicoguaro, Jul 09 '18 at 13:29
Also, there are numerous open source software projects that either solve this equation outright, or that have tutorial programs that do so. Have you looked around a bit? — Wolfgang Bangerth, Jul 09 '18 at 15:57
Maybe there is a reason for doing this on a grid, but if I were doing this, I would have thought of taking a basis of sine waves and diagonalizing the Hamiltonian matrix. I think the matrix elements would probably be expressible in terms of Erf. — , Jul 09 '18 at 22:27

nicoguaro · Accepted Answer · 2018-08-14T20:12:36.960

I saw @BenCrowell comment and wanted to check how hard it was.

If you write the variational formulation of the Schrödinger equation you would get

$$\Pi[\psi] = \frac{1}{2} \langle \nabla \psi, \nabla\psi\rangle + \langle \psi, V(x) \psi\rangle - E\langle \psi, \psi\rangle \, ,$$

and, choosing the finite expansion

$$\psi = \sum_{n=0}^N c_n \sin(n\pi x)\, ,$$

we obtain the following eigenvalue problem

$$[H]\lbrace c\rbrace = \frac{E^N}{2}\lbrace c\rbrace\, .$$

where

$$H_{mn} = \frac{\pi^2 m^2}{4}\delta_{mn} + U_{mn}\, ,$$

with

$$U_{mn} = \frac{\sqrt{2 \pi}}{8} \sum_{i=1}^3 a_i c_i f_i(m, n, b_i, c_i)\, ,$$

and

$$f_i(m, n, b, c) = {\rm erf} \left(1/2\,{\frac { \left( i \left( m-n \right) \pi\,{c}^{2} -b+1 \right) \sqrt {2}}{c}}\right){{\rm e}^{- \left( 1/2\, \left( m-n \right) \pi\,{c}^{2}+ib \right) \left( m-n \right) \pi}}-{\rm erf} \left(1/2\,{\frac { \left( i \left( m-n \right) \pi\,{c}^{2}-b \right) \sqrt {2}}{c}}\right){{\rm e}^{- \left( 1/2\, \left( m-n \right) \pi\,{c}^{2}+ib \right) \left( m-n \right) \pi}}+ \left( - {\rm erf} \left(1/2\,{\frac { \left( i \left( m-n \right) \pi\,{c}^{2} +b-1 \right) \sqrt {2}}{c}}\right)+{\rm erf} \left(1/2\,{\frac {\sqrt {2} \left( i \left( m-n \right) \pi\,{c}^{2}+b \right) }{c}}\right) \right) {{\rm e}^{ \left( -1/2\, \left( m-n \right) \pi\,{c}^{2}+ib \right) \left( m-n \right) \pi}}-{\rm erf} \left(1/2\,{\frac {\sqrt {2} \left( i \left( m+n \right) \pi\,{c}^{2}-b+1 \right) }{c}}\right){ {\rm e}^{- \left( 1/2\, \left( m+n \right) \pi\,{c}^{2}+ib \right) \left( m+n \right) \pi}}+{\rm erf} \left(1/2\,{\frac { \left( i \left( m+n \right) \pi\,{c}^{2}-b \right) \sqrt {2}}{c}}\right){ {\rm e}^{- \left( 1/2\, \left( m+n \right) \pi\,{c}^{2}+ib \right) \left( m+n \right) \pi}}-{{\rm e}^{ \left( m+n \right) \left( -1/2\, \left( m+n \right) \pi\,{c}^{2}+ib \right) \pi}} \left( {\rm erf} \left(1/2\,{\frac { \left( i \left( m+n \right) \pi\,{c}^{2}+b \right) \sqrt {2}}{c}}\right)-{\rm erf} \left(1/2\,{\frac { \left( i \left( m+n \right) \pi\,{c}^{2}+b-1 \right) \sqrt {2}}{c}}\right) \right)\, .$$

The following snippet solve the equation using this approach

from __future__ import division
import numpy as np
from scipy.linalg import eigh
from scipy.special import erf
from numpy import sqrt, pi, exp


def Uterm(m, n, a, b, c):
    aux = (erf((1j * (m - n) * pi * c ** 2 - b + 1) * sqrt(2) / c / 2) *
           exp(-(m - n) * pi * ((m - n) * pi * c ** 2 / 2 + 1j * b)) -
           erf((1j * (m - n) * pi * c ** 2 - b) * sqrt(2) / c / 2) *
           exp(-(m - n) * pi * ((m - n) * pi * c ** 2 / 2 + 1j * b)) +
           (-erf((1j * (m - n) * pi * c ** 2 + b - 1) * sqrt(2) / c / 2) +
            erf((1j * (m - n) * pi * c ** 2 + b) * sqrt(2) / c / 2)) *
            exp((m - n) * pi * (-(m - n) * pi * c ** 2 / 2 + 1j * b)) -
            erf((1j * (m + n) * pi * c ** 2 - b + 1) * sqrt(2) / c / 2) *
            exp(-(m + n) * ((m + n) * pi * c ** 2 / 2 + 1j * b) * pi) +
            erf((1j * (m + n) * pi * c ** 2 - b) * sqrt(2) / c / 2) *
            exp(-(m + n) * ((m + n) * pi * c ** 2 / 2 + 1j * b) * pi) -
            exp((-(m + n) * pi * c ** 2 / 2 + 1j * b) * (m + n) * pi) *
            (erf((1j * (m + n) * pi * c ** 2 + b) * sqrt(2) / c / 2) -
             erf((1j * (m + n) * pi * c ** 2 + b - 1) * sqrt(2) / c / 2)))
    return np.real(sqrt(2*pi)*a*c/4 * aux)


N = 11
a = [1, 2, 3]
b = [0.3, 0.6, 0.8]
c = 0.001
T = 0.5*pi**2*np.diag(range(1, N + 1))**2
U = np.zeros_like(T)
for m in range(1, N + 1):
    for n in range(1, N + 1):
        U[m - 1, n - 1] = Uterm(m, n, a[0], b[0], c) +\
                             Uterm(m, n, a[1], b[1], c) +\
                             Uterm(m, n, a[2], b[2], c)
H = T + U

vals, vecs = eigh(H)
print(np.round(np.real(vals[:11]), 6))

For $a=[1, 2, 3]$, $b=[0.3, 0.6, 08]$, $c=0.001$, the following are the results

[  4.952339  19.760811  44.43076   78.972835 123.375073 177.668876  241.822857 315.848936 399.736522 493.48025  597.12861 ]

That are really close to the ones computed using finite differences as proposed here,

[  4.932655  19.681995  44.253232  78.656791 122.880544 176.955286   240.848923 314.57281  398.115649 491.471978 594.688728]

Solving Schrodinger equation numerically

1 Answers1