A thermally driven needle-like (NL) to distorted perovskite (DP) phase transition in SrZrS3 was investigated by means of ab initio free energy calculations accelerated by machine learning. As a first step, a systematic screening of the methods to include long-range interactions in semilocal density functional theory Perdew–Burke–Ernzerhof calculations was performed. Out of the ten correction schemes tested, the Tkatchenko–Scheffler method with iterative Hirshfeld partitioning method was found to yield the best match between calculated and experimental lattice geometries, while predicting the correct order of stability of NL and DP phases at zero temperature. This method was then used in free energy calculations, performed using several approaches, so as to determine the effect of various anharmonicity contributions, such as the anisotropic thermal lattice expansion or the thermally induced internal structure changes, on the phase transition temperature (TNPDP). Accounting for the full anharmonicity by combining the NPT molecular dynamics data with thermodynamic integration with harmonic reference provided our best estimate of TNLDP = 867 K. Although this result is ∼150 K lower than the experimental value, it still provides an improvement by nearly 300 K compared to the previous theoretical report by Koocher et al. [Inorg. Chem. 62, 11134–11141 (2023)].

Chalcogenide perovskites (CPs), with composition ABX3 (A and B are metal cations and X = S, Se), have been recently identified as promising alternatives for halide perovskites in terms of chemical and thermal stability and, for a suitable choice of A and B, also a low toxicity.1–4 Several recently published experimental and theoretical studies have demonstrated the potential of these materials in a variety of applications by offering new insights into their syntheses, fundamental opto-electronic properties, and an understanding on the effect of various thermodynamic conditions.5–11 Similar to their halide counterparts, these materials also exist in several competing crystal phases differing in the arrangement of the BX6 octahedral framework.3,12 The two most prominent phases are the needle-like (NL) perovskite with NH4CdCl3-type structures (see Fig. 1) and the distorted perovskite (DP) phase with the GdFeO3 structure. Both phases form crystals with orthorhombic symmetry (space group Pnma), but they significantly differ in internal structures. In particular, the NL phase consists of mutually unbound one-dimensional chains of edge-sharing octahedra, while the DP phase comprises a three-dimensional network of corner-sharing octahedra.

FIG. 1.

Crystal structures of SrZrS3: NL phase projected onto the ac crystal plane (left) and DP phase projected onto the bc plane (right). Note that the 1D chain comprising the structure of NL is parallel to the crystal b axis. The blue, gray, and saffron colored spheres represent the Sr, Zr, and S atoms, respectively, while the black solid lines denote the supercells used in calculations consisting of 2 × 4 × 1 (NL) and 2 × 2 × 2 (DP) multiples of conventional orthorhombic cells (both with the space group Pnma).

FIG. 1.

Crystal structures of SrZrS3: NL phase projected onto the ac crystal plane (left) and DP phase projected onto the bc plane (right). Note that the 1D chain comprising the structure of NL is parallel to the crystal b axis. The blue, gray, and saffron colored spheres represent the Sr, Zr, and S atoms, respectively, while the black solid lines denote the supercells used in calculations consisting of 2 × 4 × 1 (NL) and 2 × 2 × 2 (DP) multiples of conventional orthorhombic cells (both with the space group Pnma).

Close modal

So far, most of the CP compositions, which have been experimentally synthesized, such as BaZrS3, CaZrS3, SrTiS3, BaHfS3, SrHfS3, and EuHfS3,5,13,14 have been found to preferably exist in the DP phase. However, according to the experimental report by Lee et al., SrZrS3 was found to occur in both of these aforementioned phases depending on the temperature (T) of synthesis.14 A first-order phase transformation is observed from the low T NL phase to the high T DP phase at 1253 K during heating and reverse transformation from DP to NL at 1023 K.14 Importantly, although NL is the thermodynamically stable phase at room T, it is possible to synthesize DP at the same conditions as a long-lived metastable phase.5,14 Another prominent phase transition that has been observed in SrZrS3 is the one where the DP phase further transforms quasi-continuously to the undistorted parent cubic phase (C) at ∼1200–1500 K.6,7 Such a continuous phase transition, which only causes deformation and local rearrangements within the cell without breaking and forming bonds, could be investigated using straightforward ab initio molecular dynamics (AIMD), as reported earlier.6,7 However, unlike in the DP → C transition, structural conversion from the NL to DP phase at high T can only occur through bond reorganization owing to the different octahedral connectivity of the phases (see Fig. 1). The latter implies that the phase transition involves high free energy barriers, making such a process a rare event15 not accessible by straightforward AIMD simulations with the computational power available today. Such phase transitions are therefore studied without explicitly considering the transition mechanism. Assuming thermodynamic equilibrium, the latter can be achieved by investigating the relative phase stability determined by suitably chosen thermodynamic potential, such as the Gibbs free energy (G) in the case of the NPT ensemble with a fixed number of particles (N), temperature, and pressure (P).

Various methods have been proposed to calculate G. The most popular and relatively inexpensive method is the quasi-harmonic approximation (QHA),16,17 which describes anharmonicity only partly, by including volume dependencies of vibrational frequencies. Although the QHA method has been successfully used to study thermal properties of various materials,18,19 this approach has proven to be inadequate in many cases, mainly at high T, where anharmonicity becomes important.20,21 Another way to treat anharmonicity is by including higher-order anharmonic terms in the Taylor expansion of the potential energy or non-perturbatively via a self-consistent approach. Methods such as temperature dependent effective potential (TDEP),22 self-consistent phonon theory (SCPH),23 stochastic self-consistent harmonic approximation (SSCHA),24 and the velocity auto-correlation function from molecular dynamics trajectories25 have been developed based on this approach.26–28 

In the context of free energy calculations, molecular dynamics (MD) methods are highly pertinent because they naturally involve all anharmonicity effects. One of the most popular MD-based methods to study phase transitions,29–32 without the need to define reversible transformation path in configuration space, is the thermodynamic integration with harmonic reference [Thermodynamic integration (TI)].33,34 Unfortunately, such calculations are severely limited due to their high computational cost when linked with ab initio level of theory. Recent advances in adaptive machine learning force fields31 (MLFFs) make these demanding simulations more accessible for routine use. An example of a successful combination of AIMD-based TI calculation with MLFF is the recent investigation of phase transitions in inorganic halide perovskites by Fransson et al.32 

In addition to selecting an appropriate method for calculating the free energy, another challenge is to choose a suitable functional that can correctly describe long-range dispersive interactions. Including these interactions is crucial for accurate calculation of the structure and energetics, especially in perovskites, which consist of highly polarizable atoms. This has been extensively demonstrated for various compositions of halide perovskites.35–37 Egger and Kronik,35 for instance, investigated the role of van der Waals (vdW) corrections for MAPbX3 (X = Cl, Br, I) using the PBE+TS dispersion scheme and inferred that such corrections cause a significant contraction of the unit cell compared to the one predicted by PBE, leading to better estimates of structural parameters. Similar observations for CP were also presented in our previous report,6 where correcting the results of the PBE functional for missing long-range dispersion interactions via the D3(0) method38 significantly improved structural descriptions at finite temperatures. However, as we shall show, employing the D3(0) corrections reverses the order of stability of the NL and DP phase, which is inconsistent with experimental observation. Hence, the choice of an appropriate dispersion correction method that, when linked with the given density functional approximation (here PBE), yields the best predictions for both structure and energies of well-known phases in CPs is still an open question requiring attention.

