Usage

Installation

To use QM sim, first install it using pip:

$ pip install qm-sim

To run calculations on a GPU, install the PyTorch version for your system at the PyTorch website as well

Examples

No potential

"""

Calculate eigenenergies and eigenstates of a 1D system with zero potential

"""
from qm_sim.hamiltonian import Hamiltonian
from qm_sim.nature_constants import m_e

N = 100         # Discretisation point count
L = 1e-9        # System size in meters
m = m_e         # Mass of the particle, here chosen as electron mass
n = 4           # The amount of eigenstates to find

# Set hamiltonian
H = Hamiltonian(N, L, m)

# Solve and plot the system
H.plot_eigen(n)

Quadratic potential

"""

Calculate eigenenergies and eigenstates of a 1D system with quadratic potential

"""
from qm_sim.hamiltonian import Hamiltonian
from qm_sim.nature_constants import m_e

import numpy as np

N = 1000        # Discretisation point count
L = 1e-9        # System size in meters
m = m_e         # Mass of the particle, here chosen as electron mass
n = 4           # The amount of eigenstates to find

# Set hamiltonian
H = Hamiltonian(N, L, m)

# Set potential
x, = H.get_coordinate_arrays()

def V(t):
    k = 200
    return 1/2*(k+t)*x**2
H.V = V

# Plot
H.plot_potential()
H.plot_eigen(n)

2D system

"""

Calculate and plot eigenstates in a 2D potential

"""
from qm_sim.hamiltonian import Hamiltonian
from qm_sim.nature_constants import m_e, e_0

import numpy as np

N = (200, 200)          # Discretisation point count
L = (8e-9, 8e-9)        # System size in meters
m = m_e                 # Mass of the particle, here chosen as electron mass
n = 4                   # The amount of eigenstates to find

# Set hamiltonian
H = Hamiltonian(N, L, m)

# Set potential. Here, a square well is used, with dV = 0.15 eV
V = 0.15*e_0 * np.ones(N)
V[
    N[0] // 2 - N[0] // 5 : N[0] // 2 + N[0] // 5, 
    N[1] // 2 - N[1] // 5 : N[1] // 2 + N[1] // 5
] = 0
H.V = V

# Plot potential
H.plot_potential()

# Plot eigenstate
H.plot_eigen(n)

Adiabatic evolution

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import FuncAnimation
from scipy.special import expit

from qm_sim.hamiltonian import Hamiltonian
from qm_sim.nature_constants import e_0, m_e

s = 100
N = (s, s)              # Discretisation point count
L = (8e-9, 8e-9)        # System size in meters
m = m_e                 # Mass of the particle, here chosen as electron mass
n = 4                   # The amount of eigenstates to find
steps = 100


H = Hamiltonian(N, L, m)

# Set potential. Here, a square well is used, with dV = 0.15 eV
X, Y = H.get_coordinate_arrays()

# Smoothed rectangular time-dependent potential, assymetric to avoid degeneracy 
def V(t):
    a = 0.15*e_0*expit(5*10**9*(X-2e-9 + 1e-11*t)) + 0.15*e_0*expit(-5*10**9*(X+2e-9-1e-11*t))
    b = 0.15*e_0*expit(5*10**9*(Y-2.5e-9 - 1e-11*t)) + 0.15*e_0*expit(-5*10**9*(Y+2.5e-9+1e-11*t))
    return np.where(a+b>0.15*e_0,0.15*e_0,a+b)

H.V = V

# Solve the system
energies, states = H.eigen(n)
# Choose energy to follow
E = energies[1]

# Evolve adiabatically
Energies, psi = H.adiabatic_evolution(E_n=E, t0=0, dt=1, steps=steps)

# Create animation
fig, (ax, ax2) = plt.subplots(1,2)
im = ax.matshow(psi[:,:,0]**2)
im2 = ax2.matshow(V(0).T)
ims = [im,im2]

def init():
    im.set_data(psi[:,:,0]**2)
    im2.set_data(V(0).T)
    return im

def update(j):
    ims[0].set_data(psi[:,:,j]**2)
    ims[1].set_data(V(j).T)
    return ims

