Parameter-free calculations of lattice dynamics from first principles have achieved significant progress in the past decades, with a wealth of applications in thermodynamics, phase transitions, and transport properties of materials. Current approaches to derive the interatomic force constants (IFCs) of lattice potential become challenging and sometimes infeasible when going beyond third-order anharmonicity, due to the combinatorial explosion in the number of higher-order IFCs. In this work, we present a robust and user-friendly program, Pheasy, which reliably reconstructs the prescribed potential energy surface of crystalline solids via a Taylor expansion of arbitrarily high order. Given force-displacement datasets, the program enables an efficient and accurate extraction of IFCs using advanced machine-learning algorithms, and further calculates a wide range of harmonic and anharmonic phonon related properties. We show in three prototypical examples how the obtained IFCs have been successfully applied to study anharmonic lattice dynamics and thermal transport. Through these detailed benchmarks, we have also identified the optimal approach for IFC extractions and offered general guidelines for high-fidelity lattice-dynamical simulations, addressing the large uncertainties in the IFCs extracted from existing various schemes. Overall, the Pheasy project aims to create a phonon code ecosystem that connects diverse phonon simulation platforms and offers access to the broad research community.
