qmsolver 0.1.4

A quantum mechanics solver library

QMSolver

QMSolver is a Python library for numerically solving the time-independent Schrödinger equation in one dimension. The library implements the finite differences method to compute energy eigenvalues and eigenfunctions for quantum mechanical systems with arbitrary potential energy functions.

Table of Contents

• Installation
• Usage
• Method Limitations
• Custom Potential Implementation
• SI Units
• Results Validation
  • Infinite Square Well Limit
  • Harmonic Oscillator Validation
• Development and Contributing
• License
• References

Installation

PyPI Installation

pip install qmsolver

Requirements: Python >= 3.9

Repo Installation

1. Clone the repository:

   git clone https://github.com/gbatagian/qmsolver.git
   cd qmsolver
   

2. Install the source code:

   pip install -e .
   

Usage

Finite Differences Solver (FDSolver)

The qmsolver package provides the tise module, which contains the FDSolver class for solving the time-independent Schrödinger equation using finite differences. To use it, create an FDSolver instance and inject a potential class into its potential_generator attribute. Sample potentials are provided in the potentials module. For example, the finite square well solution is demonstrated below:

from qmsolver.tise import FDSolver
from qmsolver.potentials import FiniteSquareWellPotential

solver = FDSolver(steps=2_000, x_min=-5, x_max=5, n_lowest=7)
potential = FiniteSquareWellPotential(
    x_grid=solver.x_grid, well_depth=25, well_width=2
)
solver.potential_generator = potential
solver.solve()
solver.output()
solver.plot()

****************************************

-> 7 lowest energy states:

      🔒   E(0)  =  -24.0554383661  (bound)
      🔒   E(1)  =  -21.2408973698  (bound)
      🔒   E(2)  =  -16.6240464549  (bound)
      🔒   E(3)  =  -10.3701336891  (bound)
      🔒   E(4)  =   -2.9999825644  (bound)
      🌊   E(5)  =    0.2854054853  (free) 
      🌊   E(6)  =    0.3280462337  (free) 

****************************************

⚠️ Attention: The FiniteSquareWellPotential class requires a spatial grid as input. It is recommended to provide the x_grid of the solver to ensure the potential and solver use the same grid.

This will generate a plot showing the potential (black line) and the first few bound eigenstates:

Note: The plot() method only displays bound states (i.e. $E < V_{asymptotic}$, where $V_{asymptotic}$ is the min value at the boundaries), as these are the physically meaningful solutions under the finite differences method. Free states are affected by the zero boundary conditions which result grid-dependent "artificial quantization" of the continuum energy spectrum. For more information, see the Method Limitations section

Note: You can customize the physical constants by setting the solver.h_bar and solver.m attributes to your desired values before calling solver.solve(). By default, they are set to 1 (dimensionless units). For more information on using SI units, see the SI Units section.

• After calling solver.solve(), the solver will have the E_lowest attribute containing the n-lowest eigenenergies and the Psi_lowest attribute containing the corresponding eigenfunctions.

  solver.E_lowest
  

  array([-24.05543837, -21.24089737, -16.62404645, -10.37013369, -2.99998256, 0.28540549])
  

  solver.Psi_lowest 
  

  array([[ 7.72546314e-16, -7.70210475e-15,  2.04525254e-13,
          2.68267151e-11,  7.52877452e-08,  1.30754252e-04],
      [ 1.54602276e-15, -1.54123977e-14,  4.09220681e-13,
          5.36673540e-11,  1.50586795e-07,  2.61506636e-04],
      [ 2.32136057e-15, -2.31390756e-14,  6.14256593e-13,
          8.05358477e-11,  2.25908455e-07,  3.92255285e-04],
      ...,
      [ 2.32136057e-15,  2.31390756e-14,  6.14256593e-13,
          -8.05358477e-11,  2.25908455e-07, -3.92255285e-04],
      [ 1.54602276e-15,  1.54123977e-14,  4.09220681e-13,
          -5.36673540e-11,  1.50586795e-07, -2.61506636e-04],
      [ 7.72546314e-16,  7.70210475e-15,  2.04525254e-13,
          -2.68267151e-11,  7.52877452e-08, -1.30754252e-04]],
      shape=(2000, 6))
  

