GMTHRASHpy: Python for Beams & Monte Carlo

• GMTHRASHpy is a Python-based tool with dual applications in forward-convolution fitting for molecular beams experiments and autocorrelation analysis for Monte Carlo simulations.
• Its reaction dynamics module simulates laboratory observables from guessed center-of-mass reaction models, addressing reproducibility issues tied to legacy executables.
• The Γ-method implementation delivers rigorous error estimation and autocorrelation time computation through modular workflows and customizable parameter settings.

Searching arXiv for the cited GMTHRASHpy paper and closely related work to ground the article. GMTHRASHpy is a name used in arXiv sources for two distinct Python research codes. In its primary and most specific usage, it denotes a Python-based, open-source reimplementation of the GMTHRASH forward-convolution fitting tool for crossed molecular beams experiments, created to reproduce and extend the analysis of measured time-of-flight spectra and angular intensity distributions by mapping guessed center-of-mass reaction models into the laboratory frame (Fujioka et al., 6 Oct 2025). In a separate usage, the same name is attached to a modular Python implementation of the $\Gamma$-method for Monte Carlo simulations, intended for estimating autocorrelation times, corrected uncertainties, and error propagation for primary and derived observables in correlated Markov-chain output (Palma et al., 2017). The shared label therefore spans two methodological contexts: reaction-dynamics data inversion and statistical analysis of autocorrelated simulation data.

1. Nomenclature and disciplinary scope

In these sources, the label refers to two different programs rather than to a single universally fixed software package.

Usage of the name	Scientific setting	Core function
"GMTHRASHpy: Forward Convolutions of Crossed Molecular Beams Experiments in Python" (Fujioka et al., 6 Oct 2025)	Crossed molecular beams experiments	Forward convolution fits of laboratory-frame observables
"A Python program for the implementation of the $\Gamma$-method for Monte Carlo simulations" (Palma et al., 2017)	Monte Carlo simulations	Estimation of autocorrelation times and corrected errors

The first usage is tightly coupled to crossed-beam reaction dynamics and to the legacy GMTHRASH executable. The second usage belongs to statistical error analysis for Monte Carlo data and is centered on Wolff’s $\Gamma$-method. This dual usage is a source of potential ambiguity in citation and software identification, and a precise reference to the associated arXiv record is therefore important.

2. Crossed molecular beams problem addressed by the forward-convolution code

In crossed molecular beams experiments, the directly recorded observables are time-of-flight (TOF) spectra and angular intensity distributions of reaction products at various detector angles (Fujioka et al., 6 Oct 2025). The physically interesting quantities, however, are the product center-of-mass (CM) translational energy distribution and CM scattering angle distribution, and these are not directly measured. The required procedure is forward convolution: one guesses a CM model, propagates it through the beam geometry and detector response, and compares the resulting laboratory-frame prediction with experiment.

The crossed-beams GMTHRASHpy was developed specifically for this fitting problem. Its stated purpose is to reproduce and extend the analysis of such experiments by simulating how a guessed CM reaction model would appear in the laboratory frame, where the measurements are made. The source identifies this capability as central to studies of reaction dynamics in combustion, astrochemistry, and molecular mass growth.

A central motivation for the software is the status of the original GMTHRASH code base. The original source code is described as missing, with only an executable remaining. The source then identifies two resulting limitations: limited reproducibility and flexibility, because users cannot inspect or modify the algorithm easily, and limited capability, because the executable restricts certain behaviors such as a fixed $90^\circ$ crossing angle. GMTHRASHpy was created to fill this gap with a modern, editable code base that can replicate published GMTHRASH fits while also allowing new workflows, new physics models, and different experimental geometries.

3. Forward-convolution formalism and algorithmic workflow

The crossed-beams implementation follows a standard forward-convolution workflow (Fujioka et al., 6 Oct 2025). It begins by defining the beam properties: primary and secondary beam velocities, speed ratios, angular divergences, crossing angle, and beam masses. It then samples different primary and secondary beam angles to represent finite beam divergence and collision-volume width. For each collision condition it constructs Newton diagrams, computes the set of possible product laboratory velocities, and bounds these by energy conservation in the familiar Newton-circle geometry of velocity space.

The CM product distributions are then projected into the laboratory frame. The program samples the product CM translational energy distribution and scattering-angle distribution, and integrates these over detector angle, detector aperture, ionizer length, and flight path to generate predicted TOF and angular intensity profiles. The fit is obtained by comparing these predicted laboratory-frame signals with the measured data.

The paper gives the translational energy distribution in the form

$P(E_T)=N_{ET}(E_T-E_{T,\min})^p(E_{T,\max}-E_T)^q,$

