orosmatthew/quantum-dev
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add numerov solve

Browse Source
This commit is contained in:
orosmatthew 2024-04-26 13:46:38 -04:00
parent 4f5ab5e8ee
commit ca675b83ca
1 changed files with 155 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

155
numerov_solve.py Normal file
View File

@ -0,0 +1,155 @@
"""
MO 26 April 2024
modified from: plh 26 April 2024
Template for Numerical integration using Numerov method.
Ported from SHOsolve.py
See https://physics.bu.edu/~py502/lectures4/schrod.pdf and
https://mtdevans.com/2013/07/solving-the-schrodinger-equation-with-numerovs-algorithm/
"""
# 1-D Schrodinger wave equation
# (-hbar^2/2m) d^2/dx^2 psi(x) + V(x) psi(x) = E psi(x)
#
# Spherically symmetric Psi(xvec)_L,Lz,n = R(r)_L,n Y(phi,theta)_L,Lz
# [(-hbar^2/2mu)(1/r)d^2/dr^2 + L(L+1)hbar^2/2mu-r^2 + V(r)] R(r)_L,n = E_L,n R(r)_L,n
# Define U(r) = r R(r)
# [(-hbar^2/2mu) d^2/dr^2 + L(L+1)hbar^2/2mu-r^2 + V(r)] U(r)_L,n = E_L,n U(r)_L,n
#
from typing import Callable
import numpy as np
import pylab as pl
# version 0.3
# version history
# 0.1 - ported from BohrSolve.py of 2018(ca), retains constants
# 0.2 - simplified to calculating psi(0) for a list of energies
# 0.3 - version used for 2024SP PS
# constants
hbar = 6.626070040e-34 / np.pi / 2 # J-s = Planck's constant
light_speed = 299792458 # m s-1 = speed of light
bohr_radius = 0.52917721067e-10 # m = Bohr radius
mass_proton = 1.672621898e-27 # kg = mass of proton
mass_electron = 9.10938356e-31 # kg = mass of electron
mass_electron_reduced = mass_electron / (1.0 + mass_electron / mass_proton) # kg = reduced mass of electron in H atom
electrical_constant = 9.0e9 # N-m^2/s^2 = electrical constant
charge = 1.602176634e-19 # C = electronic charge
potential_1 = -electrical_constant * charge * charge # = potential at r=1m
potential_0 = potential_1 / bohr_radius # = potential at r=a0
energy_1_ev = 13.6 # eV = ground state energy of H atom in eV
energy_1_j = energy_1_ev * charge # J = ground state energy of H atom in J = 21.88E-19 J
# initial values
intervals = 10000 # number of intervals
radius_max = 200.0 * bohr_radius # maximum radius = starting point
radius_min = bohr_radius / 500000.0 # minimum radius = ending point
# dr = a0/1000000.0 # radius step
energy = -energy_1_j # energy, turning point is where E=V(x)
energy_min = -1 * energy_1_j # min energy for scan
energy_max = -0.005 * energy_1_j # max energy for scan
energy_intervals = 1000
def bohr_potential(r):
"""Models the Bohr atom electrostatic potential"""
return potential_1 / r
def do_numerov(radial_vals: np.ndarray, energy: float, initial: tuple[float, float],
potential_func: Callable[[float], float], angular_momentum=0.0) -> np.ndarray[float]:
"""
Solves the differential equation F''(r) + g(r) F(r) = 0 using the Numerov method.
In the context of the Schrödinger equation, g(r) = k^2 = (E - V(r))*2/m
Here, the argument 'potential' points to a function definition that calculates the potential V(r).
"""
hbar_sqrd_over_2m = hbar * hbar / 2 / mass_electron_reduced
step_size = radial_vals[1] - radial_vals[0]
solution = np.zeros(len(radial_vals))
solution[0], solution[1] = initial
for i in range(2, len(radial_vals)):
g_factor = step_size * step_size * (((energy - potential_func(float(radial_vals[i]))) / hbar_sqrd_over_2m)
- (angular_momentum * (angular_momentum + 1) / radial_vals[i] ** 2)) / 12
solution[i] = ((solution[i - 1] * (2.0 - 10.0 * g_factor) - solution[i - 2] * (1.0 + g_factor))
/ (1.0 + g_factor))
# Normalize the solution so that the total probability equals 1 in the center
normalization_factor = solution.max()
solution = solution / normalization_factor
return solution
def opposite_signs(a, b):
return a * b < 0
def main():
# Main program begins here
x_vals = np.linspace(radius_max, radius_min, intervals) # generate x values
dx = abs(x_vals[1] - x_vals[0])
energies = np.linspace(energy_min, energy_max, energy_intervals).tolist()
results = [] # these will be wave function(at r near 0) for each of the energies
zeros = []
print(f'\nTesting {len(energies)} energies from {energy_min} to {energy_max} for solutions...')
for i, energy in enumerate(energies):
if i % 100 == 0:
print(f' - Testing energy: {i}/{len(energies)} = {energy}...')
initial = (0.0, 1.00) # initial value of function and slope, at large r
numerov = do_numerov(x_vals, energy, initial, bohr_potential) # does integration, returns wave function
if len(results) != 0:
if opposite_signs(results[-1], numerov[-1]):
zero = (energy + energies[i - 1]) / 2
print(f'\tFound solution at {zero}')
zeros.append(zero)
results.append(numerov[-1])
print(f'Done testing energies\n')
# guess at best solution E
# first plot the results as a function of energy
pl.figure(4)
pl.xlabel('Energy')
pl.ylabel('Wave function near r=0')
pl.title('Trial results vs Energy')
for zero in zeros:
pl.axvline(x=zero, color='red')
pl.plot(energies, results)
for n, energy_solution in enumerate(zeros):
print(f'Trying E = {energy_solution:.3e} J, {energy_solution / charge:.3} eV')
initial = (0.0, 1.0)
func_reduced = do_numerov(x_vals, energy_solution, initial, bohr_potential) # U(r) = reduced wave function
# Figure 1 = reduced radial wave function
pl.figure(1)
pl.xlabel('r')
pl.ylabel('U(r)')
pl.title('Reduced U(r) = r R(r)')
pl.plot(x_vals, func_reduced, label=f'n={n:d}')
pl.legend()
prob = dx * func_reduced.dot(func_reduced) # calculate total probability, should we x 4pi?
func_reduced = func_reduced / np.sqrt(prob) # total probability should = 1.0!
func_radial = func_reduced / x_vals # get the radial wave function R(r)
mean_x = dx * (x_vals * func_reduced).dot(func_reduced)
print(f'Expectation value <r> = {mean_x:.3e} m')
# Figure 2 = actual radial wave function R(r)
pl.figure(2)
pl.title('Radial wave function R(r) = U(r)/r')
pl.plot(x_vals, func_radial, label=f'n={n:d}')
pl.xlabel('r')
pl.ylabel('R(r)')
# pl.ylim(-1,1.1*E)
pl.legend()
# Figure 3 = probability integrand r^2 R(r)^2 = U^2
pl.figure(3)
pl.title('Probability Integrand [r^2 R(r)^2]')
pl.plot(x_vals, func_reduced * func_reduced, label=f'n={n:d}')
pl.xlabel('r')
pl.ylabel('P(r)')
# pl.ylim(-1,1.1*E)
pl.legend()
pl.show() # show graph window until closed
main()