• Additionally, the solver provides E_bound and Psi_bound attributes containing only the bound states:

  solver.E_bound
  

  array([-24.05543837, -21.24089737, -16.62404645, -10.37013369, -2.99998256])
  

  solver.Psi_bound
  

  array([[ 7.72546314e-16, -7.70210475e-15,  2.04525254e-13,
           2.68267151e-11,  7.52877452e-08],
      [ 1.54602276e-15, -1.54123977e-14,  4.09220681e-13,
           5.36673540e-11,  1.50586795e-07],
      [ 2.32136057e-15, -2.31390756e-14,  6.14256593e-13,
           8.05358477e-11,  2.25908455e-07],
      ...,
      [ 2.32136057e-15,  2.31390756e-14,  6.14256593e-13,
          -8.05358477e-11,  2.25908455e-07],
      [ 1.54602276e-15,  1.54123977e-14,  4.09220681e-13,
          -5.36673540e-11,  1.50586795e-07],
      [ 7.72546314e-16,  7.70210475e-15,  2.04525254e-13,
          -2.68267151e-11,  7.52877452e-08]], shape=(2000, 5))
  

Method Limitations

The finite difference method is well-suited for computing bound states of quantum systems, but has some limitations for free states (scattering or continuum states).

The numerical implementation imposes zero boundary conditions at the edges of the spatial grid (x_min and x_max). This effectively encloses the system within an infinite square well of width (x_max - x_min), which introduces artificial quantization of the continuum energy spectrum.

For bound states, the wavefunctions decay exponentially to zero outside the well region. Since they naturally satisfy the zero boundary conditions at the grid edges, the computed energies and wave functions are unaffected by the grid size, provided the grid extends sufficiently far to capture the exponential decaying tail.

For free states, the wavefunctions oscillate and do not decay to zero. The artificial boundary conditions at the grid edges cause reflection of the wavefunction, creating standing waves that depend on the grid length. This leads to:

• Quantization of the continuum spectrum
• Grid-dependent energy levels and wavefunctions

FDSolver is better-suited for bound state problems and provides reliable results within that domain - as long as the grid extends sufficiently beyond the potential region to minimize boundary effects. For free states, the numerical solutions behave as if the system is confined within an infinite potential well, with complete wavefunction reflection at the grid boundaries, leading to quantization of the continuum energy spectrum.

Custom Potential Implementation

While the potentials module provides several predefined potential classes, you can also implement custom potentials by inheriting from the BasePotential abstract base class. This allows you to solve the Schrödinger equation for arbitrary potential energy functions. To create a custom potential class:

1. Inherit from BasePotential: Your class must inherit from potentials.base.BasePotential
2. Implement the generate() method: This method should return a NumPy array containing the potential energy values evaluated on the spatial grid
3. Accept grid as input: The __init__ method should accept the spatial grid (x_grid) as a parameter
4. Instance attribute parameters: Any potential parameters (depths, widths, etc.) should be stored as instance attributes

Below follows an example implementation of a sinusoidal potential well:

⚠️ Important: All parameters required by the generate() method must be provided through the __init__ method and stored as instance attributes (accessed through the self namespace). The generate() method should not accept additional parameters - it should only return the potential array using the stored parameters and the grid.

import numpy as np

from qmsolver.potentials import BasePotential
from qmsolver.tise.finite_differences import FDSolver


class SinusoidalWellPotential(BasePotential):
    """
    A composite sinusoidal potential with well and barrier regions:

        V(x) = A/2 + V_well(x) + V_barrier(x)

    where:
        - V_well(x)    = -A * |sin(x)| for      |x| ≤ π          , else 0
        - V_barrier(x) =  A * |sin(x)| for  π < |x| ≤ (5/6 + 2)π , else 0
        - A/2                          for      |x| > (5/6 + 2)π , else 0
    """

    def __init__(self, x_grid: np.array, amplitude: float) -> None:
        """
        Parameters:
        - x_grid: Spatial grid points
        - amplitude: Amplitude of the sinusoidal modulation (A > 0)
        """
        self.x_grid = x_grid
        self.amplitude = amplitude

    def generate(self) -> np.array:
        """
        Generate the potential energy array.

        Returns:
            np.array: Potential energy values on the grid
        """
        return (
            np.where(
                np.abs(self.x_grid) <= np.pi,
                -self.amplitude * np.abs(np.sin(self.x_grid)),
                0,
            )
            + np.where(
                (np.abs(self.x_grid) >= np.pi)
                & (np.abs(self.x_grid) <= (5 / 6 + 2) * np.pi),
                self.amplitude * np.abs(np.sin(self.x_grid)),
                0,
            )
            + np.where(
                np.abs(self.x_grid) >= (5 / 6 + 2) * np.pi,
                self.amplitude / 2,
                0,
            )
        )