In this study, we systematically analyze the role of two important factors that contribute to accurate prediction of the phase transition temperature in SrZrS3, which are the role of long-range dispersion interactions and the importance of anharmonic effects in free energy calculations. This paper is organized as follows: the methodology used in this work is described in Sec. II, the results of screening of dispersion correction methods are presented in Sec. III A, schemes for including various types of anharmonic contributions to free energy are discussed in Sec. III B, our theoretical predictions of transition temperature are confronted with the experimental data in Sec. IV, and a summary of the most important results is given in Sec. V.

Ab initio simulations were employed within periodic density functional theory (DFT), as implemented in the Vasp software package.39–41 The Kohn–Sham equations were solved using the projected augmented wave (PAW) method of Blöchl.42 The exchange–correlation energy was treated with the Perdew–Burke–Ernzerhof (PBE) functional43 within the generalized gradient approximation. The kinetic energy cutoff to the plane-wave basis set was set to 400 eV, and the convergence threshold for the electronic energies was set to 10−6 eV. To account for long range van der Waals (vdW) interactions, various popular schemes of dispersion corrections were employed, namely, D2,44 D3 with the zero-damping function [D3(0)] and Becke–Johnson damping [D3(BJ)],38 and D445–47 of Grimme et al., Tkatchenko–Scheffler (TS) method,48 Tkatchenko–Scheffler method with iterative Hirshfeld49 partitioning (TS/HI),50,51 many-body dispersion energy52,53 with fractionally ionic model for polarizability54 method (MBD/FI) and its two-body variant (TS/FI),55 and density-dependent energy correction dDsC,56 as implemented in Vasp.50,51,55,57–59 Geometry optimizations were carried out using an external optimizer Gadget60,61 until all atomic forces were smaller than 0.005 eV/Å. The vibrational frequency calculations were performed numerically using central differences with the numerical step of 0.02 Å.

In the relaxations performed within the initial screening of dispersion correction methods (see Sec. III A), conventional unit cells containing 20 atoms and the k-points grids of 4 × 8 × 2 (NL) and 4 × 4 × 4 (DP) were employed. Consistent with our previous work,7 supercells consisting of 2 × 4 × 1 (NL) and 2 × 2 × 2 (DP) multiples of conventional unit cells containing 160 atoms each (see Fig. 1) were used in AIMD and all free energy calculations. In both cases, the Brillouin zone was sampled using a grid of 2 × 2 × 2 k-points, making the setting used in AIMD simulations consistent with that employed in relaxations with smaller structural models.

Temperature dependent properties were investigated using AIMD accelerated by machine learning force fields (MLFFs), as implemented in Vasp.62,63 These simulations were performed in the NPT ensemble at temperatures of 300, 600, 900, and 1200 K. The simulation T was kept constant by employing a Langevin thermostat64–66 with a friction coefficient of 2 ps−1 for all atoms. An external pressure of 0 GPa was maintained by using the Parrinello–Rahman barostat.67,68 The equations of motion of the ions were integrated using a velocity Verlet algorithm34 with a time step of 3 fs. A mass of 2 amu and a friction coefficient of 10 ps−1 were used for the lattice degrees of freedom.

The MLFF model for energy, forces, and stresses was trained using Bayesian regression.62,69 The training configurations were chosen automatically and on-the-fly during an AIMD run based on the similarity with other training configurations and quality of predictions of energies, forces, and stress tensor components.70 In our calculations, a separate training for each phase and T was carried out with the initial Bayesian error threshold of 0.01 eV/Å. The training typically consisted of 2 × 104 MD steps. The quality of the MLFF is evaluated in Sec. SIV of the supplementary material. Post-training, a production run of length 3 ns was performed in the NPT ensemble by employing the learned force fields. Observables such as the lattice parameters, volume, and energy were obtained as ensemble averages after discarding the data corresponding to equilibration, which was determined using the Mann–Kendall test71–73 for trend in mean and variance. MD simulations realized within the thermodynamic integration (see Sec. III B 3) were performed in the NVT ensemble using the MLFF trained in NPT MD at the same T. The thermostat setting and the integration step were the same as in NPT simulations. The length of the trajectories produced for different λ values (vide infra) was 500 ps.

In the following text, kB and are the Boltzmann and the reduced Planck constants, respectively, s represents the vector of atomic positions expressed in fractional coordinates, and h is the lattice matrix defined using the lattice vectors of the supercell (ai, i = 1, 2, 3) as h = {a1, a2, a3}. The atomic positions and cell geometries obtained in relaxations at a fixed supercell volume V = |det{h}| are indicated as sV and hV, respectively, whereby symbol V0(T) is reserved for the ground state cell volume at the temperature T. The term ⟨hP,T is used for the NPT ensemble average of h evaluated for specific values of P and T, and shP,T indicates relaxed atomic positions for the cell geometry given by ⟨hP,T. Finally, ωp(s, h) and Ael(s, h) are, respectively, the harmonic frequency of the phonon branch p and the electronic free energy computed for the structure given by s and h. Note that for a singlet state (all systems discussed in this work), Ael(s, h) corresponds to the electronic energy (E(s, h)) obtained in a QM calculation based on the Born–Oppenheimer approximation.

As discussed in Sec. I, long-range dispersion interactions are essential to accurately describe both the structure and energies of perovskite materials. Before starting time consuming free energy calculations, it is therefore useful to screen the dispersion correction methods for the best performance in static 0 K calculations. To proceed systematically, we tested all the methods presently available in Vasp that can be linked with the PBE functional (see Sec. II A). As a screening criterion, a lowest maximum error in the computed lattice geometry determined with an experimental reference14 is used. In addition, we also require that the chosen method predicts a lower potential energy for the relaxed NL structure, which is experimentally known to be stable at low values of T. This requirement is motivated by the fact that, as we shall show in Sec. III B, the thermal and entropic effects stabilize the DP structure relative to NL, i.e., only a method that predicts a lower energy for the relaxed NL can be consistent with the observed low T phase stability.

As a reference, we use sets of lattice constants a = 8.525 Å, b = 3.826 Å, c = 13.925 Å (NL), and a = 7.108 Å, b = 9.772 Å, and c = 6.741 Å (DP) determined experimentally14 at the room T. To facilitate a comparison with structural data obtained in structural relaxations, the experimental data were extrapolated to 0 K using the thermal expansion coefficients (TECs) determined in our previous AIMD study,6 which results in a = 8.509 Å, b = 3.806 Å, and c = 13.834 Å (NL) and a = 7.098 Å, b = 9.731 Å, and c = 6.704 Å (DP). We note that the result of extrapolation is almost independent of the dispersion correction used in TEC calculations.6 

