atomiongpu.m: GPU-Accelerated Atom-Ion Dynamics

• atomiongpu.m is a specialized MATLAB tool for simulating classical atom-ion collision dynamics with high throughput via GPU acceleration.
• It employs a custom ODE integrator, ode45gpu, using a scalar-optimized Runge-Kutta 4/5 algorithm that achieves significant speedups over MATLAB’s ode45.
• The script supports various trap configurations and collision ensembles, computing key observables such as scattering angles, complex lifetimes, and momentum transfers.

Searching arXiv for the specified paper and closely related GPU atom/ion simulation work. atomiongpu.m is a MATLAB user-facing script for simulating the classical dynamics of a single atom colliding with a trapped ion in the low-density regime, where one assumes the ion interacts with only one atom at a time. It is presented in “GPU-Accelerated MATLAB Software for Atom-Ion Dynamics” (Chowdhury et al., 15 Sep 2025). The script is built on a custom ODE integrator, ode45gpu, and is designed to run very large numbers of independent atom-ion trajectories—up to roughly 10 million at once—using either CPUs or GPUs. Its purpose is to automate the workflow needed to compute statistical observables relevant to cold and ultracold atom-ion experiments, including complex formation probability, scattering-angle distributions, complex lifetimes, and selected full trajectory plots (Chowdhury et al., 15 Sep 2025).

1. Definition and physical scope

atomiongpu.m models the classical dynamics of a system consisting of one ion, optionally untrapped, in a harmonic trap, or in a Paul trap; one neutral atom approaching from far away; and a low-density atom cloud, represented statistically by repeating many independent atom-ion encounters with different initial conditions (Chowdhury et al., 15 Sep 2025). The script ignores internal quantum degrees of freedom and treats the motion classically, integrating Newton’s equations for the 6 spatial degrees of freedom of the atom-ion pair. Because the ion is trapped, center-of-mass and relative motion generally do not separate, so the simulation evolves the full 12-dimensional state vector (Chowdhury et al., 15 Sep 2025).

The atom-ion interaction is assumed to be a generalized long-range-plus-short-range potential of the form

$V(r) = -\frac{C_n}{r^n} + \frac{C_m}{r^m},$

with the benchmark case using

