Quantum Mechanics Programs

Interactive Python simulations for quantum mechanics — inspired by qmsolve

These simulations run entirely in your browser using Pyodide (WebAssembly Python). For advanced 3D visualizations and GPU-accelerated computations, install qmsolve locally.

These programs use the same numerical methods as qmsolve (split-step Fourier, eigenstate solvers) but run directly in your browser. The first run downloads the Python environment (~15MB). Click "Run" to execute!

1D Harmonic Oscillator

Solve the quantum harmonic oscillator analytically. Visualize eigenstates, probability densities, and energy levels. Uses Hermite polynomials and demonstrates the equally-spaced energy spectrum E_n = (n+1/2)hbar*omega.

import numpy as np
import matplotlib.pyplot as plt
from math import factorial

# 1D Quantum Harmonic Oscillator - Eigenstates
# Inspired by qmsolve (github.com/quantum-visualizations/qmsolve)

# Hermite polynomials via recurrence: H_0=1, H_1=2x, H_{n+1}=2x*H_n - 2n*H_{n-1}
def hermite_val(n, x):
    """Evaluate physicist's Hermite polynomial H_n(x) using recurrence"""
    if n == 0:
        return np.ones_like(x)
    elif n == 1:
        return 2.0 * x
    H_prev, H_curr = np.ones_like(x), 2.0 * x
    for k in range(1, n):
        H_prev, H_curr = H_curr, 2.0 * x * H_curr - 2.0 * k * H_prev
    return H_curr

# Physical parameters (in natural units: hbar = m = omega = 1)
hbar = 1.0
m = 1.0
omega = 1.0

# Spatial grid
x = np.linspace(-6, 6, 500)

# Analytic wavefunctions: psi_n(x) = N_n * H_n(alpha*x) * exp(-alpha^2 x^2 / 2)
# where alpha = sqrt(m*omega/hbar), N_n = normalization, H_n = Hermite polynomial

alpha = np.sqrt(m * omega / hbar)

def psi_n(n, x):
    """Quantum harmonic oscillator eigenstate n"""
    H_n = hermite_val(n, alpha * x)
    N_n = (alpha / np.pi**0.5 / (2**n * factorial(n)))**0.5
    return N_n * H_n * np.exp(-alpha**2 * x**2 / 2)

def E_n(n):
    """Energy of eigenstate n"""
    return hbar * omega * (n + 0.5)

# Potential
V = 0.5 * m * omega**2 * x**2

# Plot first 6 eigenstates
n_states = 6
fig, axes = plt.subplots(1, 3, figsize=(15, 6))

# --- Left: Wavefunctions offset by energy ---
ax1 = axes[0]
colors = plt.cm.viridis(np.linspace(0.2, 0.9, n_states))

ax1.plot(x, V, 'w--', alpha=0.4, linewidth=1.5, label='V(x)')
for n in range(n_states):
    psi = psi_n(n, x)
    energy = E_n(n)
    ax1.fill_between(x, energy, energy + psi * 1.5, alpha=0.3, color=colors[n])
    ax1.plot(x, energy + psi * 1.5, color=colors[n], linewidth=1.5, label=f'n={n}, E={energy:.1f}')
    ax1.axhline(y=energy, color=colors[n], linestyle=':', alpha=0.3)

ax1.set_xlim(-5, 5)
ax1.set_ylim(-0.5, n_states * hbar * omega + 1)
ax1.set_xlabel('Position x', color='white', fontsize=11)
ax1.set_ylabel('Energy / Wavefunction', color='white', fontsize=11)
ax1.set_title('Wavefunctions (offset by energy)', fontsize=13, color='#c084fc')
ax1.legend(fontsize=7, facecolor='#1e293b', edgecolor='gray', labelcolor='white', loc='upper right')
ax1.set_facecolor('#1e293b')
ax1.tick_params(colors='white')
for spine in ax1.spines.values():
    spine.set_color('gray')
ax1.grid(True, alpha=0.2)

# --- Center: Probability densities ---
ax2 = axes[1]
for n in range(n_states):
    psi = psi_n(n, x)
    prob = np.abs(psi)**2
    ax2.plot(x, prob, color=colors[n], linewidth=1.8, label=f'|psi_{n}|^2')