solver = FDSolver(steps=2000, x_min=-4 * np.pi, x_max=4 * np.pi, n_lowest=15)
potential = SinusoidalWellPotential(x_grid=solver.x_grid, amplitude=5)
solver.potential_generator = potential
solver.solve()
solver.output()
solver.plot()

****************************************

-> 15 lowest energy states:

      🔒   E(0)  =   -3.9318637166  (bound)
      🔒   E(1)  =   -3.9000190100  (bound)
      🔒   E(2)  =   -1.9641284920  (bound)
      🔒   E(3)  =   -1.6914279178  (bound)
      🔒   E(4)  =   -0.1694428666  (bound)
      🔒   E(5)  =    0.5585681892  (bound)
      🔒   E(6)  =    1.8412337091  (bound)
      🔒   E(7)  =    2.2651012660  (bound)
      🔒   E(8)  =    2.2653605113  (bound)
      🌊   E(9)  =    2.7548549836  (free) 
      🌊   E(10) =    2.7548579298  (free) 
      🌊   E(11) =    2.7970433115  (free) 
      🌊   E(12) =    3.4976640375  (free) 
      🌊   E(13) =    3.4976741598  (free) 
      🌊   E(14) =    3.9044940884  (free) 

****************************************

SI Units

By default, FDSolver solves the Schrödinger equation in dimensionless units. However, it is also possible to perform the calculations in SI units. To use SI units:

1. Express the potential energy in Joules and the spatial grid in meters.
2. Set the h_bar and m attributes of the solver to their SI values (e.g. using scipy.constants).

Below follows an example for an electron in a finite square well:

import numpy as np
from scipy import constants
from qmsolver.tise import FDSolver
from qmsolver.potentials import FiniteSquareWellPotential

well_depth_ev = 1.0     # Well depth in electron volts
well_width_nm = 1.0     # Well width in nanometers

# Convert to SI units
well_depth_joules = well_depth_ev * constants.e     # Convert eV to Joules
well_width_meters = well_width_nm * 1e-9            # Convert nm to meters

# Spatial domain in meters
x_min_m = -3e-9
x_max_m = 3e-9

solver = FDSolver(steps=2_000, x_min=x_min_m, x_max=x_max_m, n_lowest=3)

# Set physical constants in SI units
solver.h_bar = constants.hbar  # Reduced Planck's constant in J⋅s
solver.m = constants.m_e       # Electron mass in kg

potential = FiniteSquareWellPotential(
    x_grid=solver.x_grid,
    well_depth=well_depth_joules,
    well_width=well_width_meters,
)
solver.potential_generator = potential

solver.solve()
solver.output()
solver.plot(is_dimensionless=False, scale=1e19, energy_units="J")

****************************************

-> 3 lowest energy states:

      🔒   E(0)  = -1.29761325178584494211e-19  (bound)
      🔒   E(1)  = -4.79509995005876440109e-20  (bound)
      🌊   E(2)  = 7.54877077796044010905e-21  (free) 

****************************************

Convert to eV

After solving in SI units, simply divide the eigenenergies by the electron charge to get values in eV:

E_lowest_ev = np.array(solver.E_lowest) / constants.e
print("\nEnergies in electron volts:")
for i, energy in enumerate(E_lowest_ev):
    print(f"E({i}) = {energy:.8f} eV")

Energies in electron volts:
E(0) = -0.80990649 eV
E(1) = -0.29928660 eV
E(2) = 0.04711572 eV

Results Validation

Infinite Square Well Limit

To validate the accuracy of FDSolver we can test it against a known analytical solution. A good test case is the infinite square well limit of a finite square well.

When a finite square well has a very large depth it tends to the infinite square well with analytical solutions:

$$E_n = n^2 \frac{\hbar^2 \pi^2}{2mL^2}$$

where:

• $n = 1, 2, 3, \dots$ (quantum number)
• $L$ is the well width
• $\hbar$ is the reduced Planck constant
• $m$ is the particle mass

Thus, we can solve a finite square well with very large depth (1e12 eV) using FDSolver and compare the results to the infinite square well formula. The finite well with very large depth eigenenergies should approach the infinite well eigenenergies.

import numpy as np
from scipy import constants

from qmsolver.potentials import FiniteSquareWellPotential
from qmsolver.tise import FDSolver

well_depth_ev = 1e12  # Well depth in electron volts
well_width_nm = 1.0  # Well width in nanometers

