---
title: 'atomiongpu.m: GPU-Accelerated Atom-Ion Dynamics'
url: https://www.emergentmind.com/topics/atomiongpu-m
type: topic
---

# atomiongpu.m: GPU-Accelerated Atom-Ion Dynamics

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” [2509.12381]. 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 [2509.12381].

## 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 [2509.12381]. 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 [2509.12381].

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 [2509.12381].

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 [2509.12381].

## 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 \(\Omega_{\rm rf}\) is the trap drive frequency [2509.12381].

The dynamical variables are assembled into a 12-component state vector \(\mathbf y\), with \(y_{1,2,3}\) the ion position, \(y_{4,5,6}\) the ion velocity, \(y_{7,8,9}\) the atom position, and \(y_{10,11,12}\) the atom velocity. The ODE is written as
\[
\mathbf{\dot y}=f(t,\mathbf y),
\]
with force function
\[
\begin{aligned}
&f(t,\bold{y})=(y_4,y_5,y_6,\ &
             -\frac{1}{m_i}\frac{dV}{dr}\frac{y_1-y_7}{r}-(a_x+2q_x\cos\Omega_\textrm{rf}t)\frac{\Omega_\textrm{rf}^2}{4}y_1,\ &
             -\frac{1}{m_i}\frac{dV}{dr}\frac{y_2-y_8}{r}-(a_y+2q_y\cos\Omega_\textrm{rf}t)\frac{\Omega_\textrm{rf}^2}{4}y_2,\ &
             -\frac{1}{m_i}\frac{dV}{dr}\frac{y_3-y_9}{r}-(a_z+2q_z\cos\Omega_\textrm{rf}t)\frac{\Omega_\textrm{rf}^2}{4}y_3,\ &
             y_{10},y_{11},y_{12},\ &
             -\frac{1}{m_a}\frac{dV}{dr}\frac{y_7-y_1}{r},\ &
             -\frac{1}{m_a}\frac{dV}{dr}\frac{y_8-y_2}{r},\ &
             -\frac{1}{m_a}\frac{dV}{dr}\frac{y_9-y_3}{r}),
\end{aligned}
\]
where
\[
r=r_\textrm{ai}=\sqrt{(y_1-y_7)^2+(y_2-y_8)^2+(y_3-y_9)^2}
\]
and, for the benchmark potential,
\[
\frac{dV}{dr}=\frac{4C_4}{r^5}-\frac{8C_8}{r^9}.
\]
This is the representation integrated by `ode45gpu` [2509.12381].

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 [2509.12381].

## 3. Initial conditions, ensembles, and trap models

In the benchmark setup, the atom starts on a sphere of radius \(r_0=5000a_0\) around the ion, with velocity pointing toward the origin:
\[
\begin{aligned}
\bold{y_0}=\,&(0,0,0,0,0,0,\ &
r_0\sin\theta\cos\phi,r_0\sin\theta\sin\phi,r_0\cos\theta,\ &
-v_0\sin\theta\cos\phi,-v_0\sin\theta\sin\phi,-v_0\cos\theta),
\end{aligned}
\]
with
\[
v_0=\sqrt{\frac{3k_BT}{m_a}.
\]
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 [2509.12381].

The general interface supports two collision ensembles. For `collisiontype = "head-on"`, the atom is launched toward the origin at fixed speed
\[
v_0=\sqrt{\frac{3k_BT_{\rm atom}{m_{\rm atom}.
\]
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
\[
\sqrt{\frac{k_BT_{\rm atom}{m_{\rm atom},
\]
ion velocity components are sampled from independent normal distributions with standard deviation
\[
\sqrt{\frac{k_BT_{\rm ion}{m_{\rm ion},
\]
and ion coordinates are sampled from a normal distribution with standard deviation
\[
\sqrt{\frac{k_BT_{\rm ion}{m_{\rm ion}\omega^2},
\]
where the secular frequency in \(x\) is
\[
\omega=\frac{1}{2}\Omega_\textrm{rf}\sqrt{a_x+\frac{1}{2}q_x^2}.
\]
Thus the code supports either a simplified directed-collision ensemble or a more realistic thermal one [2509.12381].

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, \(\phi\) values are evenly spaced, while \(\theta\) values are constructed using
\[
\theta=\arccos(1-2u),
\]
with \(u\) evenly spaced over a transformed interval, so each \((\theta,\phi)\) corresponds to equal area on the sphere [2509.12381].

The trap representation supports three cases. For an untrapped ion, one sets
\[
a_x=a_y=a_z=q_x=q_y=q_z=\Omega_{\rm RF}=0.
\]
For a harmonic trap,
\[
V_\textrm{trap}(x,y,z)=\frac{1}{2}m_\textrm{ion}(\omega_x^2x^2+\omega_y^2y^2+\omega_z^2z^2),
\]
with \(q_x=q_y=q_z=0\), \(\Omega_{\rm RF}=2\), and \(a_x=\omega_x^2\), \(a_y=\omega_y^2\), \(a_z=\omega_z^2\). For a Paul trap,
\[
V_\textrm{trap}(x,y,z)=\sum_{i=x,y,z}(a_i+2q_i\cos\Omega_\textrm{rf}t)\frac{m_\textrm{ion}\Omega_\textrm{rf}^2}{8}r_i^2,
\]
and the user supplies \(a_i\), \(q_i\), and \(\Omega_{\rm rf}\) [2509.12381].

These choices make the script adaptable to different atom-ion systems, trap configurations, and collision ensembles without changing the underlying solver [2509.12381].

## 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 [2509.12381].

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 [2509.12381]. 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 [2509.12381]. 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` [2509.12381].

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 \(10^6\) timesteps (`maxsteps` in the interface) to prevent pathological runs from taking too long [2509.12381]. 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 [2509.12381].

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

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 \(10^6\)-atom helium nanodroplets [2606.14300]. 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 [2509.12381].

## 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 [2509.12381].

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 \(n,m,C_n,C_m\) and
\[
V(r)=-\frac{C_n}{r^n}+\frac{C_m}{r^m};
\]
`potential = "CnDe"`, where the user supplies \(n,m,C_n,D_e\), the code computes \(C_m\) analytically, and `Re` is ignored; and `potential = "DeRe"`, where the user supplies \(n,m,D_e,R_e\), the code computes \(C_n\) and \(C_m\) analytically, and supplied `Cn`, `Cm` are ignored [2509.12381].

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`) [2509.12381].

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 \((\theta,\phi)\), and trajectory plots for custom or longest-lived cases, including distance versus time and \(xy\), \(xz\), and \(yz\) projections [2509.12381].

## 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 [2509.12381].

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 [2509.12381]. 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 [2509.12381].

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 [2509.12381].

The authors verified that `ode45gpu` and MATLAB `ode45` produce the same values for the atom’s scattering angle over 2,500 test trajectories [2509.12381]. 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 [2509.12381].

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 [2509.12381].

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 [2509.12381].

## 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 [2509.12381].

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 [2509.12381]. There is also a caveat that future MATLAB versions may not support older GPUs as well; the authors used MATLAB 2023a [2509.12381].

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 [2606.14300]. 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 [2403.04252]. The central novelty of `atomiongpu.m` lies instead in a carefully specialized MATLAB realization of RK45 for massively parallel independent classical atom-ion trajectories [2509.12381].

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 [2509.12381]. 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 [2509.12381].

Source: https://www.emergentmind.com/topics/atomiongpu-m