Figure 2 shows the relative error in lattice constants with respect to experimental reference computed for various correction methods, and the corresponding numerical values of lattice constants are listed in Table S1 of the supplementary material. The largest relative errors obtained using the uncorrected PBE functional for the NL and DP phase are 1.6% and 1.4%, respectively. The correction methods D3(0), D3(BJ), D4, dDsC, TS/HI, and TS/FI systematically improve the quality of prediction for all three lattice parameters so that the absolute maximum errors do not exceed 0.7% [D3(0)], 1.0% [D3(BJ)], 0.5% (dDsC), 1.4% (TS/FI), and 0.4% (TS/HI and D4) for NL and 0.6% [D3(0)], 0.4% [D3(BJ) and dDsC], 0.5% (TS/FI), and 0.3% (TS/HI and D4) for DP. Another reasonably well performing method is D2, which improves all parameters except c of the NL phase, with the maximum relative errors being 1.4% (NL) and 0.7% (DP). The remaining two methods degrade the quality of predictions for one (TS) or both (MBD/FI) phases as compared to the uncorrected PBE results. As evident from the comparison of the results of the MBD/FI method with those obtained using its two-body variant TS/FI, a large part of failure of the former is due to the deficiency of the physical model behind the MBD method that results from inconsistencies in the model polarizability and screening equations (see the work of Gould et al.55 for detailed discussion of the problem). The computed energy differences between the NL and DP phases obtained using different methods are also shown in Fig. 2. Positive values of ΔENLDP = E(DP) − E(NL) are predicted by all methods except D3(0) and TS. Altogether, the correction methods D3(BJ), TS/FI, dDsC, D4, and TS/HI systematically improve upon the uncorrected PBE predictions for all lattice parameters of both phases of SrZrS3 considered here and, at the same time, predict correct stability order for the NL and DP phases. The performance of the latter three methods clearly stands out and is nearly identical. However, the TS/HI method predicts the largest potential energy difference between the relaxed DP and NL phases, which, as we shall discuss in Sec. IV, is important for minimizing error in predicted transition temperature with respect to the experimental reference. We therefore opted to use the latter method in the thermodynamic calculations presented in Sec. III B. We emphasize, however, that a good performance of the TS/HI and some other methods in prediction of structure and energetics of the SrZrS3 system does not necessarily imply a good performance for any other compound. Indeed, a careful benchmarking of available correction methods for the system under study is a good practice.

FIG. 2.

Relative errors (in %) with respect to reference experimental values14 extrapolated to 0 K for the lattice constants (a, b, and c) of (left) NL and (middle) DP phases of SrZrS3 and (right) zero temperature ground state energy difference between the NL and DP phases (ΔENLDP = ENLEDP) obtained from structural relaxations using different dispersion corrections. Note that a positive value of ΔENLDP implies that the NL phase is more stable than the DP at 0 K.

FIG. 2.

Relative errors (in %) with respect to reference experimental values14 extrapolated to 0 K for the lattice constants (a, b, and c) of (left) NL and (middle) DP phases of SrZrS3 and (right) zero temperature ground state energy difference between the NL and DP phases (ΔENLDP = ENLEDP) obtained from structural relaxations using different dispersion corrections. Note that a positive value of ΔENLDP implies that the NL phase is more stable than the DP at 0 K.

Close modal

Having established the superior performance of the TS/HI method compared to other dispersion correction methods at 0 K, we now employ this approach to access the relative stability of the NL and DP phases of SrZrS3 at thermodynamic conditions given by external pressure (P) and thermodynamic temperature (T). The relevant thermodynamic potential to be analyzed is the Gibbs free energy (G).74 In particular, the positive value of the difference between G computed for the NL and DP phases (GDPGNL) indicates that the NL phase is stable, the negative value corresponds to conditions where DP is stable, while the vanishing value is achieved when the two phases are at equilibrium, and the corresponding T at a fixed P corresponds to the temperature of phase transition (TNLDP).

In order to access G in computer simulations, it is usually more convenient to determine Helmholtz free energies (A) and employ the thermodynamic relation74 
G(T,P)=A(T,V(P))+PV(P),
(1)
where V(P) is the value of volume evaluated for a particular P. In this work, P = 0 is considered, simplifying Eq. (1) to identity G(T, P) = A(T, V(P)). To obtain an insight into different contributions to ΔGNLDP, we examined a hierarchy of methods for the calculation of A.

1. Harmonic approximation

Within the semi-classical harmonic approximation (HA), the Helmholtz free energy of a crystalline system can be written as
AHA(T)=Ael(sV0(T=0K),hV0(T=0K))kBTp=13N3lnkBTωp(sV0(T=0K),hV0(T=0K)),
(2)
where N is the number of atoms in the supercell that must be chosen so as to achieve convergence of AHA per formula unit (f.u.) and summation is over all modes with non-vanishing frequency. Importantly, the frequencies are independent of the cell geometry within this approximation. Consequently, the harmonic crystal does not expand and AHA only depends on T. We note that we chose to use semi-classical statistical mechanics expressions in this work for consistency with AIMD calculations employed at certain stages of our calculations (vide infra), which also generate semi-classical ensembles. However, as we show in Sec. SII of the supplementary material, replacing Eq. (2) by its quantum mechanical analogue has only minor effect on the quality of our results.

According to our HA calculations, the NL phase is stable up to TNLDP = 824 K (see Fig. 3), while the DP phase dominates at higher temperatures. Since the difference in the electronic contribution to enthalpic difference (HNLDP) remains fixed at ΔENLDP = 0.092 eV/f.u. and the difference in the vibrational contributions to HNLDP vanishes at all temperatures (by equipartition principle, classical vibrational enthalpy of harmonic system equals 3N kBT), the thermal stabilization of the DP phase is entirely entropic in nature and, qualitatively, originates in generally lower phonon frequencies of most of vibrational modes of DP.75 Within the HA, the term −TΔSNLDP decreases with T approximately linearly from −0.034 eV/f.u. at 300 K to −0.102 eV/f.u. at 900 K.

FIG. 3.

Free energy differences between the DP and NL phase (ΔGNLDP) determined for the temperature range 0–900 K using the harmonic approximation (HA), static quasi-harmonic approximation (S-QHA), quasi-harmonic approximation performed with structures from NPT AIMD (AIMD-QHA), and fully anharmonic thermodynamic integration (TI) performed at the PBE+TS/HI level. Transition temperatures are identified as the points for which the phase equilibrium condition ΔGNLDP = 0 is fulfilled.

FIG. 3.

Free energy differences between the DP and NL phase (ΔGNLDP) determined for the temperature range 0–900 K using the harmonic approximation (HA), static quasi-harmonic approximation (S-QHA), quasi-harmonic approximation performed with structures from NPT AIMD (AIMD-QHA), and fully anharmonic thermodynamic integration (TI) performed at the PBE+TS/HI level. Transition temperatures are identified as the points for which the phase equilibrium condition ΔGNLDP = 0 is fulfilled.

Close modal

2. Quasi-harmonic approximations

Within the quasi-harmonic approximation (QHA),16 a part of anharmonicity is taken into account via the implicit dependence of harmonic frequencies on the cell geometry given by h. In the simplest approach, which we denote here as the static QHA (S-QHA),18,19 a series of structures with s and h relaxed for a range of fixed cell volumes is considered, making the computed Helmholtz free energy,
AQHA(T,V)=Ael(sV,hV)kBTp=13N3lnkBTωp(sV,hV),
(3)
explicitly volume dependent. In our calculations, a series of 13 regularly spaced values of V were considered from the interval 0.94 V0V ≤ 1.12 V0. The corresponding atomic positions and lattice geometries were relaxed, vibrational analyses were performed, and AQHA(T, V) values were determined for all relevant temperatures. For each T, the resulting AQHA(T, V) vs V dependence was then fitted by the Murnaghan76 equation of state with the ground state free energy (AQHA(T, V0(T))), ground state volume (V0(T)), bulk modulus at zero pressure (B0(T)), and its pressure derivative [B0(T)] being the fitting parameters. Finally, the cell geometry corresponding to the cell volume of V0(T) was obtained via interpolation from dependencies of lattice constant on V.

