Computer Physics Communications

The MCPlas toolbox represents a collection of MATLAB functions for the automated generation of an equation-based fluid-Poisson model for non-thermal plasmas in the multiphysics simulation software COMSOL. Following the development of the new generation of the LXCat platform, all input data are prepared in a structured and interoperable JSON format and can be supplied and validated using existing JSON schemas. The toolbox includes fully transparent, editable MATLAB source code and offers an advanced description of electron transport in addition to commonly used approaches in the plasma modelling community. It supports one-dimensional and two-dimensional modelling geometries employing Cartesian, polar and cylindrical coordinate systems. MCPlas is tested on two reference cases: DC- and RF-driven low-pressure glow discharges in argon. Comparison of MCPlas results with results obtained by employing COMSOL’s Plasma Module verifies the reliability of the plasma model implemented by MCPlas and demonstrates the significance of electron transport treatment and boundary conditions applied in the toolbox. Using the same examples, the easy handling of complex reaction kinetic models in MCPlas and the reusability of its JSON input data across different modelling platforms are illustrated. This demonstrates that MCPlas provides a transparent and reproducible workflow for the simulation of non-thermal plasmas using COMSOL.

We present C-BerryTrans, a C++ code designed for first-principles calculations of Berry-curvature-driven transverse transport properties, namely the anomalous Hall conductivity (AHC) i.e., σ_{μν}^{AHC} and anomalous Nernst conductivity (ANC) i.e., α_{μν}^{ANC}. The code directly extracts eigenvalues and momentum-matrix elements from WIEN2k calculations and evaluates the Berry curvature (Ω) using a Kubo-like formalism, thereby avoiding the interpolation errors inherent in Wannier-based approaches. To ensure computational efficiency, C-BerryTrans parallelizes Ω evaluation over k-points using OpenMP and stores band-resolved curvature data in binary format, significantly reducing memory usage. This design enables rapid post-processing of AHC and ANC over a wide range of temperature (T) and chemical potential (ω) values in a single run. The code has been benchmarked on well-studied ferromagnetic materials- Fe, Fe3Ge, Pd, Fe3Al, and Co2FeAl. For Fe, the σ_{xy}^{AHC} is obtained to be  ∼ 775 ( ∼ 744) S/cm at 0 (300) K. In case of Fe3Ge, the calculated value of σ_{xy}^{AHC} is found to be 311 S/cm at 300 K. Nextly, for Co2FeAl, the magnitude of computed value of σ_{xy}^{AHC} at 2 K is found to be  ∼ 56 S/cm. Moving further, the room temperature magnitude of α_{xy}^{ANC} for Pd is obtained to be  ∼ 0.97 AK^{-1}m^{-1}. In case of Fe3Al, the maximum magnitude of α_{xy}^{ANC} for T ≤ 500 K is computed as  ∼ 2.83 AK^{-1}m^{-1}. Lastly, for Co2FeAl, the value of α_{xy}^{ANC} is obtained to be  ∼ 0.10 AK^{-1}m^{-1} at 300 K. These results show excellent agreement with previously reported data. With its accuracy, scalability, and user-friendly workflow, C-BerryTrans provides a powerful tool for exploring Ω-driven transport phenomena and is well suited for high-throughput materials discovery. The code further enables the evaluation of Ω-derived AHC/ANC contributions along user-defined high-symmetry k-point paths. This provides valuable microscopic insight into how specific band-structure features contribute to Ω-driven AHC/ANC. Additionally, the code is equipped with a visualization module that allows analysis of k-point contributions to AHC or ANC in any material. This further enhances its capability for exploring topological materials.

We present a Monte Carlo event generator specifically developed for the study of hadronic two-photon scattering events in two-photon scattering at electron-positron colliders. The code enables the generation of events with exact leading-order QED coupling and a flat phase space decay of the hadronic state into an arbitrary number of final state particles as selected by the user. Thus, this generator is well-suited for the use of partial wave analyses tools to study the two-photon production of higher-multiplicity final states across a wide range of energies and photon virtualities. Furthermore, the code integrates both experimental and theoretical inputs on the two-photon couplings of hadrons to simulate two-photon production processes. Motivated by the investigations of the BESIII collaboration, the final states π^{+}π^{-}, π^{0}π^{0}, π^{0}η, K^{+}K^{-}, K^{0}_{S}K^{0}_{S}, ηη, and f_1(1285) → ηπ^{+}π^{-} via a^{±}_{0}(980)π^{∓} and f_{0}(500)η are currently included. The code is sufficiently flexible to easily add additional final states as well as quickly change the already included channels.