$V_{\mathrm{ai}(r_{\mathrm{ai}) = -\frac{C_4}{r_{\mathrm{ai}^4} + \frac{C_8}{r_{\mathrm{ai}^8}.$

This interaction gives rise to elastic scattering, momentum transfer, and, for some initial conditions, long-lived temporary atom-ion complexes (Chowdhury et al., 15 Sep 2025).

The script is therefore a high-throughput classical trajectory engine for trapped atom-ion collision studies. Its central design is not a new physical model, but a specialized implementation for scanning large atom-ion collision ensembles and extracting observables from them (Chowdhury et al., 15 Sep 2025).

2. Dynamical formulation

For the benchmark system—an Rb atom colliding with a trapped Yb$^+$ ion—the Hamiltonian is written as

$$H= \frac{1}{2}m_av_a^2+\frac{1}{2}m_iv_i^2-\frac{C_4}{r_\textrm{ai}^4}+\frac{C_8}{r_\textrm{ai}^8} + \sum\limits_{i=1}^{3}(a_i+2q_i\cos\Omega_\textrm{rf}t)\frac{m_\textrm{ion}\Omega_\textrm{rf}^2}{8}r_i^2 .$$

Here $m_a, v_a$ are atom mass and speed; $m_i=m_{\rm ion}, v_i$ are ion mass and speed; $r_{\rm ai}$ is atom-ion separation; $C_4$ is the induced-dipole coefficient; $C_8$ is the short-range coefficient; $a_i, q_i$ are Paul-trap parameters; and $V_{\mathrm{ai}(r_{\mathrm{ai}) = -\frac{C_4}{r_{\mathrm{ai}^4} + \frac{C_8}{r_{\mathrm{ai}^8}.$0 is the trap drive frequency (Chowdhury et al., 15 Sep 2025).

The dynamical variables are assembled into a 12-component state vector $V_{\mathrm{ai}(r_{\mathrm{ai}) = -\frac{C_4}{r_{\mathrm{ai}^4} + \frac{C_8}{r_{\mathrm{ai}^8}.$1, with $V_{\mathrm{ai}(r_{\mathrm{ai}) = -\frac{C_4}{r_{\mathrm{ai}^4} + \frac{C_8}{r_{\mathrm{ai}^8}.$2 the ion position, $V_{\mathrm{ai}(r_{\mathrm{ai}) = -\frac{C_4}{r_{\mathrm{ai}^4} + \frac{C_8}{r_{\mathrm{ai}^8}.$3 the ion velocity, $V_{\mathrm{ai}(r_{\mathrm{ai}) = -\frac{C_4}{r_{\mathrm{ai}^4} + \frac{C_8}{r_{\mathrm{ai}^8}.$4 the atom position, and $V_{\mathrm{ai}(r_{\mathrm{ai}) = -\frac{C_4}{r_{\mathrm{ai}^4} + \frac{C_8}{r_{\mathrm{ai}^8}.$5 the atom velocity. The ODE is written as

$V_{\mathrm{ai}(r_{\mathrm{ai}) = -\frac{C_4}{r_{\mathrm{ai}^4} + \frac{C_8}{r_{\mathrm{ai}^8}.$6

with force function

$V_{\mathrm{ai}(r_{\mathrm{ai}) = -\frac{C_4}{r_{\mathrm{ai}^4} + \frac{C_8}{r_{\mathrm{ai}^8}.$7

where

$V_{\mathrm{ai}(r_{\mathrm{ai}) = -\frac{C_4}{r_{\mathrm{ai}^4} + \frac{C_8}{r_{\mathrm{ai}^8}.$8

and, for the benchmark potential,

$V_{\mathrm{ai}(r_{\mathrm{ai}) = -\frac{C_4}{r_{\mathrm{ai}^4} + \frac{C_8}{r_{\mathrm{ai}^8}.$9

This is the representation integrated by ode45gpu (Chowdhury et al., 15 Sep 2025).

A notable modeling implication is that atomiongpu.m directly evolves the full trapped two-body dynamics rather than reducing the problem to relative coordinates. This follows from the nonseparability of center-of-mass and relative motion in the trapped case (Chowdhury et al., 15 Sep 2025).

3. Initial conditions, ensembles, and trap models

In the benchmark setup, the atom starts on a sphere of radius $^+$0 around the ion, with velocity pointing toward the origin: $^+$1 with

$^+$2

This benchmark corresponds to a head-on launch toward the ion from different incoming directions, with the ion initially at rest at the trap center (Chowdhury et al., 15 Sep 2025).

The general interface supports two collision ensembles. For collisiontype = "head-on", the atom is launched toward the origin at fixed speed

$^+$3

The ion starts at rest at the origin, and Tion is ignored. For collisiontype = "thermal", atom velocity components are sampled from independent normal distributions with standard deviation

$^+$4

ion velocity components are sampled from independent normal distributions with standard deviation

$^+$5

and ion coordinates are sampled from a normal distribution with standard deviation

$^+$6

where the secular frequency in $^+$7 is

$^+$8

Thus the code supports either a simplified directed-collision ensemble or a more realistic thermal one (Chowdhury et al., 15 Sep 2025).

Initial atom positions can be generated uniformly over a spherical angular region, with equal solid-angle weighting, or randomly over a spherical angular region. For uniform sampling, $^+$9 values are evenly spaced, while $$H= \frac{1}{2}m_av_a^2+\frac{1}{2}m_iv_i^2-\frac{C_4}{r_\textrm{ai}^4}+\frac{C_8}{r_\textrm{ai}^8} + \sum\limits_{i=1}^{3}(a_i+2q_i\cos\Omega_\textrm{rf}t)\frac{m_\textrm{ion}\Omega_\textrm{rf}^2}{8}r_i^2 .$$0 values are constructed using

$$H= \frac{1}{2}m_av_a^2+\frac{1}{2}m_iv_i^2-\frac{C_4}{r_\textrm{ai}^4}+\frac{C_8}{r_\textrm{ai}^8} + \sum\limits_{i=1}^{3}(a_i+2q_i\cos\Omega_\textrm{rf}t)\frac{m_\textrm{ion}\Omega_\textrm{rf}^2}{8}r_i^2 .$$1

with $$H= \frac{1}{2}m_av_a^2+\frac{1}{2}m_iv_i^2-\frac{C_4}{r_\textrm{ai}^4}+\frac{C_8}{r_\textrm{ai}^8} + \sum\limits_{i=1}^{3}(a_i+2q_i\cos\Omega_\textrm{rf}t)\frac{m_\textrm{ion}\Omega_\textrm{rf}^2}{8}r_i^2 .$$2 evenly spaced over a transformed interval, so each $$H= \frac{1}{2}m_av_a^2+\frac{1}{2}m_iv_i^2-\frac{C_4}{r_\textrm{ai}^4}+\frac{C_8}{r_\textrm{ai}^8} + \sum\limits_{i=1}^{3}(a_i+2q_i\cos\Omega_\textrm{rf}t)\frac{m_\textrm{ion}\Omega_\textrm{rf}^2}{8}r_i^2 .$$3 corresponds to equal area on the sphere (Chowdhury et al., 15 Sep 2025).

The trap representation supports three cases. For an untrapped ion, one sets

$$H= \frac{1}{2}m_av_a^2+\frac{1}{2}m_iv_i^2-\frac{C_4}{r_\textrm{ai}^4}+\frac{C_8}{r_\textrm{ai}^8} + \sum\limits_{i=1}^{3}(a_i+2q_i\cos\Omega_\textrm{rf}t)\frac{m_\textrm{ion}\Omega_\textrm{rf}^2}{8}r_i^2 .$$4

For a harmonic trap,

$$H= \frac{1}{2}m_av_a^2+\frac{1}{2}m_iv_i^2-\frac{C_4}{r_\textrm{ai}^4}+\frac{C_8}{r_\textrm{ai}^8} + \sum\limits_{i=1}^{3}(a_i+2q_i\cos\Omega_\textrm{rf}t)\frac{m_\textrm{ion}\Omega_\textrm{rf}^2}{8}r_i^2 .$$5

with $$H= \frac{1}{2}m_av_a^2+\frac{1}{2}m_iv_i^2-\frac{C_4}{r_\textrm{ai}^4}+\frac{C_8}{r_\textrm{ai}^8} + \sum\limits_{i=1}^{3}(a_i+2q_i\cos\Omega_\textrm{rf}t)\frac{m_\textrm{ion}\Omega_\textrm{rf}^2}{8}r_i^2 .$$6, $$H= \frac{1}{2}m_av_a^2+\frac{1}{2}m_iv_i^2-\frac{C_4}{r_\textrm{ai}^4}+\frac{C_8}{r_\textrm{ai}^8} + \sum\limits_{i=1}^{3}(a_i+2q_i\cos\Omega_\textrm{rf}t)\frac{m_\textrm{ion}\Omega_\textrm{rf}^2}{8}r_i^2 .$$7, and $$H= \frac{1}{2}m_av_a^2+\frac{1}{2}m_iv_i^2-\frac{C_4}{r_\textrm{ai}^4}+\frac{C_8}{r_\textrm{ai}^8} + \sum\limits_{i=1}^{3}(a_i+2q_i\cos\Omega_\textrm{rf}t)\frac{m_\textrm{ion}\Omega_\textrm{rf}^2}{8}r_i^2 .$$8, $$H= \frac{1}{2}m_av_a^2+\frac{1}{2}m_iv_i^2-\frac{C_4}{r_\textrm{ai}^4}+\frac{C_8}{r_\textrm{ai}^8} + \sum\limits_{i=1}^{3}(a_i+2q_i\cos\Omega_\textrm{rf}t)\frac{m_\textrm{ion}\Omega_\textrm{rf}^2}{8}r_i^2 .$$9, $m_a, v_a$0. For a Paul trap,

$m_a, v_a$1

and the user supplies $m_a, v_a$2, $m_a, v_a$3, and $m_a, v_a$4 (Chowdhury et al., 15 Sep 2025).

These choices make the script adaptable to different atom-ion systems, trap configurations, and collision ensembles without changing the underlying solver (Chowdhury et al., 15 Sep 2025).

4. Numerical implementation and GPU specialization

atomiongpu.m is built on ode45gpu, described as the authors’ custom MATLAB implementation of the same Dormand-Prince adaptive Runge-Kutta 4/5 algorithm used by MATLAB’s ode45. The paper states that the algorithm inside the solver is “exactly the same” as MATLAB ode45 in the sense of the RK stepping logic, but the implementation is specialized for this application (Chowdhury et al., 15 Sep 2025).

The main differences from MATLAB’s ode45 are structural. First, ode45gpu returns only final time, final state, and selected scalar observables computed on the fly, rather than the full time-dependent trajectory. This is crucial for running millions of trajectories with low memory (Chowdhury et al., 15 Sep 2025). Second, the implementation is written in a low-level scalar style because MATLAB arrayfun on GPU does not allow arbitrary vectorized MATLAB constructs. The 12 state variables are broken out explicitly, vector operations are manually unrolled, and the force evaluation is duplicated in the RK stages (Chowdhury et al., 15 Sep 2025). Third, even on CPU, a manually expanded scalar implementation avoids temporary vector construction and indexing overhead; the paper reports that rewriting the force function from a less efficient vectorized form to a scalar-expanded form gave about a 3× speedup in ode45 (Chowdhury et al., 15 Sep 2025).

The solver uses user-specified rtol and atol. It stops when simulated time reaches tmax, or step count reaches maxsteps. A safety cutoff exits if the simulation exceeds $m_a, v_a$5 timesteps (maxsteps in the interface) to prevent pathological runs from taking too long (Chowdhury et al., 15 Sep 2025). The paper states that the method is a fourth/fifth-order adaptive-timestep Runge-Kutta method; at each timestep the solver estimates the local error; if the error exceeds a threshold determined by the user’s relative and absolute tolerances, the timestep is reduced; otherwise the timestep may remain the same or increase (Chowdhury et al., 15 Sep 2025).

The GPU workflow is explicit: initial-condition arrays are created as matrices, converted to gpuArrays, constant parameters are expanded into same-shaped gpuArrays with repeated values, and arrayfun calls ode45gpu elementwise across those arrays. Each GPU thread effectively integrates one trajectory. Outputs are returned as gpuArrays, and gather transfers results back to CPU memory (Chowdhury et al., 15 Sep 2025). Because each atom-ion trajectory is independent and requires little memory if only final values are retained, this maps well to GPUs (Chowdhury et al., 15 Sep 2025).

This architecture places atomiongpu.m within a broader class of GPU-resident many-trajectory scientific workflows. A related but methodologically distinct example is a fully electromagnetic PIC workflow in PIConGPU, where GPU acceleration is combined with charge-conserving current deposition, Boris particle pushing, collisions, and atomic-kinetics sub-cycling for $m_a, v_a$6-atom helium nanodroplets (Medina et al., 12 Jun 2026). By contrast, atomiongpu.m retains a trajectory-by-trajectory ODE structure and classical binary-collision assumption, which is appropriate to its low-density atom-ion regime (Chowdhury et al., 15 Sep 2025).

5. Inputs, workflow, and outputs

The paper describes the end-to-end workflow as follows: the user edits top-level parameters in atomiongpu.m; the code interprets species, temperatures, trap, potential, and sampling setup; it constructs the full set of initial conditions for all trajectories; trajectories are distributed across CPU cores or GPUs depending on processor; each trajectory is integrated by ode45gpu; only final states and scalar observables are retained for memory efficiency; results are gathered to CPU memory; the script computes summary quantities such as complex-formation probability; it writes files including .csv raw observable arrays, .png heatmaps, trajectory plots, .txt summaries, and optional work.mat; and it can optionally rerun the longest-lived trajectory or a custom trajectory using CPU-based ode45gpu and MATLAB’s ode45 for validation and plotting (Chowdhury et al., 15 Sep 2025).

The script is designed to be adaptable to different atom-ion systems. Users can modify atom species via matom, ion species via mion, and the interaction potential model and parameters. Three potential input modes are described: potential = "CnCm", where the user supplies $m_a, v_a$7 and

$m_a, v_a$8

potential = "CnDe", where the user supplies $m_a, v_a$9, the code computes $m_i=m_{\rm ion}, v_i$0 analytically, and Re is ignored; and potential = "DeRe", where the user supplies $m_i=m_{\rm ion}, v_i$1, the code computes $m_i=m_{\rm ion}, v_i$2 and $m_i=m_{\rm ion}, v_i$3 analytically, and supplied Cn, Cm are ignored (Chowdhury et al., 15 Sep 2025).

The full input set includes species and masses (mion, matom), temperatures (Tion, Tatom), collision ensemble type (collisiontype = "head-on" or "thermal"), interaction potential model (potential = "CnCm", "CnDe", or "DeRe" together with n, m, Cn, Cm, De, Re), trap parameters (ax, ay, az, qx, qy, qz, OmegaRF), simulation control (tmax, maxsteps, r0, ntrajectories, rtol, atol), angular grid and initial-position controls (nanglesauto, ntheta, nphi, positions = "uniform" or "random", thetamin, thetamax, phimin, phimax), hardware settings (ncores, processor = "CPU" or "GPU"), full-trajectory plotting controls (longestlived, custom, t0custom, tfcustom, y0custom), and output/storage controls (saveworkspace, savecsvs, onlyonecsv) (Chowdhury et al., 15 Sep 2025).

The script creates a timestamped output directory containing .csv raw data files for observables, optional work.mat, probability.txt for complex-formation probability, outliers.txt for initial conditions of highest and lowest observable values, text printed to stdout or the MATLAB command window describing setup and elapsed time, heatmaps of observables versus initial $m_i=m_{\rm ion}, v_i$4, and trajectory plots for custom or longest-lived cases, including distance versus time and $m_i=m_{\rm ion}, v_i$5, $m_i=m_{\rm ion}, v_i$6, and $m_i=m_{\rm ion}, v_i$7 projections (Chowdhury et al., 15 Sep 2025).

6. Observables, validation, and limitations

The script computes and saves a range of scalar observables. The paper names the options for raw CSV output as "angle" for atom scattering angle, "bounces" for number of short-range collisions, "lifetime" for complex lifetime, "position" for ion displacement from origin during the first bounce, "transfer" for momentum transferred by the atom, "dist" for final atom-ion distance, "nsteps" for number of timesteps, and "KE" for final ion kinetic energy. All observables are recorded in atomic units (Chowdhury et al., 15 Sep 2025).

The code “automatically finds the probability of complex formation.” In the runtime printout example, it reports the percentage of trajectories which formed a complex. The text implies this is determined from trajectory analysis based on bounce behavior and lifetimes, though an explicit mathematical definition is not printed in the provided text (Chowdhury et al., 15 Sep 2025). The complex lifetime is explicitly defined as the amount of time between the first and last bounce, and the number of bounces is defined as the number of short-range collisions between atom and ion throughout the trajectory (Chowdhury et al., 15 Sep 2025).

For longest-lived trajectory mode, if longestlived = "yes", the script identifies the trajectory with the largest complex lifetime after the ensemble run, reruns that exact initial condition using ode45 on CPU, and produces full trajectory plots. For custom trajectory mode, if custom = "yes", the user supplies t0custom, tfcustom, and y0custom as a 12-element vector, and the code simulates and plots that specific trajectory (Chowdhury et al., 15 Sep 2025).

The authors verified that ode45gpu and MATLAB ode45 produce the same values for the atom’s scattering angle over 2,500 test trajectories (Chowdhury et al., 15 Sep 2025). At the same time, for long-lived, chaotic trajectories, rerunning exactly the same initial conditions on CPU and GPU may give different lifetimes or bounce counts. The paper attributes this to different floating-point roundoff behavior on CPU and GPU architectures and does not treat it as a bug, but as expected sensitivity in chaotic dynamics (Chowdhury et al., 15 Sep 2025).

Several limitations are explicit. atomiongpu.m describes classical dynamics only; there is no internal-state or quantum-scattering treatment. It assumes binary collisions only, valid for low-density atom clouds. The interaction potential is restricted to a generalized two-term inverse-power form, though that can be fitted to more realistic potentials. Full trajectories are not stored for ensemble runs, only final states and low-memory observables. The GPU implementation requires scalarized code and is less general than MATLAB’s high-level ODE interface. Chaotic trajectories can diverge numerically across architectures. Very large trajectory counts may hit memory limits; the authors observed issues around 25 million trajectories, and memory may become problematic above about 10 million (Chowdhury et al., 15 Sep 2025).

A common misconception would be to treat atomiongpu.m as a quantum-collision code or a many-body atom cloud simulator. The paper does not support either interpretation. It is instead a classical binary-collision simulator specialized to the low-density regime, with high-throughput ensemble statistics and optional selected-trajectory reruns (Chowdhury et al., 15 Sep 2025).

7. Performance, significance, and relation to adjacent GPU atom–ion methods

The paper reports several performance claims. On one CPU core, ode45gpu runs 22× faster than MATLAB ode45 for the atom-ion problem considered. The 22× speedup is maintained when parallelizing over up to 28 CPU cores. On GPUs, the code was benchmarked on NVIDIA Tesla K80 24GB, Tesla P100 16GB, and Tesla V100 32GB. Runtime stays roughly flat up to about 10,000 trajectories, showing effective parallel execution; beyond that threshold, runtime becomes roughly proportional to the number of trajectories, as GPU resources saturate. Using 8 GPUs, the code ran more than 10 million trajectories in under 15 hours (Chowdhury et al., 15 Sep 2025).

The software assumptions are MATLAB R2023a, Parallel Computing Toolbox, and GPU usage through MATLAB’s gpuArray and arrayfun. Hardware examples used in benchmarks include an Intel Haswell CPU node and NVIDIA Tesla K80, P100, and V100 GPUs (Chowdhury et al., 15 Sep 2025). There is also a caveat that future MATLAB versions may not support older GPUs as well; the authors used MATLAB 2023a (Chowdhury et al., 15 Sep 2025).

In technical context, atomiongpu.m occupies a distinct point in the computational landscape. It is substantially lighter-weight than production plasma frameworks such as PIConGPU, which replace explicit pairwise long-range interactions by self-consistent electromagnetic fields on a mesh and couple PIC evolution to Monte Carlo collisions and local atomic kinetics (Medina et al., 12 Jun 2026). It is also different from high-dimensional quantum collision solvers such as GPU-ported WP-CCC implementations, where the dominant cost arises from coupled-channel matrix assembly and linear solves along many impact-parameter trajectories (Abdurakhmanov et al., 2024). The central novelty of atomiongpu.m lies instead in a carefully specialized MATLAB realization of RK45 for massively parallel independent classical atom-ion trajectories (Chowdhury et al., 15 Sep 2025).

This suggests a broader methodological significance. atomiongpu.m shows that a MATLAB front end can remain viable for a narrowly defined but experimentally relevant atom-ion workflow if the solver is aggressively specialized for GPU-compatible scalar execution and low-memory observable extraction (Chowdhury et al., 15 Sep 2025). A plausible implication is that its main contribution is infrastructural: it lowers the practical barrier to ensemble-scale classical atom-ion collision studies without requiring migration to a full HPC codebase.

For cold and ultracold atom-ion research, the script’s most important role is therefore not to expand the underlying physical model, but to make large-ensemble statistical studies operationally straightforward within MATLAB, while preserving validation against MATLAB ode45 for benchmark observables such as the scattering angle (Chowdhury et al., 15 Sep 2025).