ax2.set_xlim(-5, 5)
ax2.set_xlabel('Position x', color='white', fontsize=11)
ax2.set_ylabel('Probability Density |psi|^2', color='white', fontsize=11)
ax2.set_title('Probability Densities', fontsize=13, color='#c084fc')
ax2.legend(fontsize=8, facecolor='#1e293b', edgecolor='gray', labelcolor='white')
ax2.set_facecolor('#1e293b')
ax2.tick_params(colors='white')
for spine in ax2.spines.values():
    spine.set_color('gray')
ax2.grid(True, alpha=0.2)

# --- Right: Energy level diagram ---
ax3 = axes[2]
for n in range(n_states):
    energy = E_n(n)
    ax3.hlines(energy, 0, 1, colors=colors[n], linewidth=3)
    ax3.text(1.1, energy, f'n={n}: E = {energy:.1f} hbar*omega', color=colors[n], fontsize=10, va='center')

ax3.set_xlim(-0.5, 3.5)
ax3.set_ylim(-0.5, n_states * hbar * omega + 1)
ax3.set_ylabel('Energy (hbar * omega)', color='white', fontsize=11)
ax3.set_title('Energy Levels', fontsize=13, color='#c084fc')
ax3.set_facecolor('#1e293b')
ax3.tick_params(colors='white')
ax3.set_xticks([])
for spine in ax3.spines.values():
    spine.set_color('gray')

plt.tight_layout()
plt.show()

print("Quantum Harmonic Oscillator - Eigenstates")
print("=" * 50)
print(f"Parameters: hbar={hbar}, m={m}, omega={omega}")
print(f"alpha = sqrt(m*omega/hbar) = {alpha:.4f}")
print()
for n in range(n_states):
    print(f"  n={n}: E_{n} = {E_n(n):.2f} hbar*omega")
print()
print("Equal energy spacing: Delta E = hbar * omega = 1.0")
print("Zero-point energy: E_0 = hbar * omega / 2 = 0.5")

Click Run to execute the Python code

First run will download Python environment (~15MB)

Fundamental Equations of Quantum Mechanics

The Schrodinger Equation

The fundamental equation governing quantum mechanical systems:

Time-Dependent Schrodinger Equation:

$$i\hbar \frac{\partial}{\partial t} \Psi(\mathbf{r}, t) = \hat{H} \Psi(\mathbf{r}, t)$$

Time-Independent Schrodinger Equation:

$$\hat{H} \psi(\mathbf{r}) = E \psi(\mathbf{r})$$

where $ \hat{H} = -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) $

Split-Step Fourier Method (used in tunneling & double-slit simulations):

$$\Psi(t+\Delta t) \approx e^{-i\hat{V}\Delta t/2\hbar} \cdot \mathcal{F}^{-1}\!\left[e^{-i\hbar k^2\Delta t/2m}\cdot\mathcal{F}\!\left[e^{-i\hat{V}\Delta t/2\hbar}\Psi(t)\right]\right]$$

Harmonic Oscillator

$$E_n = \hbar\omega\!\left(n + \tfrac{1}{2}\right), \quad n = 0, 1, 2, \ldots$$

Wavefunctions:

$$\psi_n(x) = \left(\frac{\alpha}{\sqrt{\pi}\,2^n n!}\right)^{1/2} H_n(\alpha x)\,e^{-\alpha^2 x^2/2}$$

Infinite Square Well

$$E_n = \frac{n^2 \pi^2 \hbar^2}{2mL^2}, \quad n = 1, 2, 3, \ldots$$

Wavefunctions:

$$\psi_n(x) = \sqrt{\frac{2}{L}}\,\sin\!\left(\frac{n\pi x}{L}\right)$$

Hydrogen Atom

Energy levels:

$$E_n = -\frac{13.6 \text{ eV}}{n^2}$$

Radial wavefunction:

$$R_{nl}(r) = N_{nl}\,e^{-r/na_0}\left(\frac{2r}{na_0}\right)^l L_{n-l-1}^{2l+1}\!\left(\frac{2r}{na_0}\right)$$

