quscope.ctem

Classical CTEM reference implementations used to validate the quantum circuits in quscope.quantum_ctem.

Kirkland Atomic Potential Module

Implements the 2D projected atomic potential from Earl J. Kirkland’s “Advanced Computing in Electron Microscopy”, 2nd Edition, Appendix C.

The potential is parameterized as:

V(r) = Σᵢ[4π²aᵢ·K₀(2πr√bᵢ)] + Σᵢ[2π^(3/2)·cᵢ/dᵢ^(3/2)·exp(-π²r²/dᵢ)]

where K₀ is the modified Bessel function of the second kind.

Parameters (aᵢ, bᵢ, cᵢ, dᵢ) are tabulated in Kirkland Appendix C for each element. The potential is scaled by 14.4 to convert to electron volts (eV).

Reference:

Kirkland, E. J. (2010). Advanced Computing in Electron Microscopy (2nd ed.). Springer. ISBN 978-1-4419-6532-5

class quscope.ctem.kirkland_potential.KirklandPotential(params_file: str | Path | None = None)[source]

Bases: object

Calculate 2D projected atomic potential using Kirkland parameterization.

This class provides methods to compute the electrostatic potential of atoms projected along the beam direction (z-axis) at 2D positions (x, y).

params_dict

Dictionary of Kirkland parameters {element: [[a], [b], [c], [d]]}

V_scaling

Scaling factor to convert potential to eV (default 14.4)

Example

>>> from quscope.ctem import KirklandPotential
>>> pot = KirklandPotential()
>>>
>>> # Calculate potential for a single carbon atom
>>> x = np.linspace(-5, 5, 100)
>>> y = np.zeros_like(x)
>>> V = pot.calculate_2d(x, y, atom_x=0, atom_y=0, Z=6)
>>>
>>> # Calculate potential for multiple atoms
>>> positions = [(0, 0, 6), (5, 0, 14), (10, 0, 29)]  # (x, y, Z)
>>> V_total = pot.calculate_multiple_atoms(x_grid, y_grid, positions)
DEFAULT_PARAMS_FILE = PosixPath('/home/docs/checkouts/readthedocs.org/user_builds/quscope/checkouts/v0.2.0/src/quscope/ctem/kirkland.json')
V_SCALING = 14.4
get_element_symbol(Z: int) → str[source]

Convert atomic number to element symbol.

Parameters:

Z – Atomic number

Returns:

Element symbol (e.g., ‘C’ for Z=6)

Raises:

ValueError – If element not in parameters dictionary

calculate_2d(x_grid: ndarray, y_grid: ndarray, atom_x: float, atom_y: float, Z: int) → ndarray[source]

Calculate 2D projected atomic potential at grid points.

This implements Kirkland’s parameterization (Appendix C, Equations 5.8-5.10):

V(r) = Σᵢ[4π²aᵢ·K₀(2πr√bᵢ)] + Σᵢ[2π^(3/2)·cᵢ/dᵢ^(3/2)·exp(-π²r²/dᵢ)]

where: - r = √[(x-x_atom)² + (y-y_atom)²] is distance from atom center - K₀ is modified Bessel function of second kind, order 0 - aᵢ, bᵢ, cᵢ, dᵢ are element-specific parameters from Kirkland Appendix C

Parameters:

  • x_grid – X coordinates (Angstroms), shape (M,) or (M, N)

  • y_grid – Y coordinates (Angstroms), shape (M,) or (M, N)

  • atom_x – X position of atom center (Angstroms)

  • atom_y – Y position of atom center (Angstroms)

  • Z – Atomic number

Returns:

Potential values in eV, same shape as x_grid

Raises:

ValueError – If element not available or grid shapes don’t match

Notes

  • For r → 0, uses limiting form: K₀(x) ≈ -ln(x/2) - γ where γ = 0.5772156649

  • For large arguments (x > 50), uses asymptotic form: K₀(x) ≈ √(π/2x)·exp(-x)

  • Final potential is scaled by 14.4 to convert to eV

calculate_multiple_atoms(x_grid: ndarray, y_grid: ndarray, atom_positions: list[Tuple[float, float, int]]) → ndarray[source]

Calculate total potential from multiple atoms.

Superposition principle: total potential is sum of individual atomic potentials.

Parameters:

  • x_grid – X coordinates (Angstroms)

  • y_grid – Y coordinates (Angstroms)

  • atom_positions – List of (x, y, Z) tuples for each atom

Returns:

Total potential in eV

Example

>>> positions = [(0, 0, 6), (5, 0, 14), (10, 0, 29)]  # C, Si, Cu
>>> V_total = pot.calculate_multiple_atoms(X, Y, positions)
get_parameters(Z: int) → Dict[str, ndarray][source]

Get Kirkland parameters for an element.

Parameters:

Z – Atomic number

