Computer Physics Communications

A new library VecAmpFit for multidimensional amplitude analyses in high-energy physics has been developed for an ongoing amplitude analysis at Belle II experiment. It includes a fitter performing likelihood calculation and explicitly-vectorized subprograms for amplitude implementation. The fitter supports explicit gradient calculation and simultaneous fitting of multiple data sets.

We present λPIC, a Python-based electromagnetic particle-in-cell framework built around a callback-centric architecture. Existing PIC codes typically tie high performance to static, pre-compiled timestep loops, hindering the implementation of custom physics, diagnostics, or output logic. λPIC breaks this coupling by exposing every stage of the loop as a named stage (hook), permitting attaching arbitrary Python functions that operate on the simulation itself, enabling custom algorithms and in-situ analysis without modifying the core algorithms. Under this flexible framework, performance-critical kernels are written in C extensions and Numba, fields and particles are stored in NumPy arrays, and MPI parallelism is paired with graph partitioning to support dynamic load balancing and non-rectangular domains. Although λPIC is designed to be general-purpose, it has special focus on intense laser-plasma interactions. Future work will extend the framework to GPU acceleration and additional physics modules including implicit solvers and nuclear physics.

The Monte Carlo wave-function (MCWF) method estimates observables by averaging over an ensemble of stochastic quantum trajectories. In standard practice, the ensemble size N is fixed before the simulation begins. We demonstrate that this strategy is statistically inadequate as a precision-control strategy: the statistical precision of the result is not a property of N alone, but of N jointly with the observable-specific variance-to-mean-squared ratio σ^2/μ^2, which is unknown prior to simulation. Concretely, on a single resonantly driven two-level atom, three observables require trajectory counts spanning a factor of 16.7 ×  to achieve the same 1% relative precision—a spread invisible to any fixed-N protocol. The consequence is systematic: fixed-N runs simultaneously under-sample high-variance observables (achieved estimated relative SEM up to 5.8 ×  the target, with no user warning) and over-sample low-variance ones (wasted compute factors up to 29 × ). The statistical root cause is the law N_opt=σ^2/(ε|μ|)^2: the required ensemble size is a function of per-observable population statistics that cannot be collapsed into a universal trajectory budget. The principled resolution is online inference: estimate σ^2 and |μ| from the trajectories themselves and stop as soon as the target precision is reached. We formalise this as an adaptive stopping criterion based on the relative standard error of the mean (RelSEM < ε), computed via Welford’s algorithm at O(1) overhead per trajectory. While tolerance-based Monte Carlo stopping is not itself new (Section 1), we show empirically that an existing jackknife-based implementation lacks a minimum-trajectory safeguard and can halt after as few as two trajectories with error exceeding 90%; LiQuID’s N_min-guarded Welford criterion avoids this failure mode. Three quantitative benchmarks on analytically tractable systems confirm that the criterion halts within  ∼ 13% of the asymptotic optimum and eliminates both failure modes of fixed-N simultaneously.

In molecular dissociation and fragmentation, some degrees of freedom that are vibrational at the saddle-point geometry correspond asymptotically to rotations or translations of the product fragments. Their correct classification as vibrational or roto-translational is essential in kinetic and statistical models that treat these degrees of freedom separately, but current assignments often rely on visual inspection or system-specific criteria. Here, we present Portico (Projection Onto Roto-Translational and Internal COordinates), a program for the automated classification of transition-state normal modes according to their asymptotic product character. First, Portico assigns the imaginary-frequency mode to the roto-translational sector by construction because it defines the reaction coordinate leading toward product separation. Second, Portico represents the transition-state normal modes in a mixed representation that combines internal coordinates selected to span the vibrational degrees of freedom of the isolated product fragments with rigid-body translational and rotational displacement vectors of those fragments. The components of each mode in this representation define a scalar weight, Ωi, which quantifies its roto-translational character. Since the total number of transition-state normal modes to be assigned to the roto-translational sector is fixed by the dimensionality of the dissociation channel, Portico completes this sector by ranking the real-frequency modes according to their Ωi values and selecting the required number of modes. The same analysis can also be applied at non-stationary geometries, where the gradient direction is projected out and plays the role of the reaction coordinate. The method applies to channels yielding atomic, linear, and nonlinear fragments and does not require propagation toward the asymptotic product region. Portico thus provides an objective and transferable tool for automated reaction analysis and statistical models of molecular fragmentation.

We present the third release of ULYSSES, a Python package for the numerical evaluation of the baryon asymmetry generated through leptogenesis. This version includes code implementing state-of-the-art density matrix equations for low-scale leptogenesis with three quasi-degenerate right-handed neutrinos. We extend the validity of the code in this scenario beyond the 100 GeV right-handed neutrino mass scale, into the regime of resonant leptogenesis, by including neutrino production rates valid in both the relativistic and non-relativistic regimes. In addition, in the high-scale vanilla scenario, we provide routines for computing ΔL=1 scattering processes, enabling full phase-space evolution of the right-handed neutrino and lepton asymmetry. A new --extended parameter interface allows users to pass model-specific inputs beyond the standard leptogenesis runcard without modifying the core infrastructure and demonstrate its use with a toy module that simultaneously solves the vanilla leptogenesis equations and the freeze-in production of dark matter. On top of these improvements, we introduce an alternative parametrisation of the Casas-Ibarra matrix, update the default neutrino oscillation parameters and report cross-checks of the new low-scale leptogenesis module against published benchmarks and independent codes. ULYSSES is publicly available on GitHub https://github.com/earlyuniverse/ulysses and pip-installable from PyPI.