Quantum numbers:

$ n = 1, 2, 3, \ldots $ — principal quantum number

$ l = 0, 1, \ldots, n{-}1 $ — angular momentum

$ m = -l, \ldots, +l $ — magnetic quantum number

Tunneling probability (rectangular barrier):

$$T \approx e^{-2\kappa a}, \quad \kappa = \frac{\sqrt{2m(V_0 - E)}}{\hbar}$$

About qmsolve

The simulations on this page are inspired by qmsolve, a Python module for solving and visualizing the Schrodinger equation. For more advanced capabilities, install it locally:

# Basic installation

pip install qmsolve

# With 3D volumetric visualization (Mayavi)

pip install qmsolve[with_mayavi]

Eigenstate Solver

Solves for eigenstates of 1D, 2D, and 3D systems using ARPACK's Lanczos method and LOBPCG. Handles single and two-particle systems.

Time-Dependent Solver

Split-step Fourier and Cayley-Crank-Nicolson methods for time evolution. Supports complex potentials and momentum-dependent interactions.

3D Visualization

Volumetric rendering via Mayavi, interactive eigenstate exploration with sliders, and GPU acceleration through CuPy/CUDA.

Quick Example (local installation)

from qmsolve import Hamiltonian, SingleParticle, TimeSimulation
from qmsolve.util.constants import eV, Å

def harmonic_oscillator(particle):
    k = 100. * eV / Å**2
    return 0.5 * k * particle.x**2

H = Hamiltonian(
    particles=SingleParticle(),
    potential=harmonic_oscillator,
    spatial_ndim=1, N=512, extent=20*Å
)

eigenstates = H.solve(max_states=30)
eigenstates.visualize(show=True)

Run It Here: qmsolve-Style Eigenstate Solver

This self-contained implementation reproduces qmsolve's eigenstate solver using finite-difference discretization of the Hamiltonian. It solves the 1D time-independent Schrödinger equation for three potentials: harmonic oscillator, double well, and Morse potential — computing eigenstates, energies, probability densities, and expectation values.

Python

script.py154 lines

import numpy as np
from math import factorial

# ============================================================
#  qmsolve-style Eigenstate Solver (self-contained)
#  Solves: -hbar^2/(2m) d^2psi/dx^2 + V(x)psi = E*psi
#  Method: Finite-difference Hamiltonian + numpy eigensolver
# ============================================================

# --- Physical constants (atomic units: hbar = m_e = e = 1) ---
hbar = 1.0
m_e  = 1.0

# --- Grid setup ---
N = 300            # grid points
x_min, x_max = -10.0, 10.0
x  = np.linspace(x_min, x_max, N)
dx = x[1] - x[0]

# --- Build kinetic energy matrix (tridiagonal) ---
# T_ij = -hbar^2/(2m) * d^2/dx^2  (finite difference)
diag     =  np.ones(N) * (-2.0)
off_diag =  np.ones(N - 1)
T = np.zeros((N, N))
for i in range(N):
    T[i, i] = -2.0
    if i > 0:     T[i, i-1] = 1.0
    if i < N - 1: T[i, i+1] = 1.0
T *= -hbar**2 / (2.0 * m_e * dx**2)

# --- Define potentials ---
def harmonic(x, omega=1.0):
    return 0.5 * m_e * omega**2 * x**2

def double_well(x, a=0.05, b=1.0):
    return a * (x**2 - b * 5)**2

def morse(x, D=5.0, alpha=0.5, x0=0.0):
    return D * (1 - np.exp(-alpha * (x - x0)))**2

potentials = [
    ("Harmonic Oscillator (omega=1)", harmonic(x)),
    ("Double Well",                   double_well(x)),
    ("Morse Potential (D=5)",         morse(x)),
]

n_states = 8  # number of eigenstates to display

print("=" * 70)
print("  qmsolve-style 1D Eigenstate Solver")
print("  Grid: N=%d points, x in [%.1f, %.1f], dx=%.4f" % (N, x_min, x_max, dx))
print("=" * 70)

for pot_name, V in potentials:
    # Build full Hamiltonian: H = T + V
    H = T + np.diag(V)