The T dependencies of lattice parameters of the NL and DP phases computed using S-QHA are shown in Fig. 4. As evident, S-QHA predicts a monotonous and approximately linear increase of all lattice parameters with T, and qualitatively, the same behavior is observed for both phases. The monotonously growing deviation of the finite-T structures from the 0 K structure considered in the HA approach causes the computed GQHA to differ from GHA and the difference increases with T. In the case of the NL phase, the computed anharmonicity (GQHAGHA) is −0.004 eV/f.u. at 300 K and −0.041 eV/f.u. at 900 K, while for the DP structure, the values of −0.004 eV/f.u. (300 K) and −0.038 eV/f.u. (900 K) were obtained. The term −TΔSNLDP decreases from −0.034 to −0.094 eV/f.u. over the interval 300–900 K, i.e., the entropic stabilization of DP is slightly weaker compared to the HA prediction. This is, however, partly compensated by the electronic contribution to ΔHNLDP, which, due to variable cell shape, decreases at 900 K from 0.092 eV/f.u. (HA) to 0.088 eV/f.u. (S-QHA). Thus, the part of anharmonicity that is taken into account within the S-QHA has a rather modest impact on the value of TNLDP, which is predicted to be 851 K, i.e., 27 K higher compared to the HA prediction (see Fig. 3).

FIG. 4.

Comparison of the thermal variation of scaled lattice parameters for (a) the NL (a′ = a/2, b′ = b, c′ = c/4) and (b) DP (a′ = a/2, b′ = b/2, c′ = c/2) phases of SrZrS3 obtained using the S-QHA (denoted by dashed lines) and NPT AIMD (denoted by solid lines) approaches at the PBE+TS/HI level.

FIG. 4.

Comparison of the thermal variation of scaled lattice parameters for (a) the NL (a′ = a/2, b′ = b, c′ = c/4) and (b) DP (a′ = a/2, b′ = b/2, c′ = c/2) phases of SrZrS3 obtained using the S-QHA (denoted by dashed lines) and NPT AIMD (denoted by solid lines) approaches at the PBE+TS/HI level.

Close modal
Clearly, the S-QHA approach assumes that all the thermally induced cell shape changes are fully determined by the dependence of h on V at zero T. This serious restriction, causing that the S-QHA is inherently unable to describe T-dependent anisotropy in thermal expansion, can be lifted by combining QHA with AIMD in the NPT ensemble, which fully accounts for anharmonicity effects on the structure. Within this approach, which we denote as AIMD-QHA, the Helmholtz free energy is determined using
AQHA(T,hP,T)=Ael(shP,T,hP,T)kBTp=13N3lnkBTωp(shP,T,hP,T),
(4)
where the ensemble averages ⟨hP,T are obtained from the MLFF-accelerated NPT AIMD described in Sec. II B.

Thermal variation of the lattice constants computed using AIMD-QHA is compared with the S-QHA results in Fig. 4. In the case of the NL phase, both methods predict qualitatively similar trends with a nearly linear increase of all lattice constants. The two methods are also in semi-quantitative agreement, whereby the AIMD-QHA predicts a slightly lower slope for the increase of a and a higher slope for the increase of b with T. At 900 K, where the deviation between structural predictions of both methods is the largest, the AIMD-QHA method yields the lattice parameters a, b, and c that differ by −0.83%, 0.44%, and 0.16% from those obtained by S-QHA, while the difference in V is −0.26%. Consequently, GQHA computed using both methods are very similar (see Table S2 of the supplementary material). A more significant difference in structure is observed for the DP phase, where, at 900 K, AIMD-QHA predicts deviations from the quasi-linear increase of a, b, c, and V obtained using S-QHA by −1.07%, 0.13%, 0.41%, and −0.58%. As discussed in detail in our previous works,6,7 the structural changes of the lattice geometry of DP with increasing T are due to a quasi-continuous transformation to the parent cubic phase (C). This is exactly a type of thermal effect that cannot be described by the S-QHA approach and leads to more significant differences in GQHA compared to the NL phase (see Table S2 of the supplementary material). Since, compared to the S-QHA, the AIMD-QHA predicts lower GQHAGHA for DP [−0.040 eV/f.u. (AIMD-QHA) vs −0.038 eV/f.u. (S-QHA)] and almost the same values for NL (−0.041 eV/f.u.), ΔGNLDP increases. Thus, somewhat surprisingly, the differences ΔGNLDP computed via AIMD-QHA are virtually identical to the HA results over the entire range of temperatures (see Table S2 of the supplementary material), which is a consequence of compensation of anharmonic contributions of the ΔHNLDP and −TΔSNLDP terms. At 900 K, for instance, the anharmonicity contributes −0.007 eV/f.u. to ΔHNLDP and 0.008 eV/f.u. to −TΔSNLDP, leading to anharmonic ΔGNLDP = 0.001 eV/f.u. Consequently, the AIMD-QHA prediction for TNLDP (826 K) is very similar to that obtained from HA (824 K); see Fig. 3.

3. Thermodynamic integration (TI)

The QHA approaches are intrinsically unable to account for anharmonicity due to thermally induced internal structure changes. In particular, the quasi-continuous transformation of DP to C (see Sec. SIII of the supplementary material) is associated with the decrease of transversal displacement of S atoms from its ideal position in the middle of the Zr–Zr link (involving Zr–S bonds from pairs of ZrS6 octahedra sharing a common atom S) from ∼0.7 to 0.0 Å.6 However, the AIMD-QHA value of free energy can be corrected for the missing anharmonicity contribution via thermodynamic integration (TI) with harmonic reference.33,77,78 In this method, Helmholtz free energy of a fully interacting system (1) driven by Hamiltonian H1 is expressed as follows:
A1=A0+ΔA01,
(5)
where A0 = AQHA(T, ⟨hP,T), determined using Eq. (4), is the free energy of reference harmonic state (0) driven by the Hamiltonian H0 and ΔA0→1 is the free energy difference between systems 0 and 1 calculated as follows:
ΔA01=01E1(λ)E0(λ)λdλ.
(6)
In Eq. (6), the term ⟨⋯⟩λ represents the average of enclosed quantity over the NVT ensemble generated by the Hamiltonian
Hλ=λH1+(1λ)H0,
(7)
and λ is the coupling parameter.

We note that state 0 can be chosen arbitrarily and does not need to correspond to the best harmonic approximation of the fully interacting system. Following Ref. 78, we defined state 0 in terms of rotationally and translationally invariant internal coordinates (here, a redundant set of interatomic distances was used) as follows. First, the atomic positions were relaxed while keeping the cell geometry fixed at ⟨hP,T obtained from NPT AIMD. Subsequently, the Hesse matrix was obtained via a finite difference method, whereby the eigenvalues of this matrix lower than the limit 1 eV Å−1 were increased to that limit while keeping the eigenvectors unaltered. As discussed in Ref. 78, this procedure avoids numerical problems due to nonphysical configurations generated by Hλ for low λ values. All these steps were performed at the DFT level using the setting described in Sec. II A. Finally, the Hesse matrix was expressed in internal coordinates, which were employed in the TI calculations; see Ref. 78 for details.

Integration in Eq. (6) was realized numerically by the method proposed in Ref. 79 that is based on the Bennett acceptance ratio approach.80 A regular grid of 5 λ points including the two integration limits was used for that purpose. The calculations of E1(λ)E0(λ)λ were performed in the NVT ensemble using the MLFF obtained in previous NPT MD runs at the same T (see Sec. II B), whereby the lengths of production runs were 450 ps for each value of λ.

