{'fraction_non_alphanumeric': 0.03590492779681969, 'fraction_numerical': 0.019497668146316795, 'mean_word_length': 4.944555778223113, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Structural Phase Transitions in SrTiO3 from Deep Potential Molecular Dynamics

Strontium titanate (SrTiO3) is regarded as an essential material for oxide electronics. One of its many remarkable features is subtle structural phase transition, driven by antiferrodistortive lattice mode, from a high-temperature cubic phase to a low-temperature tetragonal phase. Classical molecular dynamics (MD) simulation is an efficient technique to reveal atomistic features of phase transition, but its application is often limited by the accuracy of empirical interatomic potentials. Here, we develop an accurate deep potential (DP) model of SrTiO3 based on a machine learning method using data from first-principles density functional theory (DFT) calculations. The DP model has DFT-level accuracy, capable of performing efficient MD simulations and accurate property predictions. Using the DP model, we investigate the temperature-driven cubic-to-tetragonal phase transition and construct the in-plane biaxial strain-temperature phase diagram of SrTiO3. The simulations demonstrate that strain-induced ferroelectric phase is characterized by two order parameters, ferroelectric distortion and antiferrodistortion, and the ferroelectric phase transition has both displacive and order-disorder characters. This works lays the foundation for the development of accurate DP models of other complex perovskite materials.


I. INTRODUCTION
Perovskite oxides with formula ABO3 is an important class of functional materials.
In the ABO3 perovskite structure, the A site cation is twelve-fold coordinated and the B site cation is six-fold coordinated with oxygen anions, respectively. Many different cations and anions can be incorporated in this structure, resulting in a great number of functional materials that hold great promises for new electronic device applications [1][2][3]. Strontium titanate (SrTiO3, STO) as a representative ABO3 perovskite oxide has long captured the attention of condensed matter physicists, because of its rich physical properties such as two-dimensional electron gas [4][5][6], Rashba spin-orbit coupling [7][8][9], interfacial superconductivity [10,11] and multiple structural instabilities [12]. Among them, the incipient ferroelectricity is a unique structural property of STO that has attracted extensive studies [13][14][15]. At room temperatures, STO has a cubic structure with Pm-3m space group. Below 105 K, it undergoes an antiferrodistortive (AFD) structural transition, transforming to a tetragonal I4/mcm phase. The AFD structure can be described as neighboring TiO6 octahedra rotating a small angle in the opposite directions around [001] axis. The tetragonal phase is paraelectric (i.e. incipient ferroelectricity), where the quantum fluctuations and AFD instabilities suppress longrange polarization ordering. Despite the apparent paraelectric nature of STO, it exhibits many abnormal phenomena indicating its proximity to a ferroelectric phase [13,14].
The tetragonal-to-cubic phase transition in STO has been intensively studied with thermodynamic Landau-Ginsburg-Devonshire (LGD) model [23][24][25][26]. Nevertheless, the LGD model is often subjected to two well-known limitations: i). the material coefficients involved in the Landau free energy expansion are difficult to determine, and a slight variance of coefficient value may cause a drastic change of the resulting phase diagram [25]; ii). such phenomenological models cannot reveal the mechanisms of phase transition at the atomic level. First-principles density functional theory (DFT) calculations can accurately predict the atomic structures of cubic and tetragonal phases of STO and their small energy difference (~1.1 meV/atom) as well as the quantitative relationship between strain and polarization [27,28], but the expensive computational cost prevents the application of DFT from studying STO phase transition over a larger length and/or time scale (e.g. limited to a system of hundreds of atoms in the time scale of pico-second). Classical molecular dynamics (MD) simulation is a very effective technique for understanding the structural phase transition mechanisms. Several empirical interatomic potentials have been developed to simulate the cubic structure of STO [29][30][31]. However, none of them has been used to explore the cubic-to-tetragonal phase transition, likely because the subtle energy variation (~ 1 meV/atom) induced by the structural transition is beyond the energy resolution of the analytical function used to construct the interatomic potential [29]. Therefore, accurate MD simulations of phase transition in STO demands a high-accuracy potential.
In recent years, machine learning methods have been used as a powerful tool to develop the interatomic potential of crystalline and molecular materials. For instance, the Behler-Parrinello neural network [32], the Gaussian approximation potentials [33], the Select Neighbors and Parameters [34] and the SchNet [35] have achieved great success. Among them, the recently proposed Deep Potential (DP) model based on deep neural network can provide a DFT-level accurate interatomic potential by training from a large dataset of DFT energies and forces of a wide range of atomic configurations [36,37]. The well-trained deep-neural-network potential can be used to perform MD simulations without sacrificing the accuracy. In recent works, it has been demonstrated that the energy error of a DP model relative to DFT is within milli-electron-volt per atom for the ground-state configuration [38][39][40][41][42]. Whether the DP model could describe the subtle free energy variation of just ~1.1 meV/atom during tetragonal to cubic phase transition in STO remains a question.
In this work, we develop an accurate DP model for STO base on a deep learning method trained on the data from DFT calculations. The accuracy of the DP models is examined by comparing with DFT results. The temperature-driven tetragonal to cubic phase transition of STO is well captured by deep potential MD (DPMD) simulations.
Finally, the ferroelectric-paraelectric phase transition temperatures at different in-plane strains are calculated by the DP model, and the nature of ferroelectric phase transition is characterized at the atomic level. These DPMD simulation results provide atomistic insights into the microscopic ferroelectricity mechanism of STO under strain.

