Sparse regression fits crystal force constants to model heat flow and vibrations

The work

Technical · from the record

Pheasy is a Python package that reconstructs the potential energy surface of a crystal as a symmetry-reduced Taylor expansion and extracts interatomic force constants from force-displacement data by ordinary least squares, LASSO or adaptive LASSO regression. The authors fitted sixth-order lattice models for bulk silicon, monolayer WS2 and cubic SrTiO3 on ensembles of displaced supercells with forces from density functional theory, testing force prediction for silicon on an independent ensemble of 64 structures with a root mean square error of 3.239 meV/A. The extracted cubic and quartic force constants were used for self-consistent phonon calculations of temperature-dependent dispersions and for lattice thermal conductivity with three- and four-phonon scattering; for silicon the reported room-temperature value is 125.4 W/mK with three-phonon scattering and 115.4 W/mK once four-phonon scattering is included. Comparing fitting schemes on cubic SrTiO3, the authors report that dispersions from fourth-order force constants fitted together with all other orders agree with their SCHA reference while those from a scheme that fixes the harmonic terms separately deviate, and they recommend fitting all orders at once.

How AI was used

Initial harmonic force constants were obtained from a small number of randomly displaced supercells or from density functional perturbation theory, and used to build a displacement covariance matrix from which training and test ensembles of displaced supercells were drawn by quantum canonical sampling at 300 K; forces for each configuration were computed with density functional theory in Quantum ESPRESSO. The cluster expansion of the lattice potential was reduced to its irreducible parameters using space-group, permutation, translational and rotational invariance constraints via an iterative maximal-pivot Gauss-Jordan null-space construction, and the displacement polynomials were assembled into a sensing matrix whose columns were standardised to zero mean and unit variance. Force constants up to sixth order were then obtained by regressing the forces on this matrix with ordinary least squares, LASSO from scikit-learn, or an adaptive LASSO whose feature weights were initialised from the inverse of the least-squares solution, with cutoff radii and many-body restrictions applied to the high-order terms. The fitted model was evaluated by predicting forces for an independent ensemble generated the same way, and the resulting force constants were passed to self-consistent phonon and self-consistent harmonic calculations and to Boltzmann transport solvers for thermal conductivity; for temperature-dependent anharmonic force constants the ensembles were re-sampled and refitted with the harmonic terms fixed at the self-consistent phonon solutions.