where $p$ and $q$ shape the distribution and $N_{ET}$ normalizes it. The scattering-angle distribution is described as a Legendre-polynomial expansion with coefficients $C_i$ and normalization $N_{\theta_{CM}}$. These functions are the primary fit inputs used to represent the reaction dynamics.

The physical mapping from CM variables to detector observables is written through several relations. The source gives the energy-conservation limit on product motion in the CM frame as

$\Gamma$0

the laboratory/CM velocity relation as

$\Gamma$1

the detector-count relation as

$\Gamma$2

and the TOF conversion as

$\Gamma$3

The paper further writes

$\Gamma$4

and summarizes the forward-convolution integral in simplified form as

$\Gamma$5

A notable technical point is the treatment of the solid-angle Jacobian. The source emphasizes that the simpler inverse-square factor is not always sufficient and gives

$\Gamma$6

together with an additional correction

$\Gamma$7

This discussion is tied to earlier treatments by Morse, Bernstein, Zare, Suits, and Helbing. The paper also notes that if the detector angle and CM velocity are exactly aligned, the integral can become problematic under some formulations, which is part of why the exact Jacobian matters.

The contribution of each Newton diagram is weighted by physical likelihoods,

$\Gamma$8

with

$\Gamma$9

and with additional mass/time factors when different product channels are involved. This weighting is part of the attempt to preserve the physical likelihood structure of the sampled collision conditions.

4. Implementation, access, and reproducibility in the crossed-beams program

The crossed-beams GMTHRASHpy is described as Python-based, open-source, cross-platform, and installable from GitHub at github.com/kaka-zuumi/GMTHRASH, with pip installation under the package name gmthrash (Fujioka et al., 6 Oct 2025). It can be run from the command line, imported into Python, and accessed through a GUI via tkinter. The paper states that, like the original program, it can read inputs from a single compact PAN file; however, users may also modify inputs directly in Python.

These implementation choices address the reproducibility problem created by the missing original source. The program is explicitly framed as an editable replacement for an opaque executable, and its Python realization is presented as a practical means to inspect and modify the algorithm. A plausible implication is that the software shifts the forward-convolution workflow from executable-level black-box usage toward scriptable and inspectable analysis.

The source also stresses greater flexibility than the original executable. Because users are not constrained by fixed executable behavior or file-only control, the code can be adapted to new workflows, new physics models, and different experimental geometries. The explicit mention of the original executable’s fixed $\Gamma$0 crossing angle places this flexibility in direct contrast with the limitations of the legacy implementation.

5. Benchmarking against the original GMTHRASH executable

The crossed-beams code is explicitly designed to match the behavior of the original GMTHRASH executable, and the paper reports benchmarking across a wide variety of published reactions (Fujioka et al., 6 Oct 2025). The outputs are described as nearly identical in most cases. For single-channel fits, agreement is reported to be very close across test reactions. For multichannel fits, agreement is also generally close, although the paper notes small discrepancies in some cases, including certain C2 + isoprene fits.

The comparison is limited by the fact that the original program’s internal variables are hidden. As a result, benchmarking is performed only at the level of final output distributions rather than by direct internal-state comparison. This does not invalidate the comparison, but it constrains the kind of equivalence that can be demonstrated. The benchmark therefore establishes output-level reproducibility rather than algorithm-internal identity.

This output-level replication is important because the stated aim of the software is to reproduce published crossed-beam analyses. The source presents this as the basis for using GMTHRASHpy to reproduce data or fits from existing experiments while also extending the modeling framework. A plausible implication is that the code functions both as a reproduction platform for historical results and as a modifiable basis for future crossed-beam analyses.

6. GMTHRASHpy as a $\Gamma$1-method program for Monte Carlo simulations

A separate 2017 source uses the same name for a modular Python program devoted to the estimation of autocorrelation times for Monte Carlo simulations by means of the $\Gamma$2-method algorithm (Palma et al., 2017). In that usage, the problem is not reaction dynamics but the analysis of Monte Carlo output in the presence of autocorrelations. The program’s purpose is to estimate errors more reliably than naive or binning-based methods and to compute integrated autocorrelation times for both primary observables and derived quantities.

For a primary observable $\Gamma$3, measured as $\Gamma$4 on a Markov chain, the correlation matrix is given as

$\Gamma$5

with diagonal elements

$\Gamma$6

identified as autocorrelation functions. The normalized autocorrelation function is

$\Gamma$7

and the integrated autocorrelation time is defined by

$\Gamma$8

with

$\Gamma$9

The paper emphasizes the key consequence that, compared with the uncorrelated case, the variance of the mean is effectively multiplied by $90^\circ$0.

