A Novel NPT Thermodynamic Integration Scheme to Derive Rigorous Gibbs Free Energies for Crystalline Solids

Abstract

Thermodynamic Integration (TI) is the state-of-the-art computational technique for accurate Gibbs free energy predictions of solids. Conventional TI schemes start from an NVT harmonic reference and require three successive corrections to recover the Gibbs free energy of the real crystal in the NPT ensemble. However, the NVT-to-NPT correction neglects full cell flexibility. Here, we present a rigorous (and only) two-step TI scheme that operates entirely in the NPT ensemble, eliminating the need for the approximate NVT-to-NPT step. The key methodological advancement is the novel NPT reference that explicitly accounts for full cell fluctuations. The new approach is compared with the conventional one via two complementary case studies. For ice polymorphs, having simple cell-shape distributions, the new approach reproduces conventional TI results with excellent agreement. For CsPbI3, whose black phase exhibits complex cell-shape behavior, we demonstrate that our novel method provides more accurate Gibbs free energy differences than the conventional one. Moreover, the proposed framework maintains comparable computational cost while offering a simplified workflow. Overall, the new NPT TI scheme provides rigorous and direct Gibbs free energy calculations for solids.

Figures (6)

• Figure 1: Overview of the case study materials. The first case study covers three ice phases (red atoms indicate oxygen; white atoms indicate hydrogen). The second case study comprises the black and yellow phase of CsPbI3 (cyan atoms indicate Cs; purple atoms indicate I; black atoms indicate Pb).
• Figure 2: Schemes for computing Gibbs free energies via thermodynamic integration (TI). Purple arrows indicate the conventional scheme, in which the NVT harmonic reference is first corrected towards the true potential energy surface (PES). A second approximate correction is performed from isochoric to isobaric conditions and, finally, a third from low to high temperatureCheng2018. Green arrows indicate the new scheme, where the reference is a newly derived NPT harmonic PES. Subsequent corrections account for anharmonicity and temperature effects, all performed under constant-pressure conditions.
• Figure 3: Overview of the computational workflow of the new TI scheme. (i) A full geometry optimization yields the optimized atomic positions and simulation cell. (ii) At this minimum, the extended Hessian is constructed via second-order derivatives of deformed and cell coordinates. The extended Hessian defines the NPT harmonic PES, and its non-zero eigenvalues yield the reference Gibbs free energy. (iii) For several $\lambda$-values in the interval $[0,1]$ an NPT MD simulation is run in the mixed PES (Eq. \ref{['eq:mixedHam']}). The $\lambda$-TI correction is found via trapezoidal integration. Overall, two inexpensive steps yield the reference Gibbs free energy, while the computationally more demanding TI-step corrects towards the real value.
• Figure 4: Overview of the applied TI schemes.Ice-conventional: NVT reference and both corrections are computed at the final temperature. Ice-new: NPT reference and single correction are computed at the final temperature and pressure. CsPbI3-conventional: Each desired temperature yields an average cell, yielding four references and four NVT $\lambda$-TI corrections, computed at 600 K to properly sample all tilted configurations of the PbI6 octahedra. NVT temperature corrections are applied using Replica Exchange (REX) MD to ensure proper sampling at low temperatureEarl2005ParallelTempering. The NVT-to-NPT correction is applied at four different temperatures, via histograms accumulated during NPT REX MD. CsPbI3-new: A single reference at 600 K is computed, followed by a single NPT $\lambda$-TI correction. For the lower temperatures, NPT temperature corrections are applied using REX.
• Figure 5: Gibbs free energy differences and their contributions for ice. Results obtained via the conventional method and the new method, at various pressures and temperatures. The conventional method has five free energy contributions, whereas the new method has only three. Positive contributions are stacked above the zero axis and negative contributions below. In all cases, the total Gibbs free energy difference $\Delta G_\text{conv}$ and $\Delta G_\text{new}$ agree extremely well, as expected for materials with simple cell-shape distributions.