The self-organised morphogenesis of multicellular organisms emerges not only from biochemical but also from mechanical processes operating across spatial scales. Yet, computational tools for simulating tissue-scale mechanics have largely been released as restricted, study-specific codes. Here, we present NewtonBioMorphX, a GPU-accelerated, open-access comprehensive framework for simulating tissue-scale mechanical morphogenesis. Built on a finite-strain finite-element formulation on tetrahedral meshes, the software implements an unusually broad set of mechanical ingredients: isotropic and anisotropic (direction-specific) growth; multi-layer tissue mechanics; neo-Hookean hyperelasticity combined with stretch-based plasticity; Kelvin–Voigt viscoelasticity; passive anisotropic fibres emulating 3D collagen networks; anisotropic active contractile filaments emulating actin; Hertzian contact forces; and multiple orthonormal coordinate bases (Cartesian, normal–tangent, cylindrical, spherical, cone-adapted, and toroidal). Equally central is an extensive, user-friendly graphical interface spanning the entire workflow — mesh generation from geometric primitives or arbitrary mesh input, per-layer material and growth specification (including analytical functions of space and time compiled at runtime), boundary-condition assignment, real-time colour-mapped and glyph-based visualisation of scalar and tensor fields, and structured data export — allowing intricate mechanical modelling of morphogenetic processes without editing or recompiling source code. The NewtonBioMorphX framework is validated using a benchmark implementing the well-known Cook’s membrane problem, and its versatility is illustrated through parameter gradient plots, as well as twelve ready-to-run preset examples spanning diverse morphologies.

We introduce Merlino2D, an efficient fluid plasma solver (available at https://github.com/PTL-Unibo/Merlino2D) for the simulation of low-temperature gas discharges and plasma-based devices. The solver operates on structured and unstructured triangular meshes and supports robust two-dimensional simulations of stationary and time-dependent regimes. Merlino2D integrates the open-source LisbOn KInetics Boltzmann solver (LoKI-B), which computes the electron energy distribution function and associated swarm parameters. The time-dependent drift–diffusion equations coupled with Poisson’s equation are formulated as a monolithic system of differential-algebraic equations (DAEs). The system of equations is integrated in time with a fully implicit time integration scheme with adaptive time stepping and error control. This approach results in efficient handling of numerical stiffness and favorable scaling of the computational time with mesh size. In this work, the numerical methodology underlying Merlino2D is described, and the numerical accuracy of the transport schemes is verified against analytical solutions. In addition, the solver is applied to two-dimensional simulations of both DC and AC discharges, including corona discharges and surface dielectric barrier discharges (DBDs), demonstrating its capability to reproduce experimental trends and results reported in the literature.

We present EEXICC, a Monte Carlo event generator designed to simulate the production of doubly heavy baryons (Ξ_{cc}, Ξ_{bc}, and Ξ_{bb}) via e+e− annihilation. Based on nonrelativistic QCD effective theory, the generator calculates the process e+ + e− → Ξ_{QQ′} + \bar{Q}′+\bar{Q} using an improved trace technique at the amplitude level, which greatly improves numerical efficiency compared with traditional squared-amplitude methods. EEXICC uses a modular Fortran implementation for hard-process generation and a C++ interface to the bundled PYTHIA 8.317 library, enabling convenient integration into complete event simulation workflows. The program supports both weighted and unweighted event generation, and its numerical reliability has been verified against existing theoretical results. EEXICC provides a flexible and robust tool for studying the properties of doubly heavy baryons at future high-luminosity and high-energy e+e− colliders such as the CEPC and FCC-ee.

High-fidelity electromagnetic scattering computation for nonspherical and internally inhomogeneous atmospheric particles is a key link between particle microphysics and observables in measurement modeling and remote-sensing forward simulation. In practice, however, the use of high-fidelity solvers is often hindered by limited openness, insufficient modularity for large-scale parameter sweeps, high computational cost for multi-orientation reuse, and weak integration with measurement-oriented workflows. We present TransitionMatrices.jl, an open-source Julia package centered on the Invariant Imbedding T-Matrix (IITM) method. The framework provides a modular pipeline covering particle modeling, IITM-based T-matrix computation, and region-wise heterogeneity partitioning acceleration: the inscribed-sphere region is treated as a homogeneous core domain and accelerated using a logarithmic-derivative-based Mie scattering method, while the region between the inscribed and circumscribed spheres is treated as an extended inhomogeneous domain with particle-shape-dependent acceleration strategies. The package further supports validation-oriented Extended Boundary Condition Method (EBCM) cross-checking, far-field scattering products (cross sections, amplitude matrix, scattering/Mueller matrix), and coordinate transformations for efficient orientation handling. We report systematic credibility checks, including Mie consistency in spherical cases, IITM–EBCM cross-validation where applicable, and rotation-related error assessment. Finally, we demonstrate WRF-coupled polarimetric forward simulations to illustrate how particle-shape assumptions propagate into simulated atmospheric observables.

qcombo is a Python package for the symbolic evaluation of commutators between general quantum many-body operators expressed in normal-ordered form using the generalized Wick theorem. The package provides an automated and systematic framework for generating the corresponding algebraic expressions, significantly reducing the risk of human error in lengthy and complex analytical derivations. It is designed to assist the development and implementation of modern many-body methods in nuclear physics, quantum chemistry, and related fields. The functionality and workflow of the package are demonstrated through an application to the in-medium similarity renormalization group (IMSRG) method, which has been widely used for nuclear ab initio calculations. As a representative example, qcombo is employed to automatically generate the complete set of multi-reference IMSRG flow equations with operators truncated at the normal-ordered three-body level.