Computer Physics Communications

Modelling has become a third distinct line of scientific enquiry, alongside experiments and theory. Molecular dynamics (MD) simulations serve to interpret, predict and guide experiments and to test and develop theories. A major limiting factor of MD simulations is system size and in particular the difficulty in handling, storing and processing trajectories of very large systems. This limitation has become significant as the need to simulate large system sizes of the order of billions of atoms and beyond has been steadily growing. Examples include interface phenomena, composite materials, biomaterials, melting, nucleation, atomic transport, adhesion, radiation damage and fracture. More generally, accessing new length and energy scales often brings qualitatively new science, but this has currently reached a bottleneck in MD simulations due to the traditional methods of storing and post-processing trajectory files. To address this challenge, we propose a new paradigm of running MD simulations: instead of storing and post-processing trajectory files, we calculate key system properties on-the-fly. Here, we discuss the implementation of this idea and on-the-fly calculation of key system properties in the general-purpose MD code, DL_POLY. We discuss code development, new capabilities and the calculation of these properties, including correlation functions, viscosity, thermal conductivity and elastic constants. We give examples of these on-the-fly calculations in very large systems. Our developments offer a new way to run MD simulations of large systems efficiently in the future.

We present PaScaL_TDMA 2.1, a GPU-oriented release of the PaScaL_TDMA library [3] for efficiently solving large batches of distributed tridiagonal systems on modern multi-GPU platforms. Building on the original CPU-based PaScaL_TDMA formulation and the shared-memory buffering strategy introduced in PaScaL_TDMA 2.0 [2], version 2.1 reformulates the core kernels and communication path to better match the GPU execution model. CUDA threads are mapped to contiguous tridiagonal lines to achieve coalesced global-memory access, and the elimination kernels are optimized to a fully register-resident implementation to reduce memory traffic and synchronization. To lower inter-GPU overhead, the reduced-system assembly is performed via a single consolidated MPI_Alltoall exchange, and the kernel interface is restructured to eliminate descriptor transfers at launch. Benchmarks on the NURION system show that PaScaL_TDMA 2.1 reduces wall time from 0.127 s on dual-socket Intel Skylake CPUs to 9.2 ms on an NVIDIA A100 and 6.1 ms on an H100, corresponding to speedups of 14.0 ×  and 20.7 × , respectively. Strong- and weak-scaling studies quantify the performance gains from the optimization stages and demonstrate sustained scalability on multi-GPU systems. Finally, PaScaL_TDMA 2.1 is integrated into an immersed-boundary LES solver and validated through large-scale CFD simulations, including an industrial-scale cleanroom configuration with up to 128 A100 GPUs and O(10^10) degrees of freedom.

The open-source Python package su4-branching automates the derivation of comprehensive spin S and isospin T branching rules for SU(4) irreducible representations. The Wigner supermultiplet scheme underlying nuclear and hadronic physics depends critically on SU(4) symmetry. However, practical calculations of branching rules for arbitrary SU(4) irreps have been largely inaccessible to the research community. Our implementation combines Racah’s tensor contraction formula with Python modularity, enabling straightforward exploration of high-dimensional SU(4) irreps through interactive Jupyter interfaces and multiple export formats. Comprehensive validation against three independent reference frameworks (Quesne 1976; Patera 1981; Pan et al. 2024) and dimensional consistency checks (necessary and sufficient for correctness) demonstrate reliability across dimensions 1–10^7. This work enables systematic group-theoretical investigations in nuclear structure, particle physics, and quantum chemistry.

We present a cross-architecture high-order heterogeneous Navier-Stokes simulation solver, XFluids, for compressible reacting multicomponent flows on different platforms, where ‘X’ stands for multiple different devices. The multicomponent reacting flows are ubiquitous in many scientific and engineering applications, while their numerical simulations are usually time-consuming to capture the underlying multiscale features. Although heterogeneous accelerated computing is significantly beneficial for large-scale simulations of these flows, the effective utilization of various heterogeneous accelerators with different architectures and programming models in the market remains a challenge. To address this, we develop XFluids by SYCL, to perform acceleration directly targeted to different devices, without translating any source code. This solver has been open-sourced, and tested on multiple graphics processing units (GPUs) from different mainstream vendors, indicating high portability. Through various benchmark cases, including the shock tube, diffusion, autoignition, detonation, and shock-bubble interaction, the accuracy of XFluids is demonstrated, with approximately no efficiency loss compared to popular existing GPU programming models, such as CUDA and HIP. In addition, XFluids show considerable acceleration compared to other open-source multicomponent reacting flow solvers. Then, to extend the solver to multiple GPUs, the Message Passing Interface (MPI) library is employed, with the GPU-aware communication supported. With this, the weak scaling of XFluids for multi-GPU devices is larger than 95% for 1024 GPUs. Last, in order to fully exploit the computational capability of the all devices, the hybrid CPU-GPU heterogeneous simulations are achieved without changing the source code of XFluids.