A. Deep learning of interatomic potential of SrTiO3
In DP model, the potential energy E of a configuration is assumed to be a sum of each atomic energy Ei of atom i, which is mapped from a descriptor Di through an embedding network. The descriptor Di characterizes the local environment of atom i within a cutoff radius Rc. Here, the Rc is set to 6 Å. The maximum number of atoms within the Rc is set to 160 for Sr, 160 for Ti, and 480 for O, respectively. The translational, rotational, and permutational symmetry of the Di are preserved by an embedding network. The smooth edition of deep potential model was employed to remove the discontinuity introduced by the cutoff radius [37]. The size of embedding network and fitting network is (25, 50, 100) and (240, 240, 240), respectively. Residual network was employed in the training of fitting network. The loss function is defined as:

B. Generation of the training data set
In this work, we use the deep potential generator (DP-GEN) to generate a set of training data that covers a sufficient wide range of relevant configurational space efficiently [39]. DP-GEN is a concurrent learning procedure and the workflow of each iteration includes three main steps: training, exploration, and labeling (see Fig. 1).
We start with DFT relaxed ground-state structures of Pm-3m cubic and I4/mcm tetragonal phase. The lattice constants calculated by DFT with Perdew-Burke-Ernzerhof solid (PEBsol) are shown in Table I. In the first iteration, the initial training dataset contains 200 randomly perturbed structures of 2 × 2 × 2 Pm-3m cubic and 1 × 1 × 1 I4/mcm tetragonal cell. Starting with these DFT training data set, four DP models are trained starting with different values of deep neural network parameters. Then, the exploration step was performed in which one of the DP models is used for MD simulations to explore the configuration space at different temperatures. For all sampled configurations in MD trajectories, the other three DP models will predict the corresponding energies and atomic forces. The maximum deviation of the four DP model prediction of forces (σ ) can be used to formulate the criterion for labeling: where 〈… 〉 indicates the average of model predictions. The configurations with σ < are already well represented by previous training dataset (green points in Fig. 1e), which are labeled as accurate configurations. Whereas a configuration with a large model deviation σ > ℎ ℎ is highly distorted and unphysical due to poor model quality (red points in Fig. 1e), and is thus labeled as failed configuration.
Therefore, only configuration satisfying < σ < ℎ ℎ are labeled as candidates for further DFT calculations and are added to the training dataset for training in the next iteration (blue points in Fig. 1e). Here, and ℎ ℎ are set to 0.05 eV/Å and 0.20 eV/Å respectively as suggested in the original literature [39]. It is important to point out that σ can be also used as the convergence criterion for DP-GEN iterations. Here, the exploration of each system is considered converged when the percentage of accurate configurations (σ < ) is larger than 99 %. For more details of DP-GEN method, please refer to the original literature [22].

C. DFT calculations setting
The initial training dataset is obtained by performing a 10-step ab initio MD simulation for randomly perturbed structures at 50 K. After labeling candidate configurations, selfconsistent DFT should be performed. All DFT calculations were performed using a plane-wave basis set with a cutoff energy of 500 eV as implemented in the Vienna Ab initio Simulation Package (VASP) [43,44], the electron exchange-correlation potential was described using the generalized gradient approximation and PBEsol scheme [45].

D. MD simulations setting
The MD simulations were carried out using LAMMPS code with periodic boundary conditions [46]. At the exploration step, MD simulations adopt the isobaric-isothermal (NPT) ensemble with temperature set from 50 to 400 K, pressure set from 0.001 to 50 kbar, and time with 2 to 20 ps (Table S1). The Nose-Hoover thermostat and Parrinello-Rahman barostat are employed to control temperature and pressure, respectively [47,48]. After DP-GEN iteration is considered converged, the optimized DP models of STO can be used to study the structural transitions driven by temperature and strain via performing MD simulations. The NPT MD simulations starting with a ground-state tetragonal I4/mcm 10  10  10 pseudo-cubic supercell (5000 atoms) were performed with periodic boundary conditions and at the zero pressure to model the temperature and strain-driven structural transition. The time step is set to 0.001 ps. Simulations using 20 × 20 × 20 supercell give similar results. At a specified temperature, the equilibrium run is 10 ps, followed by a production run of 90 ps (see Fig. S2 [49]). The main results are not sensitive to the initial configuration (Pm-3m or I4/mcm) for MD simulation. To apply a strain of εs, the in-plane lattice constant of the supercell was fixed at a0+εsa0 Å, where a0 is the optimized in-plane lattice constant of STO at the corresponding temperature, and the out-of-plane lattice constant was allowed to fully relax.

