This page was generated from the Jupyter notebook hydrogen_wavefunction.ipynb.

Hydrogen radial wavefunction via numerov and from sympy

from typing import Any

import matplotlib.pyplot as plt
import numpy as np

from rydstate.radial import RadialKet
from rydstate.species.hydrogen import PotentialCoulombHydrogenTextBook

potential = PotentialCoulombHydrogenTextBook(l_r=5)
# analog one can also get this potential via: rydstate.species.get_potential_class("H_textbook")(l_r=5)

radial = RadialKet(
    nu=10, potential=potential, n_expected=10, dz=1e-2, integration_method="numerov", sign_convention="n_l_1"
)

from sympy.abc import r as sympy_r
from sympy.physics import hydrogen as sympy_hydrogen
from sympy.utilities.lambdify import lambdify


def calc_r_list(nu: float, l: int, x_list: np.typing.NDArray[Any]) -> np.typing.NDArray[Any]:
    if nu <= 35 or False:
        r_nl = lambdify(sympy_r, sympy_hydrogen.R_nl(nu, l, sympy_r, Z=1))
        r_list = r_nl(x_list)
    else:  # some weird sympy bug if trying to use lambdify R_nl for nu > 35
        r_list = np.zeros_like(x_list)
        for i, x in enumerate(x_list):
            r_list[i] = sympy_hydrogen.R_nl(nu, l, x, Z=1)
    return r_list


sympy_r_list = calc_r_list(radial.nu, radial.potential.l_r, radial.x_list)
sympy_u_list = sympy_r_list * radial.x_list
sympy_w_list = sympy_u_list / np.sqrt(radial.z_list)

fig, axs = plt.subplots(1, 3, figsize=(12, 5))

axs[0].plot(radial.z_list, radial.w_list, "C0-", label="Numerov")
axs[0].plot(radial.z_list, sympy_w_list, "C1--", lw=0.75, label="Sympy")
axs[0].set_xlabel(r"$z$")
axs[0].set_ylabel(r"$w(z)$")
axs[0].legend()

axs[1].plot(radial.x_list, radial.u_list, "C0-", label="Numerov")
axs[1].plot(radial.x_list, sympy_u_list, "C1--", lw=0.75, label="Sympy")
axs[1].set_xlabel(r"$r / a_0$")
axs[1].set_ylabel(r"$u(r) = r R(r)$")
axs[1].legend()

axs[2].plot(radial.x_list, radial.r_list, "C0-", label="Numerov")
axs[2].plot(radial.x_list, sympy_r_list, "C1--", lw=0.75, label="Sympy")
axs[2].set_xlabel(r"$r / a_0$")
axs[2].set_ylabel(r"$R(r)$")
axs[2].legend()

fig.tight_layout()
plt.show()