-2
$\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
    $\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$ – Mauro Vanzetto Jul 9 '18 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 Jul 9 '18 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$ – Wolfgang Bangerth Jul 9 '18 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$ – Ben Crowell Jul 9 '18 at 22:27
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 acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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