# Convert to SI units
well_depth_joules = well_depth_ev * constants.e  # Convert eV to Joules
well_width_meters = well_width_nm * 1e-9  # Convert nm to meters

# Spatial domain in meters
x_min_m = -3e-9
x_max_m = 3e-9

solver = FDSolver(steps=10_000, x_min=x_min_m, x_max=x_max_m, n_lowest=5)

# Set physical constants in SI units
solver.h_bar = constants.hbar  # Reduced Planck's constant in J⋅s
solver.m = constants.m_e  # Electron mass in kg

potential = FiniteSquareWellPotential(
    x_grid=solver.x_grid,
    well_depth=well_depth_joules,
    well_width=well_width_meters,
)
solver.potential_generator = potential

solver.solve()
solver.output()

E_lowest_ev = np.array(solver.E_lowest) / constants.e
print("\nEnergies in electron volts:")
for i, energy in enumerate(E_lowest_ev):
    infinite_square_Well_energy = (
        (i + 1) ** 2
        * constants.hbar**2
        * np.pi**2
        / (2 * solver.m * well_width_meters**2)
    ) / constants.e

    renormalized_energy = energy + well_depth_ev

    error = (
        100
        * abs(renormalized_energy - infinite_square_Well_energy)
        / infinite_square_Well_energy
    )

    print(
        f"E({i}) = {renormalized_energy:.8f} eV | E_ISW({i}): {infinite_square_Well_energy:.8f} eV | Error: {error:.8f} %"
    )

****************************************

-> 5 lowest energy states:

      🔒   E(0)  =   -0.0000001602  (bound)
      🔒   E(1)  =   -0.0000001602  (bound)
      🔒   E(2)  =   -0.0000001602  (bound)
      🔒   E(3)  =   -0.0000001602  (bound)
      🔒   E(4)  =   -0.0000001602  (bound)

****************************************

Energies in electron volts:
E(0) = 0.37585449 eV | E_ISW(0): 0.37603016 eV | Error: 0.04671698 %
E(1) = 1.50329590 eV | E_ISW(1): 1.50412065 eV | Error: 0.05483270 %
E(2) = 3.38232422 eV | E_ISW(2): 3.38427146 eV | Error: 0.05753794 %
E(3) = 6.01281738 eV | E_ISW(3): 6.01648259 eV | Error: 0.06091950 %
E(4) = 9.39514160 eV | E_ISW(4): 9.40075405 eV | Error: 0.05970214 %

Based on the comparison of the FDSolver results with the equivalent analytical solutions, we can conclude that the developed numerical method provides reliable results for bound state problems when using appropriate grid resolution. In the examined case above, the numerical results match the expected analytical values with an error less than 0.1%. It's worth noting that:

• In principle, accuracy improves with higher grid resolution
• However, accuracy depends on the specific grid spacing rather than just the total number of points. For instance, a grid with 6,000 steps might yield lower accuracy than one with 4,000 steps and a grid with 7,000 steps better from both. This behavior is problem-dependent and can vary between different potential functions.

The table below shows the average error for the first five eigenenergies across different grid resolutions for the validation test case above:

Steps	Error (%)
1,000	0.62310486
2,000	1.08952827
3,000	0.45091440
4,000	0.14786237
5,000	0.44736513
6,000	0.24246108
7,000	0.08640829
8,000	0.27421701
9,000	0.17126452
10,000	0.055941852

Harmonic Oscillator Validation

To further validate the accuracy of FDSolver, we can test it against the analytical solution for the harmonic oscillator.

The harmonic oscillator has exact analytical solutions:

$$E_n = \hbar \omega \left(n + \frac{1}{2}\right)$$

where:

• $n = 0, 1, 2, \dots$ (quantum number)
• $\omega$ is the angular frequency
• $\hbar$ is the reduced Planck constant

We can solve the harmonic oscillator potential using FDSolver and compare the numerical results to the analytical formula.

import numpy as np
from scipy import constants

from qmsolver.potentials import HarmonicOscillatorPotential
from qmsolver.tise import FDSolver

x_min_m = -10e-9
x_max_m = 10e-9

solver = FDSolver(steps=10_000, x_min=x_min_m, x_max=x_max_m, n_lowest=5)

# Set physical constants in SI units
solver.h_bar = constants.hbar
solver.m = constants.m_e

