Differentiable Gaussian Process Force Constants

The finite temperature properties of matter requires understanding thermal motion. In crystals this
can be described by representing the collective excitation of the centre of mass motion as a phonon.
Phonon frequencies (Energies) in the harmonic approximation come from the second order force con-
stants (FC) of the potential energy surface (PES). The standard approach is to use a finite displacement
method (FDM). Anharmonic contributions (required for finite thermal conductivity) require higher order
force constants. FDM-based calculations scale poorly both with the size of the system, and the order of
the force constants[1]. As an alternative approach, we can use a more sophisticated surrogate potential
energy surface model [2].
Gaussian processes (GPs) are a supervised machine learning method which can describe an arbitrary
function. GPs are naturally Bayesian (probabilistic). Our reference data is the electronic structure from
which conditions the model (training) [3, 4].
Due to the underlying Gaussian form of a GP, the model is infinitely differentiable [5, 6]. This allows
the model to be trained directly on forces (the derivative of energy), reducing the number of calculations
required for a given accuracy of PES evaluation [7-9]. This differentiation can be extended to compute
the second and the third derivative of PESs (harmonics and cubic anharmonic FCs) by using automatic
differentiation (AD). By performing linear operations between arbitrary derivative orders of the GP, the
covariance functions among PESs, forces and those FCs can be calculated.
We implement this method in the Julia language, in our GPFC.jl package. We first use our tech-
nique of anharmonic property calculations of Si, NaCl and PbTe which their atomic environments are
represented on atomic Cartesian coordinates. To impose their space group symmetry, phonon coordinate
representations are then used to describe the environment of these three materials resulting in faster con-
vergence of anharmonic properties. To further impose a three-dimensional rotation symmetry, we aim to
use an SO(3) representation (such as spherical harmonics) with our GP method to achieve accurate FCs
with less data [3, 10]. We compare our results to standard approach of FDM (in Phono(3)py.py) [2]
and cluster expansion (HiPhive.py) [11].
References
[1] A. Togo, L. Chaput and I. Tanaka, Distributions of phonon lifetimes in brillouin zones, Physical Review B 91 (2015) .
[2] A. Togo et al., Implementation strategies in phonopy and phono3py, Journal of Physics: Condensed Matter 35 (2023) 353001.
[3] A.P. Bartók et al., Gaussian approximation potentials: A brief tutorial introduction, International Journal of Quantum Chemistry 115 (2015) 1051.
[4] H. Sugisawa et al., Gaussian process model of 51-dimensional potential energy surface for protonated imidazole dimer, The Journal of Chemical Physics 153 (2020) 114101.
[5] E. Solak et al., Derivative observations in Gaussian Process Models of Dynamic Systems, NIPS’02: Proceedings of the 15th International Conference on Neural Information Processing Systems (2002) 8.
[6] A. McHutchon, Differentiating Gaussian Processes.
[7] E. Garijo del Río et al., Local Bayesian optimizer for atomic structures, Physical Review B 100 (2019) 104103.
[8] S. Kaappa et al., Global optimization of atomic structures with gradient-enhanced gaussian process regression, Physical Review B 103 (2021) .
[9] K. Asnaashari et al., Gradient domain machine learning with composite kernels: improving the accuracy of PES and force fields for large molecules, Machine Learning: Science and Technology 3 (2021) 015005.
[10] V.L. Deringer et al., Gaussian process regression for materials and molecules, Chemical Reviews 121 (2021) 10073.
[11] F. Eriksson et al., The hiphive package for the extraction of high-order force constants by machine learning, Advanced Theory and Simulations 2 (2019) .