-1
$\begingroup$

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?

$\endgroup$
4
  • 4
    $\begingroup$ 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. $\endgroup$ Commented Jul 9, 2018 at 8:58
  • $\begingroup$ 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$? $\endgroup$
    – nicoguaro
    Commented Jul 9, 2018 at 13:29
  • $\begingroup$ 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? $\endgroup$ Commented Jul 9, 2018 at 15:57
  • 1
    $\begingroup$ 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. $\endgroup$
    – user28109
    Commented Jul 9, 2018 at 22:27

1 Answer 1

3
$\begingroup$

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]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.