Calculating energy of 1D Schrödinger equation using Sympy
I was reading Szabo’s book “Modern Quantum Chemistry” and saw the exercise questions that seem to be solvable via programming. Hence I decided to give it a try and pick up sympy at the same time. I was a mathematica user back in undergrad, but have not used it ever since. Maybe sympy will just do the same trick.
Here is the question (exercise 1.18 from the book)
The Schrödinger equation (in atomic units) of an electron moving in one dimension under the influence of the potential $$-\delta(x)$$ is
$$ \left(-\frac{1}{2}\frac{d^2}{dx^2} - \delta(x)\right) | \Phi\rangle = \epsilon | \Phi\rangle $$
Use the variation method with the trial function
$$ | \tilde \Phi \rangle = N e^{-ax^2} $$
to show that $$-\pi^{-1}$$ is an upper bound to the exact group state energy (which is -0.5).
From the variational principle, we know if the normalized wavefunction satisfies the appropriate boundary condition, then the expectation of the Hamiltonian is an upper bound to the exact ground state energy. In math expressions, if
$$ \langle\tilde\Phi|\tilde\Phi \rangle = 1 $$
then
$$ \langle\tilde\Phi| \mathcal{H} | \tilde\Phi \rangle \ge \epsilon_0 $$
I will show how we can use sympy to solve this problem
from sympy import * # I dislike this way to import everything, but it seems to be common in sympy
init_printing()
x = symbols("x") # define the symbols
N, a = symbols("N a", positive=True) # a should be positive to satisfy boundary cond at infinity
phi = N * exp(-a * x ** 2) # this is our trial function
1. Normalization conditions
We will need to normalize the wavefunction
$$ \langle\tilde\Phi|\tilde\Phi \rangle = 1 $$
phi2_int = integrate(phi * phi, (x, -oo, oo)) # the integration of phi * phi. Our function is real here
phi2_int
$$\frac{\sqrt{2} \sqrt{\pi} N^{2}}{2 \sqrt{a}}$$
This expression equals to 1 from our normalization condition. Hence we can solve for $$N$$
n_cond = solve(phi2_int-1)
n_cond
$$\left [ \left { N : \frac{\sqrt[4]{2} \sqrt[4]{a}}{\sqrt[4]{\pi}}\right }\right ]$$
2. Calculate the Hamiltonian with the trial function
term1 = integrate(phi * (-0.5 * diff(diff(phi, x), x)), (x, -oo, oo))
term2 = integrate(-phi * DiracDelta(x) * phi, (x, -oo, oo))
H = term1 + term2
H
$$0.25 \sqrt{2} \sqrt{\pi} N^{2} \sqrt{a} - N^{2}$$
3. Substitute the normalization condition
H_sol = H.subs(n_cond[0])
H_sol
$$- \frac{\sqrt{2} \sqrt{a}}{\sqrt{\pi}} + 0.5 a$$
This expression still contains $$a$$. To find the minimum of this equation, we will need to solve for $$a$$
4. Minimize with respect to $$a$$
$$a$$ is minimal when $$\partial H_{sol}/\partial a = 0$$
a_sol = solve(diff(H_sol, a))[0]
a_sol
$$0.636619772367581$$
substitute $$a$$ solution into the solution for $$H$$, we get
H_sol.subs({a: a_sol}).evalf()
$$-0.318309886183791$$
which is exactly $$-\pi^{-1}$$
In summary, in this notebook, I show how we can use sympy to solve simple Schrödinger equation. Sometimes, using sympy can be unintuitive especially if the bounds of the variables are not properly set. In that case, you will get piecewise function results, and you will need to manually select the correct solutions.
I found the use of expression oo to represent infinity quite interesting and brilliant.