anim = FuncAnimation(fig, func=update, init_func=init, frames=steps, interval=50, blit=False, repeat = True, repeat_delay = 1000)
plt.show()

Temporal evolution

import numpy as np

from qm_sim.hamiltonian import Hamiltonian
from qm_sim.nature_constants import e_0, m_e


N = 200                 # Discretisation point count
L = 2e-9                # System size in meters
m = m_e                 # Mass of the particle, here chosen as electron mass
t_end = 10e-15          # Simulation time
dt = 1e-17              # Simulation data storage stepsize

# Initialize the hamiltonian
H = Hamiltonian(N, L, m, temporal_scheme="leapfrog")

# Set the potential to a quadratic potential oscilating from side to side
z, = H.get_coordinate_arrays()
H.V = lambda t: 6*z**2 + 3*z*np.abs(z)*np.sin(4e15*t)

# Plot
H.plot_temporal(t_end, dt)

Temporal evolution of a 2D system

import numpy as np

from qm_sim.hamiltonian import Hamiltonian
from qm_sim.nature_constants import e_0, m_e
from scipy.signal import convolve2d

N = (100, 100)          # Discretisation point count
L = (15e-9, 15e-9)      # System size in meters
m = m_e                 # Mass of the particle, here chosen as electron mass
t_end = 100e-15         # Simulation time
dt = 1e-16              # Simulation data storage stepsize

# Initialize the hamiltonian.
H = Hamiltonian(N, L, m, temporal_scheme="scipy-Runge-Kutta 3(2)")

# Set the potential. Use a cool elipse cross that rotates over time
a = 2e-9
b = 5e-9
X, Y = H.get_coordinate_arrays()
def V(theta: float) -> np.ndarray:
    ct, st = np.cos(theta), np.sin(theta)

    V = np.ones(N) * e_0
    V[((X*ct + Y*st) / a)**2 + ((X*st - Y*ct) / b)**2 <= 1] = 0
    V[((X*ct + Y*st) / b)**2 + ((X*st - Y*ct) / a)**2 <= 1] = 0

    # Smooth out transition 
    n = 5
    kernel = np.ones((n,n)) / n**2
    V = convolve2d(V, kernel, mode="same", boundary="wrap")
    return V
    
H.V = lambda t: V(1e14*t)

# Set initial condition to a
# superposition of third and fourth eigenstate
_, psi = H.eigen(4)
psi_0 = 1/2**0.5 * (psi[2] + psi[3])

# Plot
H.plot_potential()
H.plot_eigen(4)
H.plot_temporal(t_end, dt, psi_0)

Hydrogen atom

# NOTE: this example uses skimage, available with:
# `pip install scikit-image`
from skimage import measure

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from qm_sim.hamiltonian import Hamiltonian
from qm_sim.nature_constants import e_0, 𝜀_0, m_e, m_p

N = (60, 60, 60)            # Discretisation point count
L = (1e-9, 1e-9, 1e-9)      # System size in meters
m = m_e * m_p / (m_e + m_p) # Mass of the particle, here the reduced mass of the system

nx, ny = 2, 2               # Used for the plot
n = nx * ny                 # Number of orbitals to find

# Initialize the hamiltonian.
H = Hamiltonian(N, L, m)

# Define the potential: -e^2 / (4*pi*𝜀_0*r)
x, y, z = H.get_coordinate_arrays()
r = (x**2 + y**2 + z**2)**0.5
H.V = -e_0**2 / (4 * np.pi * 𝜀_0 * r)

# Find a couple eigenvalues
energies, orbitals = H.eigen(n)
orbitals = np.abs(orbitals**2)

# Plot a iso-surface of each orbital
plt.figure()
plots = []
for i in range(n):
    iso_val = np.mean(orbitals[i])
    vertices, faces, _, _ = measure.marching_cubes(orbitals[i], iso_val, spacing=L)
    ax = plt.subplot(ny, nx, i+1, projection='3d')
    plots.append(ax.plot_trisurf(vertices[:, 0], vertices[:,1], faces, vertices[:, 2],
                    cmap='Spectral', lw=1))
    plt.title(f"E = {energies[i] :.2e}")
plt.tight_layout()
plt.show()