Since the MLFF represents an approximation to true DFT potential energy landscape, our MLFF-based TI calculations are biased to some extent. To correct for the error introduced by MLFF, free energy perturbation theory81–83 (FEPT) was used, whereby the corrected Helmholtz free energy (A1,DFT) was obtained from the MLFF value (A1,MLFF) as follows:
A1,DFT=A1,MLFFkBTlnexp[ΔE(s,h)/kBT]MLFF.
(8)
In Eq. (8), the term ⟨⋯⟩MLFF is the NVT ensemble average generated by the Hamiltonian corresponding to MLFF and ΔE(s, h) = EDFT(s, h) − EMLFF(s, h) is the potential energy difference between DFT and MLFF computed for the configuration with atomic positions s and lattice geometry h. In the calculations of exp[ΔE(s,h)/kBT]MLFF, a second-order cumulant expansion82 employing 250 configurations randomly chosen from the trajectory generated in TI with λ = 1 was used.

Taking full anharmonicity into account leads to stabilization of both the NL and DP phases. Nevertheless, the term GTIGHA computed for NL is greater in absolute value [−0.048 eV/f.u. (NL) vs −0.043 eV/f.u. (DP)]. Thus, anharmonicity stabilizes the NL relative to the DP phase. This can be intuitively understood because the former structure consists of one-dimensional ZrS3 chains that are not connected by covalent bonds and are therefore less accurately described by the harmonic model than the fully connected 3D framework of the DP structure (see Fig. 1). This additional stabilization leads to an increase in the predicted transition temperature from 826 K obtained in AIMD-QHA calculations to 867 K (see Fig. 3). Altogether, our best estimate of the anharmonicity effect on TNLDP, obtained by comparing the HA and TI values, is 43 K.

According to the experiment,14 both the NL and the DP phases can exist at ambient conditions, but the NL → DP phase transition occurs when T of the former is increased to 1250 K. This process can be reversed, whereby the DP phase is fully transformed into NL when heated to 1020 K.14 The presence of hysteresis in the experimental thermochemical data indicates that the observed phase transitions followed irreversible transformation paths, and hence, the reported transition temperatures should be interpreted only as upper and lower bounds for TNLDP determined for ideal reversible process considered in our work. Our best estimate of TNLDP = 867 K obtained in fully anharmonic treatment of atomic and lattice degrees of freedom therefore underestimates experimental transition T by at least 150 K. Although this error might seem large, our result still represents a significant improvement compared to the previous theoretical report of Koocher et al.75 based on self-consistent QHA calculations performed at the PBEsol84 level, which underestimated TNLDP by at least 430 K. We note, however, that the seemingly large error in theoretical predictions is mainly due to the very similar free energies of the NL and DP phases (see Table S2 of the supplementary material) and only small differences in the slopes of their T dependencies, which makes theoretical results very sensitive on the choice of the simulation methodology. In our work, we have shown that various types of anharmonic effects, such as thermal expansion of cell volume or thermally induced cell shape and internal structure changes, have only a relatively modest effect (within 50 K) on predicted TNLDP. Thus, the choice of the electronic structure method, directly affecting both the electronic energy and harmonic frequencies, turns out to be the most important factor in accurate calculations. To this end, we have shown that just the dispersion correction linked with the same functional (PBE) leads to variations in ΔENLDP between −0.121 (PBE+TS) and 0.099 eV/f.u. (PBE+D2), which, assuming that all the thermal and entropic effects are similar to those explicitly explored in this work using the PBE+TS/HI method, would translate into TNLDP between 0 and 900 K. Using the same assumption, we estimate that an ideal electronic structure method should yield ΔENLDP between 0.110 and 0.150 eV/f.u., which might serve, together with a good accuracy in structural predictions, as a criterion in the future screening studies involving multiple functionals and/or post-DFT treatments. The fact that the lower bound of the ideal ΔENLDP is just 0.017 eV/f.u. higher than the value of 0.093 meV/f.u. obtained by the PBE+TS/HI identified in this work as the best for describing the structure of both phases, illustrates the sensitivity of the problem.

A systematic investigation on the thermally driven NL to DP phase transition in SrZrS3 was performed by means of free energy calculations at various levels of complexity. The quality of various dispersion correction schemes for the PBE functional [D2, D3(0), D3(BJ), D4, TS, TS/HI, TS/FI, MBD-FI, and dDsC] was evaluated to determine their ability to predict the ground state geometries and potential energy differences of the NL and DP phases in agreement with available experimental data. Based on the results of our screening, the TS/HI correction was identified as the best, predicting correct zero-T relative stability of DP and NL phases and the lattice geometries with error within 0.4% for both phases. PBE+TS/HI was subsequently used in thermodynamic calculations.

Free energy calculations were performed to explore the relative stability of the NL and DP phases of SrZrS3 at various temperatures and to determine the transition T. In order to investigate the effect of different types of anharmonicity, a hierarchy of methods was used in the calculations. A simple harmonic approximation that ignores the thermal expansion of the lattice provided us with a reference to determine the magnitude of the anharmonicity. The value of TNLDP obtained using this approach is 824 K. The static quasi-harmonic approximation takes into account a part of anharmonicity via implicit dependence of harmonic frequencies on the cell volume determined for a series of relaxed structures. Applying this approach led to an increase of TNLDP to 851 K.

The assumption of the S-QHA that the lattice geometry changes are fully fixed by volume can be lifted by applying the quasi harmonic approximation for relaxed structures with lattice geometries obtained as ensemble averages from NPT AIMD. This part of anharmonicity has a more significant impact on the structure of the DP phase, which is primarily due to the tendency of DP to transform into parent cubic phase when T is increased.6,7 Hence, unsurprisingly, this effect decreased free energy of DP more than that of NL, shifting TNLDP to lower value (826 K) compared to the S-QHA result.

Since both the S-QHA and AIMD-QHA approach calculations are performed for structures with relaxed atomic positions, they inherently fail to account for the anharmonicity associated with thermally induced internal structure changes. Starting from the AIMD-QHA results, this additional contribution can be taken into account using thermodynamic integration with a harmonic reference. Compared to AIMD-QHA, TI leads to stronger stabilization of the NL structure, which is related to the fact that it consists of 1D chains of ZrS3 that are not bound by covalent bonds and are therefore less accurately described by the harmonic model than the fully connected 3D framework of the DP phase. Consequently, the transition temperature calculated at the TI level is shifted to a higher value of 867 K as compared to the AIMD-QHA prediction.

In comparison with the experiment,14 our best theoretical prediction based on TI is at least 150 K lower, but still represents a significant improvement over the previous theoretical report of Koocher et al.75 that underestimated TNLDP by at least 430 K. The main reason that complicates accurate theoretical predictions is the similarity in stability of both phases and in the slope of corresponding G vs T curves. A comparison of our TI and HA result suggests that the anharmonicity effect on TNLDP is rather modest (43 K), and therefore, the choice of the electronic structure method is a key factor that needs to be improved in future work. Based on our results, we estimate that an ideal electronic structure method should yield ΔENLDP between 0.110 and 0.150 eV/f.u. This proposed range of values can serve as one of the criteria in the future screening studies involving multiple functionals and/or post-DFT treatments.