Returns:

Dictionary with keys ‘a’, ‘b’, ‘c’, ‘d’ containing parameter arrays

Example

>>> params = pot.get_parameters(6)  # Carbon
>>> print(f"a parameters: {params['a']}")

Multislice CTEM Simulator - Classical Implementation

This module implements the multislice algorithm for thick specimen CTEM simulation following Kirkland’s “Advanced Computing in Electron Microscopy” 2nd Edition, Chapter 7.

The multislice method divides the specimen into thin slices and alternately applies: 1. Transmission function: t_n(x,y) = exp(iσV_n(x,y)) 2. Propagation function: P(k) = exp(-iπλk²Δz)

This is the classical implementation using standard FFT (not quantum QFT).

References

Kirkland, E.J. (2010). Advanced Computing in Electron Microscopy (2nd ed.). Springer. Chapter 7: Multislice. - Figure 7.2: Magnitude of wave function at different thicknesses - Figure 7.3: Intensity and phase vs thickness curves - Figure 7.4: BF phase contrast images - Table 7.2: Mean intensity validation values

class quscope.ctem.multislice_simulator.MultisliceSimulator(image_size: float, pixels: int, beam_energy: float, slice_thickness: float = 2.0, kirkland_params_path: str | None = None)[source]

Bases: object

Classical multislice CTEM simulator for thick specimens.

The multislice algorithm simulates electron scattering through thick specimens by alternating between transmission (phase grating approximation) and propagation (Fresnel diffraction) through thin slices.

image_size

Real-space image size in Angstroms

Type:

float

pixels

Number of pixels (must be power of 2)

Type:

int

beam_energy

Electron beam energy in eV

Type:

float

slice_thickness

Thickness of each slice in Angstroms

Type:

float

wavelength

Relativistic electron wavelength (Å)

Type:

float

sigma

Interaction parameter σ (rad/(eV·Å))

Type:

float

dx

Pixel size in real space (Å)

Type:

float

Example

>>> sim = MultisliceSimulator(
...     image_size=40.0,
...     pixels=256,
...     beam_energy=200e3,
...     slice_thickness=2.0
... )
>>> results = sim.simulate_thickness_series(
...     atoms,
...     thicknesses=[40, 80, 200],
...     defocus=0
... )
kirkland_potential_2d(x_grid: ndarray, y_grid: ndarray, atom_x: float, atom_y: float, Z: int) → ndarray[source]

Calculate 2D projected atomic potential using Kirkland parameterization.

V(x,y) = Σᵢ 4π²aᵢ K₀(2πr√bᵢ) + Σᵢ π^(3/2) cᵢ/dᵢ^(3/2) exp(-π²r²/dᵢ)

where K₀ is the modified Bessel function of the second kind.

Parameters:

  • x_grid – 2D array of x coordinates (Å)

  • y_grid – 2D array of y coordinates (Å)

  • atom_x – Atom x position (Å)

  • atom_y – Atom y position (Å)

  • Z – Atomic number

Returns:

2D array of potential values in eV·Å

References

Kirkland Table 4.1: Scattering factor parameterization

get_atoms_in_slice(atoms: List[Dict], z_min: float, z_max: float) → List[Dict][source]

Get atoms within a z-range [z_min, z_max).

Parameters:

  • atoms – List of atom dictionaries with keys: - ‘position’: [x, y, z] in Angstroms - ‘Z’: Atomic number - ‘element’: Element symbol (optional)

  • z_min – Minimum z coordinate (Å)

  • z_max – Maximum z coordinate (Å)

Returns:

List of atoms within the z-range

calculate_slice_transmission(atoms_in_slice: List[Dict], slice_thickness: float) → ndarray[source]

Calculate transmission function for a slice.

t_n(x,y) = exp(iσV_n(x,y)Δz)

where V_n is the projected potential for slice n.

Parameters:

  • atoms_in_slice – List of atoms in this slice

  • slice_thickness – Thickness of slice in Angstroms

Returns:

Complex transmission function array

calculate_propagator(dz: float) → ndarray[source]

Calculate Fresnel free-space propagator.

P(k) = exp(-iπλk²Δz)

where k² = k_x² + k_y² is the squared spatial frequency.

Parameters:

dz – Propagation distance (slice thickness) in Angstroms

Returns:

Complex propagator array in reciprocal space

References

Kirkland Eq. 7.2: Fresnel propagation

simulate_thickness(atoms: List[Dict], thickness: float, defocus: float = 0, Cs: float = 0) → Dict[str, ndarray][source]

Simulate CTEM image for a single specimen thickness.

Algorithm: 1. Initialize ψ = 1 (plane wave) 2. For each slice:

  1. Apply transmission: ψ → t_n(x,y) × ψ

  2. FFT: ψ → ψ_k

  3. Apply propagator: ψ_k → P(k) × ψ_k

  4. IFFT: ψ_k → ψ

  1. Apply CTF (objective lens aberrations)

  2. Calculate intensity: I = |ψ|²