# Create harmonic oscillator potential active over entire grid
omega = 1e14  # Angular frequency (rad/s)
m = constants.m_e
potential = HarmonicOscillatorPotential(
    x_grid=solver.x_grid,
    spring_constant=m * omega**2,
    grid_active_range=1,  # Active over entire grid
)
solver.potential_generator = potential

solver.solve()
solver.output()

E_lowest_ev = np.array(solver.E_lowest) / constants.e
print("\nEnergies in electron volts:")
s = 0
for i, energy in enumerate(E_lowest_ev):
    analytical_energy = (constants.hbar * omega * (i + 0.5)) / constants.e
    renormalized_energy = energy + np.abs(np.min(solver.potential)) / constants.e
    error = 100 * abs(renormalized_energy - analytical_energy) / analytical_energy
    s += error
    print(
        f"E({i}) = {renormalized_energy:.8f} eV | E_HO({i}): {analytical_energy:.8f} eV | Error: {error:.8f} %"
    )

print(f"\nAverage error: {s/len(E_lowest_ev):.8f} %")

****************************************

-> 5 lowest energy states:

      🔒   E(0)  = -4.50196327746029571420e-19  (bound)
      🔒   E(1)  = -4.39650614125021047605e-19  (bound)
      🔒   E(2)  = -4.29104905059774577557e-19  (bound)
      🔒   E(3)  = -4.18559200550213846300e-19  (bound)
      🔒   E(4)  = -4.08013500595898489953e-19  (bound)

****************************************

Energies in electron volts:
E(0) = 0.03291056 eV | E_HO(0): 0.03291060 eV | Error: 0.00010800 %
E(1) = 0.09873173 eV | E_HO(1): 0.09873179 eV | Error: 0.00006480 %
E(2) = 0.16455287 eV | E_HO(2): 0.16455299 eV | Error: 0.00007344 %
E(3) = 0.23037398 eV | E_HO(3): 0.23037418 eV | Error: 0.00008948 %
E(4) = 0.29619506 eV | E_HO(4): 0.29619538 eV | Error: 0.00010800 %

Based on the comparison of the FDSolver results with the analytical harmonic oscillator solutions, the numerical method demonstrates excellent accuracy for this potential as the numerical results match the expected analytical values with an error less than 0.0001% (6th decimal digit accuracy).

The accuracy in this case is much higher than the infinite square well limit validation that took place in the previous section, because the harmonic oscillator potential definition does not approximate an ideal physical system (unlike the very deep finite square well, which was a limit approximation). In the lowest energy range of the harmonic oscillator, where boundary effects are almost non-existent, the numerical description is exceptionally accurate.

The table below shows the average error for the first five eigenenergies across different grid resolutions for the harmonic oscillator test case above:

Steps	Error (%)
1,000	0.00889065
2,000	0.00222036
3,000	0.00098649
4,000	0.00055481
5,000	0.00035504
6,000	0.00024654
7,000	0.00018112
8,000	0.00013867
9,000	0.00010956
10,000	0.00008874

Development and Contributing

Based on the comparison of the FDSolver results with the analytical harmonic oscillator solutions, the numerical method demonstrates excellent accuracy for this potential as the numerical results match the expected analytical values with an error less than 0.001%.

Development and Contributing

Build and Test

1. Build the package: make build
2. Install in development mode: make install-dev
3. Run the test suite: make test

Virtual Environment Setup

1. Python >= 3.9 should be installed
2. Pipenv should be installed (if not: pip install pipenv)
3. Create the virtual environment: make venv
4. Run the test suite: make test

⚠️ Important: For development, it is recommended to build and test from within the virtual environment.

Additional Development Commands

• Format code: make reformat (runs black and isort)
• Run tests with coverage: make coverage
• Create the virtual environment: make venv

License

This project is licensed under the MIT License - see the LICENSE file for details.

Important Notes:

• This is a personal project developed in my free time
• No warranty is provided
• For serious scientific work, please verify results independently and consult domain experts

References

1. Computational Quantum Mechanics
   Joshua Izaac, Jingbo Wang
   Springer, Chapter 9.6: The Direct Matrix Method

2. Solving the time-dependent Schrödinger equation using finite difference methods
   R. Becerril, F.S. Guzmán, A. Rendón-Romero, S. Valdez-Alvarado
   Revista Mexicana de Física E, Vol. 54, No. 2, pp. 120-132, 2008

3. A Python Program for Solving Schrödinger's Equation in Undergraduate Physical Chemistry
   Matthew N. Srnec, Shiv Upadhyay, Jeffry D. Madura
   Journal of Chemical Education, Vol 94/Issue 6, 2017