Studying the interactions between energetic particles (EPs) and Alfvén eigenmodes (AEs) is essential for understanding alpha particle transport in burning plasma conditions within future fusion reactors. To numerically simulate EP-driven Alfvén instabilities, we have developed GMEC[1,2] (Gyrokinetic-Magnetohydrodynamics Energetic-particle Code) on the central processing unit (CPU) platform. However, long-time and large-scale simulations require higher grid resolution and more particles, demanding substantial computational resources. Therefore, applying more efficient numerical algorithms and optimization strategies on advanced computing architectures becomes critical. Computer scientists have designed the graphics processing unit (GPU) to implement single instruction multiple thread (SIMT) with fewer control units and higher throughput than CPUs. The GPU is optimized for data-parallel tasks and has emerged as a powerful platform for scientific computing. Based on NVIDIA’s Compute Unified Device Architecture (CUDA), this work presents cuGMEC, a high-performance gyrokinetic-MHD hybrid code developed from scratch for GPU platforms, incorporating newly added equations and physical terms. After testing on NVIDIA GPUs, the results show favorable scaling and acceleration, achieving up to 20 times speedup compared to the Intel Xeon Gold 6348 processor. Several benchmarks with theories and other codes have also been conducted to verify the gyrokinetic-MHD hybrid code.

Molecular dynamics is widely used to study various phenomena, such as diffusion, shock wave propagation, and plasma dynamics. A wide range of software packages supports the expanding scope of molecular dynamics applications. However, the quality of simulations depends on force field approximations, ranging from simple models to direct quantum solutions. Recently, machine learning approaches for constructing accurate interatomic potentials have received significant attention. In MDcraft, we integrate these advances into a scalable, physically accurate framework. MDcraft is a comprehensive, modern molecular dynamics platform. It offers a high-level Python API (Application Programming Interface) with a user-friendly, script-based interface. The core simulation algorithms are implemented in C++ to ensure robustness and computational efficiency. MDcraft is built for high-performance computing on modern clusters and supports dynamic domain decomposition and load balancing via the Message Passing Interface (MPI) for scalable parallelization. Additionally, MDcraft leverages multithreading within nodes through standard C++ parallelism, enabling efficient use of heterogeneous architectures. We demonstrate the code's capabilities through several examples, including the shock response in aluminum, the shock Hugoniot in argon, and the cold curve of copper.

Optical Bloch equations and rate equations serve as powerful tools to model light-matter interactions from textbook-like two-level atoms to the complex internal dynamics of molecules. A particular challenge in this context is posed by molecular laser cooling, where many dozens or hundreds of levels need to be taken into account for a comprehensive modeling. Here, we present MoleCool, a numerically efficient Python toolbox to implement and solve the corresponding differential equation systems. We illustrate both the capabilities of the toolbox and some of the intricacies of molecular laser cooling by educational examples, which range from simple Rabi oscillations to spontaneous and coherent cooling schemes for various currently studied or considered molecular species. This includes, in particular, a comprehensive modeling of laser cooling dynamics with full hyperfine structure resolution in radioactive radium monofluoride (RaF), as well as studies of other complex species such as barium monofluoride (BaF) and ytterbium monohydroxide (YbOH).

Density functional theory (DFT) remains the most widely used electronic structure method. Although exact in principle, in practice, it relies on approximations to the exchange-correlation (XC) functional, which is known to be a unique functional of the electron density. Despite 50 years of active research, existing XC approximations remain far from general purpose chemical accuracy of various thermochemical and materials properties. In that light, the inverse DFT problem, of finding the exact XC potential corresponding to an accurate groundstate density, offers an insightful tool to understand the nature of the XC functional as well as aid in the development of more accurate functionals. However, solving the inverse DFT problem is fraught with several numerical challenges, such as non-uniqueness or spurious oscillations in the solution and non-convergence. We present invDFT as an open-source framework to address the outstanding challenges in inverse DFT and computed XC potentials solely from a target density. We do so by use of a systematically convergent differential finite-element basis—higher-order finite-elements for the Kohn-Sham orbitals and linear finite-elements for the XC potential—which together render the inverse DFT problem well-posed. We also employ necessary asymptotic corrections to the target density to avoid any unphysical oscillations in the resulting XC potential. We also employ several numerical and high-performance computing (HPC) advances that affords both efficiency and parallel scalability, on CPU-GPU hybrid architectures. We demonstrate the accuracy and scalability of invDFT using accurate full-configuration interaction (FCI) densities as well as model densities, ranging up to 100 electrons and spanning both weakly and strongly correlated molecules.