Parameters:

  • atoms – List of atom dictionaries

  • thickness – Total specimen thickness in Angstroms

  • defocus – Defocus value in Angstroms (default: 0)

  • Cs – Spherical aberration coefficient in Angstroms (default: 0)

Returns:

  • ‘intensity_image’: 2D intensity array

  • ’exit_wave’: Complex exit wave function

  • ’mean_intensity’: Mean intensity value

  • ’n_slices’: Number of slices used

Return type:

Dictionary with keys

simulate_thickness_series(atoms: List[Dict], thicknesses: List[float], defocus: float = 0, Cs: float = 0, verbose: bool = True) → Dict[float, Dict[str, ndarray]][source]

Simulate CTEM images for multiple specimen thicknesses.

This is useful for generating thickness series plots like Kirkland Figures 7.2-7.4 and validating against Table 7.2.

Parameters:

  • atoms – List of atom dictionaries

  • thicknesses – List of thickness values in Angstroms

  • defocus – Defocus value in Angstroms (default: 0)

  • Cs – Spherical aberration coefficient in Angstroms (default: 0)

  • verbose – Print progress messages (default: True)

Returns:

Dictionary mapping thickness → simulation results

Example

>>> results = sim.simulate_thickness_series(
...     atoms,
...     thicknesses=[40, 80, 128, 200, 256, 400, 512],
...     defocus=0
... )
>>> for t, res in results.items():
...     print(f"{t:.1f} Å: I_mean = {res['mean_intensity']:.3f}")

Weak Phase Object Approximation (WPOA) CTEM Simulator

Implements classical CTEM image simulation using the weak phase object approximation from Earl J. Kirkland’s “Advanced Computing in Electron Microscopy”, 2nd Edition, Chapter 5.

The WPOA assumes the transmission function is:

t(x,y) = exp(iσV(x,y))

where: - σ = interaction parameter (depends on beam energy) - V(x,y) = projected atomic potential (from KirklandPotential)

The simulation pipeline: 1. Calculate projected potential V(x,y) for all atoms 2. Calculate transmission function t(x,y) = exp(iσV) 3. Fourier transform to reciprocal space 4. Apply objective lens transfer function (CTF with aberrations) 5. Inverse Fourier transform back to real space 6. Calculate intensity I = |ψ|²

Reference:

Kirkland, E. J. (2010). Advanced Computing in Electron Microscopy (2nd ed.). Springer. Chapter 5, Figures 5.11-5.12

class quscope.ctem.wpoa_simulator.WPOASimulator(image_size: float = 50.0, pixels: int = 512, beam_energy: float = 200000.0, params_file: Path | None = None)[source]

Bases: object

Classical CTEM simulator using Weak Phase Object Approximation.

This simulator implements the WPOA method for conventional transmission electron microscopy (CTEM) image simulation. It handles: - Relativistic wavelength calculation - Interaction parameter σ calculation - Transmission function with weak phase object approximation - Objective lens transfer function with spherical aberration - Coherent bright field image formation

image_size

Physical size of image in Angstroms

pixels

Number of pixels (assumed square image)

beam_energy

Electron beam energy in eV

wavelength

Relativistic electron wavelength in Angstroms

sigma

Interaction parameter in rad/eV

dx

Pixel size in Angstroms

X, Y

2D coordinate grids in Angstroms

potential_calculator

KirklandPotential instance for atomic potentials

Example

>>> from quscope.ctem import WPOASimulator
>>>
>>> # Initialize simulator for 200 keV electrons
>>> sim = WPOASimulator(
...     image_size=50.0,  # 50 Angstroms
...     pixels=512,
...     beam_energy=200e3  # 200 keV
... )
>>>
>>> # Define atom positions: (x, y, Z)
>>> atoms = [
...     (-20, 0, 6),   # C at -20 Å
...     (-10, 0, 14),  # Si at -10 Å
...     (0, 0, 29),    # Cu at 0 Å
...     (10, 0, 79),   # Au at 10 Å
...     (20, 0, 92),   # U at 20 Å
... ]
>>>
>>> # Simulate CTEM image
>>> results = sim.simulate_image(
...     atom_positions=atoms,
...     defocus=700.0,      # 700 Å underfocus
...     Cs=1.3e7,           # 1.3 mm spherical aberration
...     alpha_max=10.37     # 10.37 mrad aperture
... )
>>>
>>> # Access results
>>> intensity = results['intensity']
>>> transmission = results['transmission']
>>> potential = results['potential']
M0C2 = 511000.0
HC = 12.2639
calculate_transmission_function(atom_positions: List[Tuple[float, float, int]]) → Tuple[ndarray, ndarray][source]