E. DP-GEN iteration process
The distributions of σ at different temperatures in the first five iterations are shown in Fig. S1 [49], and the detailed percentage values of the accurate, candidate, and failed configurations in each iteration are shown in Table S1 [49]. In the first iteration, starting with two systems of 2  2  2 Pm-3m supercell (40 atoms) and 1  1 (20 atoms iterations. All settings of the exploration strategy are listed in Table S1 [49].

A. Accuracy of DP model
The whole deep learning process is parameter-free, which ensures that the generality and transferability of DP models. We compare the energies and atomic forces calculated using the DP model and DFT calculations for the configurations in the final training dataset. As shown in Fig. 2, the mean absolute error of energy (ΔE DP-DFT ) and atomic force (ΔF DP-DFT ) between DP and DFT is 0.322 meV/atoms and 0.025 eV/Å,    The results are shown in Table III, Fig. 6. It is indicated that the tetragonal phase has a longrange AFD ordering pattern at low temperature, and that the long-range feature vanishes at higher temperature where we observed short-range AFD ordering patterns appearing dynamically in space at 200 K (see Fig. 6c). According to the variation of lattice constant and φ, the transition from tetragonal to cubic phase takes place when the temperature is raised over 190 K. The DP model predicted behavior of temperature driven phase transition in STO is in qualitative agreement with experiment, but the transition temperature is slightly higher than the experimental value of 105 K [51]. In addition, DPMD predicts a linear behavior of a and b-axis lattice constants at low temperature, whereas the saturated temperature dependence of the lattice parameter below 50 K was observed in the experiment [52]. The overestimate of transition temperature and linear behavior of lattice constants at low temperature is likely a consequence of classical limitation of standard MD simulation. This is due to quantum effect plays an important role at low temperatures [14,15], but it is neglected in the standard MD simulations. Adding the statistical quantum features into MD could eliminate the disagreement between experiment and DPMD simulation [53][54][55].

C. Strain effect
It is well known that STO is an incipient ferroelectric so that the ferroelectric transition does not appear under strain-free state even when temperature decreases to 0 K. It is generally believed the ferroelectric order is suppressed by quantum fluctuation [14,15] and preceding AFD distortion [12,56]. No matter what type of the reasons it was, that the situation may change dramatically in the strained state. Many theoretical works and experiments reported epitaxial strain could induce a polar or ferroelectric phase in STO films [17,[24][25][26]28,50]. Here, we predict the in-plane biaxial straintemperature phase diagram of STO using DPMD simulations. Fig. 7 shows the phase diagram of ferroelectric (FE) and AFD phase as a function of the temperature and in-plane biaxial strain for STO by DPMD. Here, we use two order parameters, the average polarization (P) and the AFD order parameter (φ), to distinguish the pure FE, FE, AFD, and paraelectric phase: the pure FE phase is defined as P ≠ 0 and φ = 0, the FE phase is defined as P ≠ 0 and φ ≠ 0, the AFD phase is defined as P = 0 and φ ≠ 0, whereas paraelectric phase is defined as P = 0 and φ = 0. It is seen that FE phase become most stable phase when the compressive (only when large than 0.2 %) and tensile strain was applied. At the low temperatures, the compressive strain can induce the out-of-plane P and AFD structure with rotation axis along [001], whereas the tensile strain induces the in-plane P and AFD structure with rotation axis along [100] and [010] (structural details see Fig. S3 and S4 [49]). Then the FE phase transfer to AFD structural phase when the temperature increases. It is clear to see that the ferroelectric transition temperature increases with the in-plane strain. When further increasing temperature, the AFD structure phase will transform to the paraelectric phase with P and φ both decreased to zero. The DP model shows that the ferroelectric transition temperature is ~110 K at 0.8 % compressive strain and ~210 K at 0.8 % tensile strain. DPMD simulations clarify that the wide region of possible ferroelectric transition range in the temperature-strain phase diagram predicted by LDG model is the AFD phase [11,13]. In addition, it is interesting to note that the previously reported pure FE phase (P ≠ 0 and φ = 0) was not observed in this our phase diagram [10,11]. These results are in line with theoretical and latest experiment studies that ferroelectric phase transition in STO under strain has the mixture of displacive and order-disorder characteristics [17,57,58]. Similar characteristics of phase transition are also found in tensilely strained STO (see Fig. S5 [49]).

IV. CONCLUSIONS
In summary, we developed DP model to describe dynamic and energetic properties, and structural transition in STO with the accuracy of DFT at the PBEsol level. After confirming the accuracy of the DP models by comparison of lattice constant, free energy and other properties with DFT calculation results, we use the DP model to simulate temperature driven tetragonal to cubic phase transition at the finite temperature.
Then, we predict the in-plane biaxial strain-temperature phase diagram of STO using DP model and get the atomistic insights into the ferroelectric and AFD structural transition. We clarify that ferroelectric (FE) phase at low temperature under strain is associated with AFD. With the temperature increasing, such FE phase transform into AFD phase. Moreover, it is demonstrated that ferroelectric phase transition in STO under strain has characteristics of both displacive and order-disorder phase transition.
These simulations provide new insights into the behavior of the ferroelectric phase transition of STO. Ultimately, the machine learning method provides significant promise for the successful development of interatomic potentials for more complex perovskite materials.