Computer Physics Communications

The GW method is a widely employed approach for calculating quasiparticle energies and band structures of materials. The static Coulomb hole plus screened exchange (COHSEX) approximation, which corresponds to the static limit of the GW self-energy, provides a computationally efficient alternative. In systems with spin-orbit coupling (SOC), it is conventional practice to treat SOC as a perturbation and apply GW corrections separately to obtain quasiparticle energies incorporating SOC effects. However, this approach may be inadequate for accurately describing these effects due to its initial neglect of spin degrees of freedom. To address this limitation, we implement the spinor COHSEX and GW formalism in KSSOLV, a MATLAB toolbox for electronic structure calculations using Kohn-Sham density functional theory. Our implementation enables meticulous treatment of spinor in COHSEX and GW calculations. We validate the correctness and computational efficiency of our approach by calculating quasiparticle energies for various molecular and periodic solids with different spin-orbit coupling strengths.

We present SPARC-atomSFE, a spectral finite-element package for accurate and efficient atomic structure calculations within the framework of Kohn-Sham density functional theory. The package supports both all-electron and norm-conserving pseudopotential calculations across a broad hierarchy of exchange-correlation approximations, spanning local, semilocal, and nonlocal functionals, within a spin-unpolarized, non-relativistic, spherically symmetric atomic framework. The nonlocal functionals include hybrid functionals and the many-body random phase approximation (RPA); for hybrid functionals, we implement both the generalized Kohn–Sham and optimized effective potential (OEP) approaches, while RPA is treated within the OEP framework. SPARC-atomSFE also includes support for fractional orbital occupations and charged atoms. Spatial discretization is based on an adaptive grid with element nodes distributed according to the Legendre–Gauss–Lobatto scheme, high-order C^0-continuous Lagrange polynomial basis functions, and Gauss–Legendre quadrature for numerical integration. We present systematic convergence studies and identify the computational parameters required to achieve target accuracies. We validate the accuracy of SPARC-atomSFE through representative calculations spanning the various exchange-correlation approximations, obtaining results that are in very good agreement with values reported in the literature. We further demonstrate two representative applications: accuracy testing of pseudopotentials for advanced exchange-correlation, and machine learning of the exact exchange OEP potential.

Accurate simulations of combustion phenomena require solving stiff ordinary differential equations (ODEs) governing chemical reactions, which often constitute a substantial fraction of the total simulation time. However, realistic chemical reaction mechanisms typically involve numerous species and exhibit high levels of mathematical stiffness. While implicit solvers with adaptive step sizes are essential for handling stiff systems, their implementation on Graphics Processing Units (GPUs) faces significant challenges due to thread divergence arising from adaptive time-stepping and error control mechanisms. As a result, existing GPU-accelerated chemical ODE solvers have predominantly focused on explicit methods or simple implicit methods, leaving a significant gap in efficient GPU-accelerated implicit solvers for stiff ODEs requiring high-order accuracy. Here we introduce GPSODE (GPU-accelerated Parallel Stiff ODE solver), a CUDA-centric package (kernels and device-side scheduling in CUDA; Python/C utilities in the repository) that implements an enhanced state-detection-based SEULEX (SD-SEULEX) algorithm specifically optimized for GPU architectures. The SD-SEULEX method enables the software to dynamically adjust the set of active grid points to mitigate GPU thread divergence while maintaining the numerical accuracy and stability of the classical SEULEX method. Floating-point work, state tagging, and index compaction for clustering run on the GPU; the host launches kernels and reads compact counters. Validation tests demonstrate excellent agreement with reference solutions for zero-dimensional homogeneous reactors using H2, CH4, and C7H16 mechanisms, as well as one-dimensional hydrogen-air premixed flames. Performance evaluations show that SD-SEULEX achieves a two-fold performance improvement compared to the baseline SEULEX algorithm. When integrated into OpenFOAM, GPSODE achieves a 172x speedup in the chemical solver compared to a single CPU core, with multi-GPU scaling tests showing 18.4x speedup in the chemical solver and 10.1x overall acceleration versus a 32-core CPU node. C-linkable entry points are provided for integration with OpenFOAM and other CFD stacks.

Phonon calculations are pivotal for understanding the thermodynamic and dynamical properties of materials, yet density functional theory (DFT) approaches often incur prohibitive computational costs for large-scale systems such as moiré superlattices. We introduce NEPHONON, an open-source high-performance package designed to quickly calculate phonon properties of extended systems with near-DFT accuracy. By leveraging the computational efficiency of machine-learned Neuroevolution potentials, NEPHONON enables the rapid generation of second-order force constants even in systems with thousands of atoms. Beyond standard phonon band structures, density of states, and group velocity, NEPHONON implements the simulation of isofrequency phonon surfaces and inelastic neutron scattering spectra S(Q,E). Furthermore, the code provides access to full phonon eigenvectors, facilitating advanced analyses of topological chiral phonons and visualization of vibrational modes. Through benchmarks on twisted bilayer phosphorene, graphene, and copper, we demonstrate the package's capabilities, confirming that NEPHONON is a powerful tool for rapidly exploring complex phononic phenomena in materials.