Calculate transmission function for given atom positions.

Steps: 1. Calculate total projected potential V(x,y) from all atoms 2. Calculate phase shift: φ(x,y) = σ·V(x,y) 3. Calculate transmission: t(x,y) = exp(i·φ)

Parameters:

atom_positions – List of (x, y, Z) tuples where: - x, y: atom coordinates in Angstroms - Z: atomic number

Returns:

Complex transmission function, shape (pixels, pixels) potential: Total projected potential in eV, shape (pixels, pixels)

Return type:

transmission

Example

>>> atoms = [(0, 0, 6), (5, 0, 14)]  # C and Si
>>> t, V = sim.calculate_transmission_function(atoms)
>>> np.abs(t)  # Magnitude (always 1 in WPOA)
array([[1., 1., ...], ...])
>>> np.angle(t)  # Phase shift
array([[0.12, 0.15, ...], ...])
objective_lens_transfer_function(kx: ndarray, ky: ndarray, defocus: float, Cs: float, alpha_max: float | None = None) → ndarray[source]

Calculate objective lens transfer function (CTF).

Implements Kirkland Equation 5.27:

H(k) = exp(-i·χ(k))

where the phase aberration is:

χ(k) = π·λ·k²·(½·Cs·λ²·k² - Δf)

Parameters:

  • kx – X spatial frequencies in 1/Angstrom, shape (M, N)

  • ky – Y spatial frequencies in 1/Angstrom, shape (M, N)

  • defocus – Defocus value Δf in Angstroms (positive = underfocus)

  • Cs – Spherical aberration coefficient in Angstroms

  • alpha_max – Optional aperture semi-angle in radians

Returns:

Transfer function H(k), complex array shape (M, N)

Notes

  • Positive defocus = underfocus (weaker lens)

  • Negative defocus = overfocus (stronger lens)

  • Scherzer defocus: Δf = -1.2·(Cs·λ)^(1/2)

  • If alpha_max provided, applies hard aperture cutoff

Example

>>> kx = np.fft.fftfreq(512, d=0.1)
>>> ky = np.fft.fftfreq(512, d=0.1)
>>> KX, KY = np.meshgrid(kx, ky)
>>> H = sim.objective_lens_transfer_function(
...     KX, KY, defocus=700, Cs=1.3e7, alpha_max=0.01037
... )
simulate_image(atom_positions: List[Tuple[float, float, int]], defocus: float = 700.0, Cs: float = 13000000.0, alpha_max: float | None = None, return_wavefunction: bool = False) → Dict[str, ndarray][source]

Simulate CTEM image using weak phase object approximation.

Full simulation pipeline: 1. Calculate transmission function t(x,y) 2. Fourier transform: ψ(k) = FFT[t(x,y)] 3. Apply lens CTF: ψ’(k) = ψ(k)·H(k) 4. Inverse transform: ψ(x,y) = IFFT[ψ’(k)] 5. Calculate intensity: I(x,y) = |ψ(x,y)|²

Parameters:

  • atom_positions – List of (x, y, Z) atom coordinates and atomic numbers

  • defocus – Defocus in Angstroms (positive = underfocus)

  • Cs – Spherical aberration in Angstroms (e.g., 1.3e7 = 1.3 mm)

  • alpha_max – Aperture semi-angle in milliradians (converted internally)

  • return_wavefunction – If True, include complex wavefunction in results

Returns:

  • ‘intensity’: Image intensity I = |ψ|², shape (pixels, pixels)

  • ’transmission’: Complex transmission function

  • ’potential’: Projected atomic potential in eV

  • ’positions’: Atom positions used

  • ’psi’: Complex wavefunction (if return_wavefunction=True)

Return type:

Dictionary containing

Example

>>> # Reproduce Kirkland Figure 5.12
>>> atoms = [(-20,0,6), (-10,0,14), (0,0,29), (10,0,79), (20,0,92)]
>>> results = sim.simulate_image(
...     atom_positions=atoms,
...     defocus=700.0,      # 700 Å underfocus
...     Cs=1.3e7,           # 1.3 mm
...     alpha_max=10.37     # 10.37 mrad
... )
>>> intensity = results['intensity']
>>> # Intensity range should be ~0.72 to 1.03 (Kirkland value)
get_simulation_parameters() → Dict[str, float][source]

Get current simulation parameters.

Returns:

  • beam_energy: in eV

  • wavelength: in Angstroms

  • sigma: interaction parameter in rad/eV

  • pixel_size: in Angstroms

  • image_size: in Angstroms

  • pixels: number of pixels

Return type:

Dictionary with simulation parameters

Example

>>> params = sim.get_simulation_parameters()
>>> print(f"Wavelength: {params['wavelength']:.5f} Å")
Wavelength: 0.02508 Å