Because finite data make the tail of the correlation function noisy, the method uses a truncation window $90^\circ$1 rather than summing to infinity. The source describes this as a bias–variance tradeoff: truncation introduces a systematic bias,

$90^\circ$2

but summing too far into the noisy tail prevents proper variance reduction. The paper follows Wolff’s prescription and uses an automatic criterion based on an estimated decay time $90^\circ$3 and a tuning factor $90^\circ$4. The optimal window $90^\circ$5 is identified through the sign change of

$90^\circ$6

If no sign change occurs up to $90^\circ$7, the window is set to $90^\circ$8, and the final autocorrelation sum is computed up to $90^\circ$9.

The program also treats derived quantities $P(E_T)=N_{ET}(E_T-E_{T,\min})^p(E_{T,\max}-E_T)^q,$0 through a Taylor expansion around the exact mean values. The effective variance is written as

$P(E_T)=N_{ET}(E_T-E_{T,\min})^p(E_{T,\max}-E_T)^q,$1

and the integrated autocorrelation time for the derived quantity is

$P(E_T)=N_{ET}(E_T-E_{T,\min})^p(E_{T,\max}-E_T)^q,$2

The paper states that practical estimators are provided through replica and finite-difference derivatives, allowing the treatment of nonlinear functions without requiring explicit analytic error propagation. Support for $P(E_T)=N_{ET}(E_T-E_{T,\min})^p(E_{T,\max}-E_T)^q,$3 statistically independent replica is also included, with replica means $P(E_T)=N_{ET}(E_T-E_{T,\min})^p(E_{T,\max}-E_T)^q,$4 defined for separate files or subdivided long runs. The source notes that replica lengths must still be large compared to the autocorrelation time.

7. Architecture, workflows, and validation of the $P(E_T)=N_{ET}(E_T-E_{T,\min})^p(E_{T,\max}-E_T)^q,$5-method program

In the Monte Carlo usage, the software is packaged as a Python library and command-line tool and is described as modular and extensible (Palma et al., 2017). The package organization includes the modules ioutils, analysis, plots, configuration, session, unew, and supporting components. The main analysis classes are PrimaryAnalysis and DerivedAnalysis, which generate an AnalysisData object containing fields including value, rep_value, rep_mean, deviation, w_opt, t_max, dvalue, ddvalue, variance, naive_err, tau_int, dtau_int, rho, drho, and qval.

Graphical output is produced with matplotlib, including integrated autocorrelation time versus window $P(E_T)=N_{ET}(E_T-E_{T,\min})^p(E_{T,\max}-E_T)^q,$6, normalized autocorrelation function, histogram of replica values, and distribution of the underlying data. The paper lists the required Python packages as numpy, scipy, matplotlib, docopt, voluptuous, PyYAML, tqdm, and colorama. Installation is described through pip install ., python setup.py install, and pip install --user ., and the package is stated to be distributed under the MIT license.

Two principal input workflows are supported: direct file input from the command line and directory-based input in which the program scans a directory for data files. The software also supports arbitrary numbers of files, splitting files into multiple replica, command-line parameters, and YAML configuration. A YAML configuration can specify the data directory or replica files, file patterns, indices or ranges, the number of replica $P(E_T)=N_{ET}(E_T-E_{T,\min})^p(E_{T,\max}-E_T)^q,$7, the $P(E_T)=N_{ET}(E_T-E_{T,\min})^p(E_{T,\max}-E_T)^q,$8 parameter for windowing, which primary observables to analyze, parameters for derived functions, the Python module name, and target functions. This organization is presented as a practical basis for scripting, reproducing, and embedding analyses into larger workflows.

The paper’s operational test uses the 2D Ising model at criticality on square lattices and compares Metropolis with Wolff single-cluster updates. The analyzed quantities are the energy density $P(E_T)=N_{ET}(E_T-E_{T,\min})^p(E_{T,\max}-E_T)^q,$9 and the derived scaling quantity

$p$0

where $p$1 is the susceptibility. The source reports that the implementation reproduces the results of the older MATLAB implementation UWERR exactly. It further reports that the mean values from the two algorithms agree at the $p$2 level, that the Wolff algorithm exhibits much shorter autocorrelation times than Metropolis, that integrated autocorrelation time versus window $p$3 shows a plateau near $p$4, that the Metropolis dynamic critical exponent is near $p$5, and that, for the cluster algorithm, a rescaled autocorrelation time is needed for comparison with prior literature. The paper contrasts the $p$6-method with simpler binning techniques, arguing that the $p$7-method exploits correlated data more effectively and yields more reliable error estimates.