We present DiffMoM.jl, an open-source Julia package for differentiable electric-field integral-equation (EFIE) method-of-moments (MoM) analysis and surface-impedance inverse design of single-layer reactive metasurfaces with Rao–Wilton–Glisson (RWG) basis functions. The implementation exposes a reproducible forward–adjoint pipeline in which impedance sensitivities are assembled analytically and scalar radiation objectives are differentiated without differentiating through the linear solver. Quadratic objectives require one adjoint solve, directivity-ratio objectives require two, and full central finite-difference gradients require \(O(P)\) additional forward solves for \(P\) design parameters. Verification combines energy-balance checks, mesh-sweep adjoint/finite-difference agreement, and an independent Bempp-cl comparison. Across the tested meshes, the energy-balance error remains within 2%, the mesh-sweep gradient error remains below 6 x 10^-5, and beam-centric Bempp-cl agreement stays within 0.38 dB on main-beam level and 0.47 dB on sidelobe suppression across seven reactive cases. A 4λ x 4λ, 144-cell beam-steering benchmark raises the target-angle gain by about 23 dB relative to a perfect electric conductor within the stated dense-solver scope.

Condensed matter compounds typically form crystals, which break the rotational and translational invariance of space but remain invariant under a discrete set of symmetry operations. Understanding the effects allowed by this symmetry breaking, as well as the constraints imposed by the crystal structure, is a crucial problem in condensed matter physics. Here, we present a Python package for determining the symmetry-restricted forms of tensors describing physical properties of crystals, focusing particularly on magnetic materials. The primary focus is on response tensors; however, the program can also describe equilibrium properties and other physical properties, such as magnetic interactions. The program can describe the symmetry using the conventional magnetic space groups, as well as using the spin groups that describe the non-relativistic limit. Additional functionality includes the treatment of quantities projected onto a particular site and expansions in the magnetic order parameter. The code can be used either from the command line or via a Python API.

The rapidly increasing computational demand of large-scale density functional theory (DFT) simulations poses significant performance and scalability challenges for traditional CPU-based supercomputers. Recent advances in GPU-like accelerator architectures, featuring massive on-chip parallelism and high memory bandwidth, are reshaping the computational landscape of large-scale first-principles calculations. In this work, we port the plane-wave DFT software PWDFT to a GPU-like deep computing unit (DCU) heterogeneous platform on the Sugon supercomputer and achieve significant acceleration for hybrid functional calculations through low-rank approximations and high-performance computing optimizations. We integrate the adaptively compressed exchange (ACE) approximation with the projected commutator direct inversion of the iterative subspace (PCDIIS) scheme, thereby reducing the computational time per self-consistent step to only about four times that of a corresponding PBE calculation under the benchmark settings. We further exploit DCU hardware acceleration through optimized linear algebra and FFT libraries, as well as communication-aware kernel optimizations. We demonstrate that plane-wave hybrid DFT calculations for a 2,048-atom silicon system can be performed using only 32 DCUs, and that the framework efficiently scales to employ 4096 DCUs for simulations of a 4,096-atom system. These advances establish an efficient computational framework for large-scale plane-wave hybrid functional calculations at the thousand-atom scale.

We introduce DiFfRG (Discretisation Framework for functional Renormalisation Group flows), a comprehensive computational C++ framework for solving functional Renormalisation Group flows in very general truncation schemes. Its central features are threefold: Firstly, the use of Finite Element Methods (FEM) for efficient, easy to set up and quantitatively reliable computation of field dependences. Secondly, the (simultaneous) setup of large, fully momentum-dependent vertex expansions. Thirdly, efficient time-discretisation methods incorporating insights from studies of solving theories which exhibit spontaneous symmetry breaking, going hand in hand with shocks and an exponential increase of the information flow velocity in field space. The framework provides a Mathematica package for automatic code generation of flow equations, finite and zero temperature integration routines with support for GPU hardware and extensive parallelisation capabilities. Detailed examples and tutorials are provided and discussed herein which showcase and introduce the framework to the user. We illustrate the capabilities of DiFfRG with four examples, with the complete codes fully included: finite temperature O(N) theory, a Quark-Meson model, SU(3) Yang-Mills theory and four-Fermi flows in the QCD phase diagram.

Trigonal solid-state defects are often subjects of spontaneous symmetry breaking driven by the E⊗e Jahn-Teller effect, reflecting strong electron-phonon coupling. These systems, particularly paramagnetic defect qubits in solids are central for quantum technology applications, where accurate knowledge of their fine-structure parameters – shaped by the complex interplay of spin-orbit and electron-phonon interactions – is essential. We introduce the jahn - teller - dynamics package, a Python code that implements the first-principles approach of [Phys. Rev. X 8, {021063} (2018)] to accurately compute the spin-orbit-phonon entanglement in trigonal defects utilizing the output from density functional theory calculations (DFT) to predict fine-structure parameters of zero-phonon lines (ZPLs), including Zeeman shifts under external magnetic fields. We demonstrate its capabilities on negatively charged Group-IV–vacancy (G4V) defects in diamond: SiV-, GeV-, SnV-, PbV- and the neutral N3V0 defect in diamond, and the CH3O0 methoxy radical. Additionally, we implement a generic electron-phonon code that is capable entangling an arbitrary amount of (i) vibration modes and (ii) electronic levels by (iii) arbitrarily high order of vibronic interaction terms. Exemplarily, we demonstrate the Jahn-Teller multimode problem on the CH3O0 methoxy radical and the Pseudo Jahn-Teller case on the C4H4+ butatrien cation.

We present KnudsenModel, an open-source Python implementation of the Knudsen model for computing evaporative cooling of microscopic liquid jets in vacuum. The code builds on widely used scientific Python packages and provides a transparent and flexible framework that can be applied to both cylindrical jets and spherical droplets. A fixed-mass radial discretization is employed to consistently account for the coupled evolution of mass loss and density, ensuring physically accurate predictions of both temperature and jet or droplet diameter. We analyze the numerical behavior of the method, including convergence properties and comparisons with literature calculations based on alternative discretization schemes, revealing that, in particular for water droplets, the present formulation captures a shallow minimum in the diameter that reflects the interplay between evaporation-induced mass loss and water's density anomaly—a feature absent in previous results.