The occurrence of thermal transport phenomena is widespread in nature and plays a pivotal role in the functionality of diverse electronic and thermo-electric energy-conversion devices. The traditional first-principles theory governing the thermal and thermodynamic properties of insulators relies on the perturbative treatment of interatomic potential and ad hoc displacement of atoms within supercells. However, the limitations of these approaches for highly anharmonic and weakly bonded materials, along with discrepancies arising from not considering explicit finite temperature effects, highlight the necessity for a well-defined quasiparticle approach to the lattice vibrations. To address these limitations, we present a comprehensive numerical framework in this study, designed to compute the thermal and thermodynamic properties of crystalline semiconductors and insulators. The self-consistent phonon renormalization method we have devised reveals phonons as quasiparticles, diverging from their conventional characterization as bare normal modes of lattice vibration. The extension of renormalization effects to third- and fourth-order interatomic force constants (IFCs) is also implemented and demonstrated. For the comprehensive physical insights, we employed an iterative solution of the Peierls-Boltzmann transport equation (PBTE) to determine thermal conductivity and carry out Helmholtz free energy calculations, treating anharmonic effects up to the fourth order. We utilize our numerical framework to showcase its applicability through an investigation of phonon dispersion, phonon linewidth, anharmonic phonon scattering rates, and temperature-dependent lattice thermal conductivity in both highly anharmonic materials (NaCl and AgI) and weakly anharmonic materials (cBN and 3C-SiC). Our study reveals that neglecting higher-order phonon scattering processes, particularly four-phonon interactions, is not viable for materials with strong anharmonicity in their interatomic potential. Meanwhile, renormalization demonstrates a negligible impact on materials characterized by weak anharmonicity. We also investigate the effects of fourth-order anharmonicity on the Helmholtz free energy and pressure-dependent thermal conductivity of the NaCl crystal. The theoretical and computational framework developed in this work will help to understand the physical insight of phonons and phonon-driven thermal and thermodynamic properties of materials and offer valuable guidance for the strategic development of efficient thermal management techniques.

RHEA-X is an open-source solver for the direct numerical simulation of complex multiphysics turbulent flows on modern accelerated architectures. It reconciles physical accuracy with computational efficiency, enabling simulations ranging from ideal gases and particle-laden turbulence to high-pressure transcritical fluids and external flow aerodynamics. At its core lies a modular C++ architecture built on MPI and OpenACC that combines low-dissipation spatial discretizations, strong-stability-preserving time integration, advanced thermodynamic formulations up to real-fluid equations of state and complex correlation-based transport-property models. In addition, RHEA-X incorporates immersed boundary methods and Lagrangian point particles to solve external and particle-laden flows, respectively. Through persistent GPU data residency, GPU-aware MPI communication, and portable design, RHEA-X achieves reproducible and scalable performance across diverse systems and flow configurations. Validation against canonical and multiphysics benchmarks demonstrates its ability to capture near-wall turbulence, vortex-dominated separation, and pseudo-boiling thermodynamics with accuracy consistent with reference data and published results. Performance analyses on multiple supercomputing platforms demonstrate that 97% of ideal strong-scaling behavior is maintained up to 64 GPUs, while weak-scaling efficiencies remain above 90%, and time-to-solution reductions exceeding one order of magnitude compared to CPU execution (up to approximately 70 × ). Both CPU and GPU runs exhibit consistent relative scalability, supporting deployment from academic clusters towards exascale systems. RHEA-X yields an open, extensible environment for high-fidelity computational fluid dynamics, where physical rigor, computational performance, and scientific reproducibility meet within a single, coherent framework.

Although various methods exist for modelling non-spherical particles in DEM, particles’ shapes are usually treated as immutable. However, particles often change shape gradually, e.g., due to abrasion or accrued plastic deformation. This manner of shape evolution has largely been neglected in DEM, even though it can significantly influence bulk-scale behaviour. The following introduces an extendable framework for modelling the gradual and permanent evolution of particle shapes in DEM, focusing on abrasion. By extending the LAMMPS rigid-body implementation, a comprehensive novel wear model is adapted to simulate the abrasion of arbitrarily shaped particles. Abradable particles are represented as hollow shells of discrete spheres which compute forces via standard pair interactions. These spheres form the nodes of a triangulated mesh used to calculate well-defined local surface areas and normals. Following an impact exceeding a material yield criterion, spheres are displaced inwards along their associated normals. The result is a reduction in volume and a permanent change in shape. Each abraded particle’s moment of inertia is then recomputed to resolve future rigid-body dynamics. Thus, particle-level changes in shape affect the system’s bulk dynamics, which in turn informs subsequent abrasion. Results exhibit shape evolution in agreement with a variety of abrasion scenarios in literature and showcase the consequent effect on their bulk dynamics. By linking microscale abrasion mechanisms to macroscale system behaviour, the present research has widespread applications in both natural processes and industrial particle-handling systems. Furthermore, the outlined framework can be readily adapted to other sources of mechanistic particle shape evolution in DEM.