Efficient systematic scheme to construct second-principles lattice dynamical models

We start from the polynomial interatomic potentials introduced by Wojdeł et al. [J. Phys.: Condens. Matter 25, 305401 (2013)] and take advantage of one of their key features—namely, the linear dependence of the energy on the potential's adjustable parameters—to devise a scheme for the construction of first-principles-based (second-principles) models for large-scale lattice-dynamical simulations. Our method presents the following convenient features. The parameters of the model are computed in a very fast and efficient way, as it is possible to recast the fit to a training set of first-principles data into a simple matrix diagonalization problem. Our method selects automatically the interactions that are most relevant to reproduce the training-set data, by choosing from a pool that includes virtually all possible coupling terms, and produces a family of models of increasing complexity and accuracy. We work with practical and convenient cross-validation criteria linked to the physical properties that will be relevant in future simulations based on the new model, and which greatly facilitate the task of identifying a potential that is simultaneously simple (thus computationally light), very accurate, and predictive. We also discuss practical ways to guarantee that our energy models are bounded from below, with a minimal impact on their accuracy. Finally, we demonstrate our scheme with an application to ferroelastic perovskite ${\mathrm{SrTiO}}_3$, which features many nontrivial lattice-dynamical features (e.g., a phase transition driven by soft phonons, competing structural instabilities, highly anharmonic dynamics) and provides a very demanding test.

Figure 1

Sketch illustrating how the choice of a supercell for the generation of the TS data implicitly imposes a spatial cutoff for the interatomic interactions. The elemental unit cell of the crystal is colored in dark grey, while the supercell used for the TS calculations is in light grey. The elemental cell contains two atoms. The atoms in the supercell are displayed in strong colors, while their periodic images are shaded. The interactions labeled as “Int. 1” and “Int. 2” are obviously different, but they cannot be resolved on the basis of data obtained for this supercell. Note that the longest-range “Int. 2” connects atoms that are also linked by “Int. 1” due to the periodic supercell repetition.

Figure 2

Sketch of the STO perovskite structure. (a) shows the elemental cell for the high-symmetry cubic ($Pm\overline{3}m$) structure. (b) sketches the antiphase ${\mathrm{O}}_6$ rotations that characterize the AFD distortion driving the structural phase transition in STO.

Figure 3

Eigenvalues obtained from the diagonalization of the harmonic part of the phonon energy $E_{\mathrm{p}}$, at two specific $q$ points $\mathrm{\Gamma }$ and $R$ (see text). The solid lines indicate the results obtained from the fitted model as a function of the number of terms included, while the exact DFT-computed results are given by dashed lines.

Figure 4

Energy evolution as a function of time corresponding to the typical MD run used to construct our TSs of DFT data. This particular case corresponds to the run used to build TS@10 (see text). The run starts from the RS, which is taken to be the zero of energy. Eventually, the simulated system finds its way to the ground-state structure and fluctuates around it.

Figure 5

Optimized goal function, as a function of the number of terms, for the various models considered in this work. Red and black lines correspond to FM and EHM models, respectively. Solid, dashed, and dotted lines correspond to models fitted to the TS@10, TS@300, and TS@10+300, respectively. See text for details.

Figure 6

Behavior of the goal function evaluated using fitted models and test-set configurations, as described in the text. The colors and line types are as in Fig. 5.

Figure 7

Behavior of the goal function evaluated using fitted models and test-set configurations, as described in the text. The colors and line types are as in Fig. 5.

Figure 8

Performance of the six models discussed in this work as regards the prediction of various ground state [(a)–(d)] and MD [(e)–(j)] related quantities. To obtain the data in (e)–(j), we use our model potentials to compute the energy of the configurations visited in the MD trajectories computed from DFT at various temperatures; we report the average difference between the DFT energy and that obtained from the model.

Figure 9

Comparison of the DFT-computed energy for MD trajectories at various temperatures and the energies obtained by evaluating our optimum model potentials for the corresponding structures. Except in the 10 and 500 K cases, the deviations in the energy are all but invisible in the scale of the figure, emphasizing the good overall accuracy of our fitted models.

Figure 10

Values of the three parameters that our fitting procedure selects as most relevant to construct the EHM fitted to TS@10. The figure shows their evolution as a function of the total number of terms included in the model. Note that all three terms are linear in strain and quadratic in the atomic distortions; hence, they all have the same units, as indicated in the figure.

Figure 11

Values of the parameters obtained for our FM fitted to TS@10 and ordered as a function of the longest interatomic separation associated to the corresponding coupling. Note that this figure includes information about coefficients corresponding to different orders of our Taylor series, which thus have different units; hence, we indicate arbitrary units (a.u.) and stress that this figure is to be taken only as a qualitative illustration of the spatial decay of the interactions.

Figure 12

Sketch of some of the most important interactions identified by our automatic fitting procedure. The atomic displacements corresponding to a representative term of the SAT are shown; strains are indicated with dark triangles. The circled numbers indicate the order of specific displacement-difference terms. Note that all the sketched interactions play a role to control the energetics of AFD modes and polar distortions, as described in the text.

Figure 13

Computed evolution of the AFD order parameter, as a function of temperature, for our TS@10 and TS@10+300 models. We also show the result for a pressure-corrected version of the FM model fitted to TS@10 (see text).