See the supplementary material for additional information on the effects of dispersion corrections on structure and energies of the phases, comparison between semi-classical and quantum free energies as discussed in the main text, comparison of thermal effects on structures predicted by PBE, D3(0), and TS/HI methods, and a description of the quality of machine learned force fields used in the simulations.

This work was supported by the European Union’s Horizon 2020 Research and Innovation Program under Grant Agreement No. 810701 and by the Slovak Research and Development Agency under Grant Agreement No. APVV-19-0410. N.J. acknowledges partial support from Project USCCCORD (ZoNFP: NFP313020BUZ3), cofinanced by the European Regional Development Fund within the Operational Programme Integrated Infrastructure. This research was realized using the computational resources procured in the national project National competence centre for high performance computing within the Operational Programme Integrated infrastructure (Project Code: 311070AKF2).

The authors have no conflicts to disclose.

Namrata Jaykhedkar: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Visualization (equal); Writing – original draft (equal). Roman Bystrický: Writing – review & editing (equal). Milan Sýkora: Funding acquisition (lead); Project administration (equal); Writing – review & editing (equal). Tomáš Bučko: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Writing – original draft (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
S.
Niu
,
H.
Huyan
,
Y.
Liu
,
M.
Yeung
,
K.
Ye
,
L.
Blankemeier
,
T.
Orvis
,
D.
Sarkar
,
D. J.
Singh
,
R.
Kapadia
, and
J.
Ravichandran
, “
Bandgap control via structural and chemical tuning of transition metal perovskite chalcogenides
,”
Adv. Mater.
29
,
1604733
(
2017
).
2.
A.
Swarnkar
,
W. J.
Mir
,
R.
Chakraborty
,
M.
Jagadeeswararao
,
T.
Sheikh
, and
A.
Nag
, “
Are chalcogenide perovskites an emerging class of semiconductors for optoelectronic properties and solar cell?
,”
Chem. Mater.
31
,
565
575
(
2019
).
3.
K. V.
Sopiha
,
C.
Comparotto
,
J. A.
Márquez
, and
J. J. S.
Scragg
, “
Chalcogenide perovskites: Tantalizing prospects, challenging materials
,”
Adv. Opt. Mater.
10
,
2101704
(
2022
).
4.
M.-G.
Ju
,
J.
Dai
,
L.
Ma
, and
X. C.
Zeng
, “
Perovskite chalcogenides with optimal bandgap and desired optical absorption for photovoltaic devices
,”
Adv. Energy Mater.
7
,
1700216
(
2017
).
5.
R.
Bystrický
,
S. K.
Tiwari
,
P.
Hutár
,
L.
Vančo
, and
M.
Sýkora
, “
Synthesis of sulfide perovskites by sulfurization with boron sulfides
,”
Inorg. Chem.
61
,
18823
18827
(
2022
).
6.
N.
Jaykhedkar
,
R.
Bystrický
,
M.
Sýkora
, and
T.
Bučko
, “
How the temperature and composition govern the structure and band gap of Zr-based chalcogenide perovskites: Insights from ML accelerated AIMD
,”
Inorg. Chem.
62
,
12480
12492
(
2023
).
7.
N.
Jaykhedkar
,
R.
Bystrický
,
M.
Sýkora
, and
T.
Bučko
, “
Understanding the structure-band gap relationship in SrZrS3 at elevated temperatures: A detailed NPT MD study
,”
J. Mater. Chem. C
10
,
12032
12042
(
2022
).
8.
Y.-Y.
Sun
,
M. L.
Agiorgousis
,
P.
Zhang
, and
S.
Zhang
, “
Chalcogenide perovskites for photovoltaics
,”
Nano Lett.
15
,
581
585
(
2015
).
9.
Y.
Nishigaki
,
T.
Nagai
,
M.
Nishiwaki
,
T.
Aizawa
,
M.
Kozawa
,
K.
Hanzawa
,
Y.
Kato
,
H.
Sai
,
H.
Hiramatsu
,
H.
Hosono
, and
H.
Fujiwara
, “
Extraordinary strong band-edge absorption in distorted chalcogenide perovskites
,”
Sol. RRL
4
,
1900555
(
2020
).
10.
W.
Meng
,
B.
Saparov
,
F.
Hong
,
J.
Wang
,
D. B.
Mitzi
, and
Y.
Yan
, “
Alloying and defect control within chalcogenide perovskites for optimized photovoltaic application
,”
Chem. Mater.
28
,
821
829
(
2016
).
11.
S.
Perera
,
H.
Hui
,
C.
Zhao
,
H.
Xue
,
F.
Sun
,
C.
Deng
,
N.
Gross
,
C.
Milleville
,
X.
Xu
,
D. F.
Watson
,
B.
Weinstein
,
Y.-Y.
Sun
,
S.
Zhang
, and
H.
Zeng
, “
Chalcogenide perovskites—An emerging class of ionic semiconductors
,”
Nano Energy
22
,
129
135
(
2016
).
12.
J. A.
Brehm
,
J. W.
Bennett
,
M. R.
Schoenberg
,
I.
Grinberg
, and
A. M.
Rappe
, “
The structural diversity of ABS3 compounds with d0 electronic configuration for the b-cation
,”
J. Chem. Phys.
140
,
224703
(
2014
).
13.
R.
Lelieveld
and
D. J. W.
IJdo
, “
Sulphides with the GdFeO3 structure
,”
Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem.
36
,
2223
2226
(
1980
).
14.
C.-S.
Lee
,
K.
Kleinke
, and
H.
Kleinke
, “
Synthesis, structure, and electronic and physical properties of the two SrZrS3 modifications
,”
Solid State Sci.
7
,
1049
1054
(
2005
).
15.
M.
Tuckerman
and
B.
Smit
, “
Free energy calculations
,” in
Statistical Mechanics
(
Oxford University Press
,
2010
), pp.
155
158
.
16.
S.
Baroni
,
P.
Giannozzi
, and
E.
Isaev
, “
Density-functional perturbation theory for quasi-harmonic calculations
,” in
Theoretical and Computational Methods in Mineral Physics: Geophysical Applications
(
Walter de Gruyter GmbH, Berlin/Munich/Boston
,
2010
), pp.
39
57
.
17.
A.
Togo
and
I.
Tanaka
, “
First principles phonon calculations in materials science
,”
Scr. Mater.
108
,
1
5
(
2015
).
18.
P.
Nath
,
J. J.
Plata
,
D.
Usanmaz
,
R.
Al Rahal Al Orabi
,
M.
Fornari
,
M. B.
Nardelli
,
C.
Toher
, and
S.
Curtarolo
, “
High-throughput prediction of finite-temperature properties using the quasi-harmonic approximation
,”
Comput. Mater. Sci.
125
,
82
91
(
2016
).
19.
S.
Schryver
and
A.
Lamichhane
, “
Temperature-driven structural phase transitions in CsPbBr3
,”
Solid State Commun.
371
,
115237
(
2023
).
20.
J.
Li
and
P.
Rinke
, “
Atomic structure of metal-halide perovskites from first principles: The chicken-and-egg paradox of the organic-inorganic interaction
,”
Phys. Rev. B
94
,
045201
(
2016
).
21.
P.
Lazar
,
T.
Bučko
, and
J.
Hafner
, “
Negative thermal expansion of ScF3: Insights from density-functional molecular dynamics in the isothermal-isobaric ensemble
,”
Phys. Rev. B
92
,
224302
(
2015
).
22.
O.
Hellman
,
P.
Steneteg
,
I. A.
Abrikosov
, and
S. I.
Simak
, “
Temperature dependent effective potential method for accurate free energy calculations of solids
,”
Phys. Rev. B
87
,
104111
(
2013
).
23.
K.
Esfarjani
and
Y.
Liang
, “
Thermodynamics of anharmonic lattices from first principles
,” in
Nanoscale Energy Transport
(
IOP Publishing
,
2020
), pp.
7–1
7–35
.
24.
L.
Monacelli
,
R.
Bianco
,
M.
Cherubini
,
M.
Calandra
,
I.
Errea
, and
F.
Mauri
, “
The stochastic self-consistent harmonic approximation: Calculating vibrational properties of materials with full quantum and anharmonic effects
,”
J. Phys.: Condens. Matter
33
,
363001
(
2021
).
25.
A.
Carreras
,
A.
Togo
, and
I.
Tanaka
, “
Dynaphopy: A code for extracting phonon quasiparticles from molecular dynamics simulations
,”
Comput. Phys. Commun.
221
,
221
234
(
2017
).
26.
K.
Tolborg
,
J.
Klarbring
,
A. M.
Ganose
, and
A.
Walsh
, “
Free energy predictions for crystal stability and synthesisability
,”
Digital Discovery
1
,
586
595
(
2022
).
27.
I.
Errea
,
M.
Calandra
, and
F.
Mauri
, “
Anharmonic free energies and phonon dispersions from the stochastic self-consistent harmonic approximation: Application to platinum and palladium hydrides
,”
Phys. Rev. B
89
,
064302
(
2014
).
28.
L.
Monacelli
and
N.
Marzari
, “
First-principles thermodynamics of CsSnI3
,”
Chem. Mater.
35
,
1702
1709
(
2023
).
29.
S. G.
Moustafa
,
A. J.
Schultz
, and
D. A.
Kofke
, “
Very fast averaging of thermal properties of crystals by molecular simulation
,”
Phys. Rev. E
92
,
043303
(
2015
).
30.
S. G.
Moustafa
,
A. J.
Schultz
, and
D. A.
Kofke
, “
Harmonically assisted methods for computing the free energy of classical crystals by molecular simulation: A comparative study
,”
J. Chem. Theory Comput.
13
,
825
834
(
2017
).
31.
R.
Jinnouchi
,
F.
Karsai
, and
G.
Kresse
, “
Making free-energy calculations routine: Combining first principles with machine learning
,”
Phys. Rev. B
101
,
060201
(
2020
).
32.
E.
Fransson
,
J.
Wiktor
, and
P.
Erhart
, “
Phase transitions in inorganic halide perovskites from machine-learned potentials
,”
J. Phys. Chem. C
127
,
13773
13781
(
2023
).
33.
J. G.
Kirkwood
, “
Statistical mechanics of fluid mixtures
,”
J. Chem. Phys.
3
,
300
313
(
2004
).
34.
D.
Frenkel
and
B.
Smit
, “
Free energies and phase equilibria
,” in
Understanding Molecular Simulation
, 2nd ed. (
Academic Press
,
San Diego
,
2002
).
35.
D. A.
Egger
and
L.
Kronik
, “
Role of dispersive interactions in determining structural properties of organic–inorganic halide perovskites: Insights from first-principles calculations
,”
J. Phys. Chem. Lett.
5
,
2728
2733
(
2014
).
36.
H.
Beck
,
C.
Gehrmann
, and
D. A.
Egger
, “
Structure and binding in halide perovskites: Analysis of static and dynamic effects from dispersion-corrected density functional theory
,”
APL Mater.
7
,
021108
(
2019
).
37.
N.
Pandech
,
T.
Kongnok
,
N.
Palakawong
,
S.
Limpijumnong
,
W. R. L.
Lambrecht
, and
S.
Jungthawan
, “
Effects of the van der Waals interactions on structural and electronic properties of Ch3NH3(Pb,Sn)(I,Br,Cl)3 halide perovskites
,”
ACS Omega
5
,
25723
25732
(
2020
).
38.
S.
Grimme
,
S.
Ehrlich
, and
L.
Goerigk
, “
Effect of the damping function in dispersion corrected density functional theory
,”
J. Comput. Chem.
32
,
1456
1465
(
2011
).
39.
G.
Kresse
and
J.
Hafner
, “
Ab initio molecular dynamics for open-shell transition metals
,”
Phys. Rev. B
48
,
13115
13118
(
1993
).
40.
G.
Kresse
and
J.
Hafner
, “
Norm-conserving and ultrasoft pseudopotentials for first-row and transition elements
,”
J. Phys.: Condens. Matter
6
,
8245
8257
(
1994
).
41.
G.
Kresse
and
J.
Hafner
, “
Ab initio molecular dynamics for liquid metals
,”
Phys. Rev. B
47
,
558
561
(
1993
).
42.
P.
Blöchl
, “
Projector augmented-wave method
,”
Phys. Rev. B
50
,
17953
(
1994
).
43.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
3868
(
1996
).
44.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
, “
A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu
,”
J. Chem. Phys.
132
,
154104
(
2010
).
45.
E.
Caldeweyher
,
C.
Bannwarth
, and
S.
Grimme
, “
Extension of the D3 dispersion coefficient model
,”
J. Chem. Phys.
147
,
034112
(
2017
).
46.
E.
Caldeweyher
,
S.
Ehlert
,
A.
Hansen
,
H.
Neugebauer
,
S.
Spicher
,
C.
Bannwarth
, and
S.
Grimme
, “
A generally applicable atomic-charge dependent London dispersion correction
,”
J. Chem. Phys.
150
,
154122
(
2019
).
47.
E.
Caldeweyher
,
J.-M.
Mewes
,
S.
Ehlert
, and
S.
Grimme
, “
Extension and evaluation of the D4 London-dispersion model for periodic systems
,”
Phys. Chem. Chem. Phys.
22
,
8499
8512
(
2020
).
48.
A.
Tkatchenko
and
M.
Scheffler
, “
Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data
,”
Phys. Rev. Lett.
102
,
073005
(
2009
).
49.
P.
Bultinck
,
C.
Van Alsenoy
,
P. W.
Ayers
, and
R.
Carbó-Dorca
, “
Critical analysis and extension of the Hirshfeld atoms in molecules
,”
J. Chem. Phys.
126
,
144111
(
2007
).
50.
T.
Bučko
,
S.
Lebègue
,
J.
Hafner
, and
J. G.
Ángyán
, “
Improved density dependent correction for the description of London dispersion forces
,”
J. Chem. Theory Comput.
9
,
4293
4299
(
2013
).
51.
T.
Bučko
,
S.
Lebègue
,
J. G.
Ángyán
, and
J.
Hafner
, “
Extending the applicability of the Tkatchenko-Scheffler dispersion correction via iterative Hirshfeld partitioning
,”
J. Chem. Phys.
141
,
034114
(
2014
).
52.
A.
Tkatchenko
,
R. A.
DiStasio
,
R.
Car
, and
M.
Scheffler
, “
Accurate and efficient method for many-body van der Waals interactions
,”
Phys. Rev. Lett.
108
,
236402
(
2012
).
53.
A.
Ambrosetti
,
A. M.
Reilly
,
R. A.
DiStasio
,
A.
Tkatchenko
, and
A.
Tkatchenko
, “
Long-range correlation energy calculated from coupled atomic response functions
,”
J. Chem. Phys.
140
,
18A508
(
2014
).
54.
T.
Gould
and
T.
Bučko
, “
C6 coefficients and dipole polarizabilities for all atoms and many ions in rows 1–6 of the periodic table
,”
J. Chem. Theory Comput.
12
,
3603
3613
(
2016
).
55.
T.
Gould
,
S.
Lebègue
,
J. G.
Ángyán
, and
T.
Bučko
, “
A fractionally ionic approach to polarizability and van der Waals many-body dispersion calculations
,”
J. Chem. Theory Comput.
12
,
5920
5930
(
2016
).
56.
S. N.
Steinmann
and
C.
Corminboeuf
, “
A generalized-gradient approximation exchange hole model for dispersion coefficients
,”
J. Chem. Phys.
134
,
044117
(
2011
).
57.
T.
Bučko
,
J.
Hafner
,
S.
Lebègue
, and
J. G.
Ángyán
, “
Improved description of the structure of molecular and layered crystals: Ab initio DFT calculations with van der Waals corrections
,”
J. Phys. Chem. A
114
,
11814
11824
(
2010
).
58.
T.
Bučko
,
S.
Lebègue
,
J.
Hafner
, and
J. G.
Ángyán
, “
Tkatchenko-Scheffler van der Waals correction method with and without self-consistent screening applied to solids
,”
Phys. Rev. B
87
,
064110
(
2013
).
59.
T.
Bučko
,
S.
Lebègue
,
T.
Gould
, and
J. G.
Ángyán
, “
Many-body dispersion corrections for periodic systems: An efficient reciprocal space implementation
,”
J. Phys.: Condens. Matter
28
,
045201
(
2016
).
60.
T.
Bučko
,
J.
Hafner
, and
J. G.
Angyan
, “
Geometry optimization of periodic systems using internal coordinates
,”
J. Chem. Phys.
122
,
124508
(
2005
).
61.
T.
Bučko
, “
Transition state optimization of periodic systems using delocalized internal coordinates
,”
Theor. Chem. Acc.
137
,
164
(
2018
).
62.
R.
Jinnouchi
,
F.
Karsai
, and
G.
Kresse
, “
On-the-fly machine learning force field generation: Application to melting points
,”
Phys. Rev. B
100
,
014105
(
2019
).
63.
R.
Jinnouchi
,
F.
Karsai
,
C.
Verdi
,
R.
Asahi
, and
G.
Kresse
, “
Descriptors representing two- and three-body atomic distributions and their effects on the accuracy of machine-learned inter-atomic potentials
,”
J. Chem. Phys.
152
,
234102
(
2020
).
64.
D. J.
Evans
, “
Computer ‘experiment’ for nonlinear thermodynamics of Couette flow
,”
J. Chem. Phys.
78
,
3297
3302
(
1983
).
65.
W. G.
Hoover
,
A. J. C.
Ladd
, and
B.
Moran
, “
High-strain-rate plastic flow studied via nonequilibrium molecular dynamics
,”
Phys. Rev. Lett.
48
,
1818
1820
(
1982
).
66.
D. S.
Lemons
and
A.
Gythiel
, “
Paul Langevin’s 1908 paper ‘On the theory of Brownian motion’ [‘sur la théorie du mouvement Brownien,’ C. R. Acad. Sci. (Paris) 146, 530–533 (1908)]
,”
Am. J. Phys.
65
,
1079
1081
(
1997
).
67.
M.
Parrinello
and
A.
Rahman
, “
Crystal structure and pair potentials: A molecular-dynamics study
,”
Phys. Rev. Lett.
45
,
1196
1199
(
1980
).
68.
M.
Parrinello
and
A.
Rahman
, “
Polymorphic transitions in single crystals: A new molecular dynamics method
,”
J. Appl. Phys.
52
,
7182
7190
(
1981
).
69.
C. M.
Bishop
,
Pattern Recognition and Machine Learning
(
Springer
,
2006
).
70.
R.
Jinnouchi
,
J.
Lahnsteiner
,
F.
Karsai
,
G.
Kresse
, and
M.
Bokdam
, “
Phase transitions of hybrid perovskites simulated by machine-learning force fields trained on the fly with Bayesian inference
,”
Phys. Rev. Lett.
122
,
225701
(
2019
).
71.
M.
Kendall
and
A.
Stuart
,
The Advanced Theory of Statistics: Design, Analysis and Time-Series
(
Charles Griffin
,
1966
), Vol.
3
.
72.
M.
Kendall
and
J. D.
Gibbons
,
Rank Correlation Methods
, 5th ed. (
A Charles Griffin Title
,
1990
).
73.
S. K.
Schiferl
and
D. C.
Wallace
, “
Statistical errors in molecular dynamics averages
,”
J. Chem. Phys.
83
,
5203
5209
(
1985
).
74.
M.
Tuckerman
and
B.
Smit
, “
The isobaric ensembles
,” in
Statistical Mechanics
(
Oxford University Press
,
2010
), pp.
155
158
.
75.
N. Z.
Koocher
and
J. M.
Rondinelli
, “
Effect of octahedral connectivity on the negative thermal expansion of SrZrS3
,”
Inorg. Chem.
62
,
11134
11141
(
2023
).
76.
F. D.
Murnaghan
, “
The compressibility of media under extreme pressures
,”
Proc. Natl. Acad. Sci. U. S. A.
30
,
244
247
(
1944
).
77.
M.
Jorge
,
N.
Garrido
,
A.
Queimada
,
I.
Economou
, and
E.
Macedo
, “
Effect of the integration method on the accuracy and computational efficiency of free energy calculations using thermodynamic integration
,”
J. Chem. Theory Comput.
6
,
1018
1027
(
2010
).
78.
J.
Amsler
,
P. N.
Plessow
,
F.
Studt
, and
T.
Bučko
, “
Anharmonic correction to adsorption free energy from DFT-based MD using thermodynamic integration
,”
J. Chem. Theory Comput.
17
,
1155
1169
(
2021
).
79.
J.
Amsler
,
P. N.
Plessow
,
F.
Studt
, and
T.
Bučko
, “
Anharmonic correction to free energy barriers from DFT-based molecular dynamics using constrained thermodynamic integration
,”
J. Chem. Theory Comput.
19
,
2455
2468
(
2023
).
80.
C. H.
Bennett
, “
Efficient estimation of free energy differences from Monte Carlo data
,”
J. Comput. Phys.
22
,
245
268
(
1976
).
81.
R. W.
Zwanzig
, “
High-temperature equation of state by a perturbation method. I. Nonpolar gases
,”
J. Chem. Phys.
22
,
1420
1426
(
1954
).
82.
C.
Chipot
and
A.
Pohorille
, “
Calculating free energy differences using perturbation theory
,” in
Free Energy Calculations
(
Springer
,
Berlin, Heidelberg
,
2007
).
83.
A.
Pohorille
,
C.
Jarzynski
, and
C.
Chipot
, “
Good practices in free-energy calculations
,”
J. Phys. Chem. B
114
,
10235
10253
(
2010
).
84.
J. P.
Perdew
,
A.
Ruzsinszky
,
G. I.
Csonka
,
O. A.
Vydrov
,
G. E.
Scuseria
,
L. A.
Constantin
,
X.
Zhou
, and
K.
Burke
, “
Restoring the density-gradient expansion for exchange in solids and surfaces
,”
Phys. Rev. Lett.
100
,
136406
(
2008
).

Supplementary Material