How to calculate atomic potentials in SigSpace

Up until now, we have used callback functions to define potentials, such as rectangular or box-shaped potentials. These are useful for demonstrating wave function properties and building conceptual understanding. However, for more realistic problems, most potentials originate from point-like charges. Atomic potentials, for instance, follow a $ 1/r $ dependence due to Coulomb's law.

From a numerical perspective, such potentials pose challenges because they approach infinity as $ r = 0 $. While it is possible to avoid placing a grid point directly at $ r = 0 $, the extreme steepness near this point makes it difficult to achieve stable and accurate calculations.

Fortunately, analytical solutions exist for $ -Z/r $ potentials (often referred to as Hydrogen-like atoms). For these potentials, the ground state energy is given by $ E_0 = - 0.5 \cdot Z^2 $ where Z is the atomic number. For hydrogen ($Z=1$), the ground state energy simplifies to $ E_0 = -0.5 $. This result is consistent for both 2D and 3D cases.

Let’s attempt to reproduce this result using SigSpace by defining a potential function:

In [1]:
import numpy as np
import sigspace as sg

def Vfun(pts):
    return -1/np.sqrt(np.sum(np.square(pts), axis=1))

N_list = [20, 40, 60, 80]
E_list = []

for N in N_list:
    H = sg.Hamiltonian(sdims=2) + Vfun
    
    solver = sg.Solver(H)
    solver.grid((-15, 15, N))
    solver.print_info = True
    sol = solver.run()
    
    E_list.append(sol.energy.value)


sg.quickplot((N_list, E_list), layout={"xlabel": "Grid points [#]", "ylabel": "Energy [Ha]"})
Iter: 94, energy -0.583, residual 0.009878, duration 0.001340s
Iter: 375, energy -0.875, residual 0.009816, duration 0.077994s
Iter: 585, energy -1.068, residual 0.009903, duration 0.114318s
Iter: 780, energy -1.201, residual 0.009941, duration 0.135880s

The initial results are not very accurate. While the calculation starts near the expected value, increasing the grid resolution causes the energy to decrease further, effectively worsening the results. This occurs because a higher grid resolution amplifies the impact of the sharp peak at the center of the potential, resulting in steep gradients that are difficult for the solver to handle.

To address this issue, SigSpace provides atom objects that represent the Coulombic term directly. These objects are specifically designed to handle such potentials more effectively and can be added directly to the Hamiltonian. We can update the code as follows:

In [2]:
import sigspace as sg

N_list = [20, 40, 60, 80]
E_list = []

for N in N_list:
    H = sg.Hamiltonian(sdims=2) + sg.Atom(1, (0, 0))
    
    solver = sg.Solver(H)
    solver.grid((-15, 15, N))
    solver.print_info = True
    sol = solver.run()
    
    E_list.append(sol.energy.value)


sg.quickplot((N_list, E_list), layout={"xlabel": "Grid points [#]", "ylabel": "Energy [Ha]"})
Iter: 167, energy -0.532, residual 0.009804, duration 0.002076s
Iter: 478, energy -0.503, residual 0.009939, duration 0.021301s
Iter: 783, energy -0.501, residual 0.009973, duration 0.076016s
Iter: 1084, energy -0.501, residual 0.009997, duration 0.189022s

Instead of using a callback function, sg.Atom(1, (0, 0)) can be used to insert a hydrogen atom at the center (0, 0). The solver employs an internal representation to efficiently manage the infinite peak at the center of the Coulomb potential, ensuring stable and accurate calculations.

The resulting graph demonstrates the effectiveness of this approach and highlights the limitations of adding Coulombic potentials via callback functions. While it is still possible to combine Atoms with callback functions, using sg.Atom for Coulomb potentials is strongly recommended for better performance and accuracy.

In [3]:
sg.quickplot(sol)

A Note on 1D Potentials

In principle, sg.Atom can be used for problems in all dimensions. However, one-dimensional (1D) cases are unique, and attempting to directly use the code from higher dimensions in 1D would result in an error.

In 2D and higher dimensions, the wave function has a finite energy for a $ 1/r $ potential. In 1D, however, the energy diverges for such potentials, making it necessary to impose a limit on the potential to achieve stable results.

A limit can be specified in the Hamiltonian using H.coulomb_limit = 10. This setting caps the potential to a maximum value of 10, as described by the following formula:

\begin{equation} V = \frac{-Z}{\sqrt{r^2 + \frac{1}{C^2}}} \end{equation}

For instance, at $ r = 0 $ the potential is $ V = -Z \cdot C $.

In [4]:
import sigspace as sg

H = sg.Hamiltonian(sdims=1) + sg.Atom(1, 0)
H.coulomb_limit = 10

solver = sg.Solver(H)
solver.grid(-10, 10, 201)
sol = solver.run()

sg.quickplot(sol)

And to plot the corresponding potential:

In [5]:
import numpy as np

xw = np.linspace(-10, 10, 1000)

sg.quickplot((xw, H.potential(xw)), layout={"xlabel": "x", "ylabel": "Potential"})