The detection of gravitational waves from binary black hole and neutron star mergers by ground-based interferometers, as well as the evidence for a gravitational wave background from pulsar timing array experiments, has marked a new era in astrophysics and cosmology. These experiments also have great potential for discovering new physics through gravitational wave detection. One of the most motivated sources of gravitational waves that can be realized only within a beyond-the-Standard-Model framework is first-order phase transitions. In this work we release PT2GWFinder, a Mathematica package designed to compute phase transition parameters and the gravitational wave power spectrum for an arbitrary scalar theory exhibiting a first-order phase transition, in scenarios where a single scalar acquires a vacuum expectation value. PT2GWFinder performs the phase tracing, computes the bounce profile and action using FindBounce, calculates the relevant temperatures and phase transition parameters, and finally evaluates the gravitational wave spectrum. Additionally, it offers a user-friendly interface with DRalgo, which enables the computation of the dimensionally reduced effective potential in the high-temperature regime. This work includes a user manual and two models that demonstrate the capability and performance of PT2GWFinder. As a supplement, for one of these models we obtain the bounce solution and action analytically in the thin-wall approximation and demonstrate excellent agreement with the numerical approach.

Differentiable programming is revolutionizing computational science by enabling automatic differentiation (AD) of numerical simulations. While first-order gradients are well-established, second-order derivatives (Hessians) for implicit functions in finite-element-based differentiable physics remain underexplored. This work bridges this gap by deriving and implementing a framework for implicit Hessian computation in PDE-constrained optimization problems. We leverage primitive AD tools (Jacobian-vector product/vector-Jacobian product) to build an algorithm for Hessian-vector products and validate the accuracy against finite difference approximations. Four benchmarks spanning linear/nonlinear, 2D/3D, and single/coupled-variable problems demonstrate the utility of second-order information. Results show that the Newton-CG method with exact Hessians accelerates convergence for nonlinear inverse problems (e.g., traction force identification, shape optimization), while the L-BFGS-B method suffices for linear cases. Our work provides a foundation for integrating second-order implicit differentiation into differentiable physics engines, enabling faster and more reliable PDE-constrained optimization.

Variational quantum algorithms (VQAs) have emerged as promising approaches for solving partial differential equations on near-term quantum hardware, offering hybrid optimization schemes that are robust to noise. A key challenge lies in balancing decomposition complexity and circuit depth: reducing the number of unitary terms lowers measurement cost but increases circuit depth, thereby amplifying hardware-induced noise. This work benchmarks three VQA-based solvers for the Poisson equation-Permutation Operator (PO), Sparse Decomposition (SD), and Fourier Diagonalization (FD) methods-all built on a common cost function. The FD method, in particular, leverages the Quantum Fourier Transform to diagonalize circulant Laplacians, enabling a decomposition with only one term per spatial dimension. The trade-off between gate count and noise resilience is systematically analyzed under realistic hardware constraints, including limited qubit connectivity and native gate sets. Empirical results on IBM quantum hardware indicate that while the FD method achieves superior cost scaling and requires fewer measurements, its deeper circuits are more sensitive to hardware noise. In contrast, the SD method provides better noise resilience with shallower circuits, but at the expense of a larger number of measurements. To facilitate reproducibility and further research, an open-source Qiskit-based library is provided for solving Poisson equations with mixed boundary conditions, supporting both IonQ simulators and IBM devices.

Monte Carlo (MC) simulations are a cornerstone of scientific computing in fields like medical physics, but their complexity often poses significant usability challenges. Setting up simulations requires intricate configuration and 3D geometry definition, which are error-prone and time-consuming tasks. Furthermore, the integration of MC tools into modern, Python-centric scientific workflows for analysis and AI can be difficult. This work addresses these challenges for the penRed MC code by introducing a comprehensive framework designed to enhance its accessibility, usability, and integration. We present pyPenred, a high-performance Python module that exposes the complete capabilities of penRed within the Python ecosystem. Built with pybind11, it allows computationally intensive particle transport to be handled by optimized C++ binaries while enabling seamless control and analysis in Python. Performance benchmarks show a minimal overhead of only 1–2% for locally compiled versions compared to native C++ execution. To simplify geometry creation and simulation setup, we developed a dedicated Blender plug-in. This integrated graphical environment supports constructing models with both quadric surfaces and triangular meshes, and provides an intuitive interface for defining materials, sources, and tallies. Finally, we have established robust cross-platform compatibility through continuous integration, automatically distributing pre-compiled binaries and pip-installable Python wheels for Linux, Windows, and macOS. Collectively, these contributions transform penRed from a specialized code-centric tool into an integrated and user-friendly simulation platform, lowering the barrier to advanced MC simulations and fostering tighter integration with contemporary data science workflows.

In this work, we report a new version of the HOS-Ocean code which is an open-source solver for deterministic nonlinear ocean wave propagation. The numerical model makes use of the so-called High-Order Spectral method, which ensures high efficiency and accuracy. This new release includes i) additional physical features such as spatially varying current and bathymetry together with dedicated models to account for wave-breaking, and ii) numerical developments with parallelization of the code through MPI as well as a more user-friendly code (pre-built binaries available, and simplified building procedure and operation). HOS-Ocean v2.1 is released as open-source, available from GitLab, developed and distributed under the terms of GNU General Public License (GPLv3). Along with the source code, detailed documentation under Sphinx format is available.

We present a novel platform for higher-order spectral analysis of time series data in Python. The theory, utility, and applications of such analyses are summarized. Direct estimation of coherence (n = 2), bicoherence (n = 3), and tricoherence (n = 4) spectra are given for test signals; higher-order (n > 4) spectra are inferred at single points in polyfrequency space. Quantification of uncertainty for nonstationary processes is considered, and applications to nonlinear dynamics research are given.