# Solve eigenvalue problem (only need lowest states)
    eigenvalues, eigenvectors = np.linalg.eigh(H)

# Extract lowest n_states
    energies = eigenvalues[:n_states]
    states   = eigenvectors[:, :n_states]

# Normalize wavefunctions
    for j in range(n_states):
        norm = np.sqrt(np.trapezoid(states[:, j]**2, x))
        states[:, j] /= norm

print("\n" + "-" * 70)
    print("  Potential: %s" % pot_name)
    print("-" * 70)
    print("  %-6s  %12s  %12s  %12s  %12s" % (
        "State", "Energy", "E_exact*", "<x>", "<x^2>"))
    print("  " + "-" * 58)

for j in range(n_states):
        psi = states[:, j]
        prob = psi**2

# Expectation values
        x_mean  = np.trapezoid(x * prob, x)
        x2_mean = np.trapezoid(x**2 * prob, x)

# Exact energy for harmonic oscillator
        if "Harmonic" in pot_name:
            E_exact = hbar * 1.0 * (j + 0.5)
            print("  n=%-4d  %12.6f  %12.6f  %12.6f  %12.6f" % (
                j, energies[j], E_exact, x_mean, x2_mean))
        else:
            print("  n=%-4d  %12.6f  %12s  %12.6f  %12.6f" % (
                j, energies[j], "---", x_mean, x2_mean))

if "Harmonic" in pot_name:
        print("\n  * E_exact = hbar*omega*(n+1/2) for harmonic oscillator")
        max_err = max(abs(energies[j] - hbar*(j+0.5)) for j in range(n_states))
        print("  * Max |E_numerical - E_exact| = %.2e" % max_err)

# Uncertainty principle check
    print("\n  Uncertainty principle check (sigma_x * sigma_p >= hbar/2):")
    for j in range(min(4, n_states)):
        psi = states[:, j]
        prob = psi**2
        x_mean  = np.trapezoid(x * prob, x)
        x2_mean = np.trapezoid(x**2 * prob, x)
        sigma_x = np.sqrt(x2_mean - x_mean**2)

# Compute <p^2> via -hbar^2 * integral(psi * d^2psi/dx^2)
        d2psi = np.gradient(np.gradient(psi, dx), dx)
        p2_mean = -hbar**2 * np.trapezoid(psi * d2psi, x)
        sigma_p = np.sqrt(abs(p2_mean))

product = sigma_x * sigma_p
        print("    n=%d: sigma_x=%.4f, sigma_p=%.4f, product=%.4f (>= %.4f)" % (
            j, sigma_x, sigma_p, product, hbar / 2))

# --- Wavefunction properties at the boundaries ---
print("\n" + "=" * 70)
print("  Wavefunction Decay at Boundaries")
print("=" * 70)
V_ho = harmonic(x)
H_ho = T + np.diag(V_ho)
evals, evecs = np.linalg.eigh(H_ho)
for j in range(4):
    psi = evecs[:, j]
    norm = np.sqrt(np.trapezoid(psi**2, x))
    psi /= norm
    print("  n=%d: |psi(x_min)|=%.2e, |psi(x_max)|=%.2e  (should be ~0)" % (
        j, abs(psi[0]), abs(psi[-1])))

# --- Transition dipole moments ---
print("\n" + "=" * 70)
print("  Transition Dipole Moments <m|x|n> (Harmonic Oscillator)")
print("=" * 70)
print("  Selection rule: <m|x|n> != 0 only if m = n +/- 1")
print()
print("  %6s" % "", end="")
for n in range(6):
    print("  n=%-6d" % n, end="")
print()
for mi in range(6):
    psi_m = evecs[:, mi] / np.sqrt(np.trapezoid(evecs[:, mi]**2, x))
    print("  m=%-4d" % mi, end="")
    for ni in range(6):
        psi_n = evecs[:, ni] / np.sqrt(np.trapezoid(evecs[:, ni]**2, x))
        dipole = np.trapezoid(psi_m * x * psi_n, x)
        print("  %8.4f" % dipole, end="")
    print()

