PARPHOM: PARallel PHOnon calculator for Moiré systems

The introduction of a twist between two layers of two-dimensional materials has opened up a new and exciting field of research known as twistronics. In these systems, the phonon dispersions show significant renormalization and enhanced electron-phonon interactions as a function of the twist angle. However, the large system size of the resulting moiré patterns in these systems makes phonon calculations computationally challenging. In this paper, we present PARPHOM, a powerful code package designed to address these challenges. PARPHOM enables the generation of force constants, computation of phononic band structures, and determination of density of states in twisted 2D material systems. Moreover, PARPHOM provides essential routines to investigate the finite temperature dynamics in these systems and analyze the chirality of the phonon bands. This paper serves as an introduction to PARPHOM, highlighting its capabilities and demonstrating its utility in unraveling the intricate phononic properties of twisted 2D materials.

💡 Research Summary

**
The manuscript introduces PARPHOM, a parallel software framework designed to compute phonon properties of twisted two‑dimensional (2D) materials—so‑called moiré systems—where the unit cells can contain thousands to tens of thousands of atoms. Traditional density‑functional‑theory (DFT) phonon codes are limited to a few hundred atoms because the construction of force‑constant matrices and the subsequent diagonalization of the dynamical matrix become prohibitively expensive in memory and CPU time. PARPHOM circumvents these limitations by exploiting two key ideas: (i) a massively parallel generation of force constants using classical force fields within LAMMPS, and (ii) a ScaLAPACK‑based parallel diagonalization of the dynamical matrix.

Force‑constant generation
The workflow starts from a relaxed atomic configuration obtained with a classical potential in LAMMPS. A Python script (get_force_constant.py) displaces each atom in the moiré unit cell by a small amount (default Δr = 10⁻⁴ Å) along ±x, ±y, ±z, computes the forces on all atoms for each displacement, and evaluates the second‑derivative force constants via central finite differences (Eq. 2). Because each atom’s six displacements are independent, the task is “embarrassingly parallel.” The implementation distributes the n atoms over Nₚ MPI ranks, each rank performing its subset of displacements and writing the resulting 4‑dimensional tensor Φαβ(j0, j′0) (size n × n × 3 × 3) to an HDF5 file using parallel‑h5py. Unlike PHONOPY, PARPHOM does not reduce the number of displacements by symmetry; the authors argue that the parallel speed‑up outweighs the modest savings from symmetry reduction. After all displacements are completed, the code iteratively enforces two fundamental symmetries: (a) index (transposition) symmetry Φαβ(jl, j′l′)=Φβα(j′l′, jl) and (b) the acoustic sum rule ∑j′l′Φαβ(jl, j′l′)=0. The iteration proceeds until convergence or a maximum number of cycles is reached. Rotational point‑group symmetries are not imposed because low‑angle relaxed moiré structures typically lose D₆/D₃ symmetry, leaving only time‑reversal symmetry as identified by Spglib.

Dynamical matrix construction and diagonalization
Given the force‑constant tensor, the dynamical matrix Dαβ(jj′, q) is assembled for each user‑specified q‑point using Eq. 7–8. The authors note that LAMMPS already incorporates the minimum‑image convention, so the phase factor f(q) automatically includes contributions from neighboring cells, eliminating the need for an explicit sum over supercells. The dynamical matrix is distributed in a block‑cyclic layout across the MPI grid and diagonalized with ScaLAPACK’s PZHEEVX routine, which can compute either the full spectrum or a selected energy window. Users may request eigenvalues only or both eigenvalues and eigenvectors; the latter are stored in HDF5 as a compound datatype (real and imaginary parts). The parallel diagonalization scales efficiently to hundreds of cores, as demonstrated by benchmark comparisons with PHONOPY and PhonoLAMMPS.

Group velocities and chirality
Group velocities vqν=∇q ων are essential for transport and chirality analysis. PARPHOM offers two strategies: (1) an analytic derivative of the dynamical matrix with respect to q (Eq. 12), and (2) a central finite‑difference approach using D(q ± Δq/2) (Eq. 13). The analytic method requires the vector rj(0) − rj′(l) and the phase factor, while the finite‑difference method uses a small q‑step Δq = 10⁻⁴ Å⁻¹. The velocity matrix elements are then obtained via Eq. 14, which avoids the numerical instability that plagues simple finite‑difference of ω(q) near band crossings. The code also computes a “chirality” indicator by examining the rotational character of each eigenmode, enabling studies of non‑reciprocal phonon transport in twisted bilayers.

Temperature‑dependent phonon DOS
Beyond zero‑temperature spectra, PARPHOM can evaluate temperature‑dependent phonon density of states (DOS) using mode‑projected velocity autocorrelation functions (VACF). The user may choose Gaussian, Lorentzian, or linear‑triangulation smearing to construct the DOS from the discrete eigenfrequencies. This capability bridges the gap between static harmonic calculations and finite‑temperature experimental observables such as Raman linewidths.

Software architecture and usability
PARPHOM is released under GPL‑v3 on GitHub (https://github.com/qtm‑iisc/PARPHOM). It is written in a hybrid of Fortran (core parallel routines) and Python (pre‑ and post‑processing). Required dependencies include NumPy, mpi4py, h5py, and the serial LAMMPS Python wrapper. Build scripts (Makefile.intel, Makefile.aocc) link against MPI, parallel HDF5, and ScaLAPACK, supporting both Intel and AMD compiler suites. A supplemental Python package, pyphutil, provides convenient functions for reading the HDF5 output, plotting phonon band structures, DOS, and group velocities. The output file follows a clear hierarchy: each q‑point has a group containing datasets /qvec, /eigenvalues, optional /evec, and optional /vel, facilitating downstream analysis.

Performance benchmarks
The authors benchmark PARPHOM on a twisted bilayer graphene supercell containing ~10 000 atoms. Force‑constant generation scales nearly linearly with the number of MPI ranks, achieving a speed‑up of ~80× on 128 cores compared with the serial PHONOPY approach. The subsequent ScaLAPACK diagonalization of a dynamical matrix of dimension 30 000 × 30 000 (3 × 10 000 degrees of freedom) completes in under 5 minutes on 256 cores, whereas a conventional LAPACK routine would be infeasible due to memory constraints. Memory usage remains bounded because the force‑constant tensor is stored in parallel HDF5 and only the necessary block of the dynamical matrix resides in each process’s memory.

Conclusions and outlook
PARPHOM delivers a complete, scalable solution for phonon calculations in moiré superlattices, addressing the two most critical bottlenecks: (i) the generation of massive force‑constant matrices and (ii) the diagonalization of correspondingly large dynamical matrices. By leveraging classical force fields, embarrassingly parallel finite‑difference displacements, and ScaLAPACK’s high‑performance linear algebra, the package makes it possible to explore twist‑angle‑dependent phonon renormalization, temperature effects, and chirality in systems previously out of reach. Future extensions could incorporate electron‑phonon coupling, anharmonic perturbation theory, or machine‑learning‑based force‑constant prediction, further expanding the scientific impact of this framework.