print("\n  Note: Non-zero elements appear only on the super/sub-diagonal,")
print("  confirming the selection rule Delta_n = +/- 1.")
print("\n  Solver complete. %d eigenstates computed for %d potentials." % (
    n_states, len(potentials)))

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Run It Here: Fortran Eigenstate Solver

A Fortran implementation of the same finite-difference eigenstate solver — solving the 1D Schrödinger equation for a quantum harmonic oscillator using the Jacobi eigenvalue algorithm (no LAPACK required).

Fortran

program.f90212 lines

program qmsolve_fortran
  implicit none
  integer, parameter :: N = 100       ! grid points
  integer, parameter :: n_states = 6  ! eigenstates to show
  double precision, parameter :: hbar = 1.0d0, mass = 1.0d0, omega = 1.0d0
  double precision, parameter :: x_min = -8.0d0, x_max = 8.0d0
  double precision :: dx, x(N), V(N), H(N,N), eigvec(N,N)
  double precision :: energies(N), psi(N), prob(N)
  double precision :: x_mean, x2_mean, sigma_x, norm, E_exact, err
  integer :: i, j, k, n_rot

dx = (x_max - x_min) / dble(N - 1)

! Build grid and potential
  do i = 1, N
    x(i) = x_min + (i - 1) * dx
    V(i) = 0.5d0 * mass * omega**2 * x(i)**2
  end do

! Build Hamiltonian matrix: H = T + V
  H = 0.0d0
  do i = 1, N
    H(i, i) = hbar**2 / (mass * dx**2) + V(i)
    if (i > 1)   H(i, i-1) = -hbar**2 / (2.0d0 * mass * dx**2)
    if (i < N)   H(i, i+1) = -hbar**2 / (2.0d0 * mass * dx**2)
  end do

! Initialize eigenvector matrix to identity
  eigvec = 0.0d0
  do i = 1, N
    eigvec(i, i) = 1.0d0
  end do

! Jacobi eigenvalue algorithm (for symmetric matrices)
  call jacobi_eigen(H, eigvec, N, n_rot)

! Extract eigenvalues (diagonal of H after Jacobi)
  do i = 1, N
    energies(i) = H(i, i)
  end do

! Sort eigenvalues and eigenvectors (simple selection sort)
  do i = 1, n_states + 2
    k = i
    do j = i + 1, N
      if (energies(j) < energies(k)) k = j
    end do
    if (k /= i) then
      ! Swap eigenvalues
      err = energies(i)
      energies(i) = energies(k)
      energies(k) = err
      ! Swap eigenvectors
      do j = 1, N
        err = eigvec(j, i)
        eigvec(j, i) = eigvec(j, k)
        eigvec(j, k) = err
      end do
    end if
  end do

write(*,'(A)') '======================================================================'
  write(*,'(A)') '  Fortran qmsolve-style 1D Eigenstate Solver'
  write(*,'(A,I4,A,F6.2,A,F6.2,A,F6.4)') '  Grid: N=', N, &
    ', x in [', x_min, ',', x_max, '], dx=', dx
  write(*,'(A)') '======================================================================'
  write(*,'(A)') ''
  write(*,'(A)') '  Potential: Quantum Harmonic Oscillator (omega=1)'
  write(*,'(A)') '----------------------------------------------------------------------'
  write(*,'(A6,A14,A14,A14,A14)') 'State', 'E_numerical', 'E_exact', 'Error', '<x>'

do j = 1, n_states
    ! Get wavefunction and normalize
    do i = 1, N
      psi(i) = eigvec(i, j)
    end do
    norm = 0.0d0
    do i = 1, N
      norm = norm + psi(i)**2 * dx
    end do
    norm = sqrt(norm)
    psi = psi / norm

! Probability density
    prob = psi**2

! Expectation values
    x_mean = 0.0d0
    x2_mean = 0.0d0
    do i = 1, N
      x_mean  = x_mean  + x(i) * prob(i) * dx
      x2_mean = x2_mean + x(i)**2 * prob(i) * dx
    end do

E_exact = hbar * omega * (dble(j - 1) + 0.5d0)
    err = abs(energies(j) - E_exact)

write(*,'(A2,I2,A2,F12.6,F14.6,ES14.2,F14.6)') &
      'n=', j-1, '  ', energies(j), E_exact, err, x_mean
  end do

! Uncertainty principle
  write(*,'(A)') ''
  write(*,'(A)') '  Uncertainty Principle: sigma_x * sigma_p >= hbar/2 = 0.5000'
  write(*,'(A)') '  --------------------------------------------------------'
  do j = 1, min(4, n_states)
    psi = eigvec(:, j) / sqrt(sum(eigvec(:, j)**2) * dx)
    x_mean  = sum(x * psi**2) * dx
    x2_mean = sum(x**2 * psi**2) * dx
    sigma_x = sqrt(x2_mean - x_mean**2)
    write(*,'(A,I1,A,F8.4)') '    n=', j-1, ':  sigma_x = ', sigma_x
  end do

write(*,'(A)') ''
  write(*,'(A)') '  Solver complete.'

contains

subroutine jacobi_eigen(A, V, n, nrot)
    ! Jacobi eigenvalue algorithm for real symmetric matrix
    implicit none
    integer, intent(in) :: n
    double precision, intent(inout) :: A(n,n), V(n,n)
    integer, intent(out) :: nrot
    integer :: i, j, ip, iq, max_iter
    double precision :: tresh, theta, tau_val, t, s, c_val
    double precision :: g, h, sm, aip, aiq
    double precision :: b(n), z(n), d(n)
    max_iter = 100

do ip = 1, n
      b(ip) = A(ip, ip)
      d(ip) = b(ip)
      z(ip) = 0.0d0
    end do

nrot = 0
    do i = 1, max_iter
      sm = 0.0d0
      do ip = 1, n - 1
        do iq = ip + 1, n
          sm = sm + abs(A(ip, iq))
        end do
      end do
      if (sm < 1.0d-12) return

if (i < 4) then
        tresh = 0.2d0 * sm / dble(n * n)
      else
        tresh = 0.0d0
      end if

do ip = 1, n - 1
        do iq = ip + 1, n
          g = 100.0d0 * abs(A(ip, iq))
          if (i > 4 .and. abs(d(ip)) + g == abs(d(ip)) &
              .and. abs(d(iq)) + g == abs(d(iq))) then
            A(ip, iq) = 0.0d0
          else if (abs(A(ip, iq)) > tresh) then
            h = d(iq) - d(ip)
            if (abs(h) + g == abs(h)) then
              t = A(ip, iq) / h
            else
              theta = 0.5d0 * h / A(ip, iq)
              t = 1.0d0 / (abs(theta) + sqrt(1.0d0 + theta**2))
              if (theta < 0.0d0) t = -t
            end if
            c_val = 1.0d0 / sqrt(1.0d0 + t**2)
            s = t * c_val
            tau_val = s / (1.0d0 + c_val)
            h = t * A(ip, iq)
            z(ip) = z(ip) - h
            z(iq) = z(iq) + h
            d(ip) = d(ip) - h
            d(iq) = d(iq) + h
            A(ip, iq) = 0.0d0

do j = 1, ip - 1
              g = A(j, ip); h = A(j, iq)
              A(j, ip) = g - s * (h + g * tau_val)
              A(j, iq) = h + s * (g - h * tau_val)
            end do
            do j = ip + 1, iq - 1
              g = A(ip, j); h = A(j, iq)
              A(ip, j) = g - s * (h + g * tau_val)
              A(j, iq) = h + s * (g - h * tau_val)
            end do
            do j = iq + 1, n
              g = A(ip, j); h = A(iq, j)
              A(ip, j) = g - s * (h + g * tau_val)
              A(iq, j) = h + s * (g - h * tau_val)
            end do
            do j = 1, n
              g = V(j, ip); h = V(j, iq)
              V(j, ip) = g - s * (h + g * tau_val)
              V(j, iq) = h + s * (g - h * tau_val)
            end do
            nrot = nrot + 1
          end if
        end do
      end do

do ip = 1, n
        b(ip) = b(ip) + z(ip)
        d(ip) = b(ip)
        z(ip) = 0.0d0
      end do
    end do
  end subroutine jacobi_eigen

end program qmsolve_fortran

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server