Ab initio determination of phase stabilities of dynamically disordered solids: rotational C2 disorder in Li2C2

Johan Klarbring; Stanislav Filippov; Ulrich Häussermann; Sergei I. Simak

PMC · DOI:10.1038/s41598-026-43795-z·March 12, 2026

Ab initio determination of phase stabilities of dynamically disordered solids: rotational C2 disorder in Li2C2

Johan Klarbring, Stanislav Filippov, Ulrich Häussermann, Sergei I. Simak

PDF

Open Access

TL;DR

This paper explains how a material's structure changes with temperature using advanced simulations, focusing on a specific compound's phase transition.

Contribution

A novel ab initio method is introduced to determine phase stabilities in dynamically disordered solids.

Findings

01

The free energy difference between ordered and disordered phases of Li2C2 is calculated using AIMD thermodynamic integration.

02

The entropy stabilizing the dynamically disordered cubic phase is captured through stress behavior analysis.

Abstract

The temperature-induced orthorhombic to cubic phase transition in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} $Li_{2} C_{2}$ \end{document} is a prototypical example of a solid to solid phase transformation between an ordered phase, which is well described within the phonon theory, and a dynamically disordered phase with rotating molecules, for which the standard phonon theory is not applicable. The transformation in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek}…

Linked entities

Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.

Chemicals2

C2

Li2C2

Figures5

Click any figure to enlarge with its caption.

Illustration of the deformation path (parametrized by $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ ) from the ordered orthorhombic (left) to the dynamically disordered cubic (right) phase with the supercells used in this work. Li and C atoms are colored green and brown, respectively. In the cubic phase a snapshot of the positions of C atoms are taken from an AIMD simulation, in order to

Thermal expansion of the *a*, *b* and *c* lattice parameters and volume of the orthorhombic structure of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Li}_2\hbox {C}_2$$\end{document}$ . The thermal expansion is evaluated as $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \set

Results of the stress-strain thermodynamic integration on the deformation path from the ordered orthorhombic phase at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ = 0 to the dynamically disordered cubic phase at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek

Self-correlation function $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{self}(t)/C_{self}(0)$$\end{document}$ for a select number of points along the deformation path between the orthorhombic and cubic phases at 600 K.

Internal energy change $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta U$$\end{document}$ , along the deformation path from the ordered orthorhombic phase at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{documen

Funding1

—Linköping University

Keywords

Phase transformationThermodynamic integrationDynamical disorderAb initio molecular dynamicsChemistryMaterials sciencePhysics

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsSuperconductivity in MgB2 and Alloys · Boron and Carbon Nanomaterials Research · MXene and MAX Phase Materials

Full text

Introduction

The ability to accurately predict and understand solid-solid phase transformations is a longstanding goal in theoretical solid state physics. Predicting and ultimately controlling the external conditions (temperature (T), pressure (P) or applied fields) at which these phase transformations occur is of great importance for a host of applications. Indeed, the transformations may be detrimental to performance, or may be the phenomenon on which a device could be based, as eg. in mechanocaloric materials for solid-state cooling^1,2^ and shape-memory alloys^3,4^.

In order to predict the critical temperature of the phase transformations, accurate determination of the relevant free energies is required. First-principles electronic structure methods based on density functional theory (DFT) have been very successful in predicting pressure induced phase transformations at low temperatures, where the free energy is well approximated by the enthalpy.

At relatively low T (as compared to melting) many solids are well described within the harmonic approximation. That is, the atomic vibrations in the system are well represented by a set of non-interacting harmonically oscillating normal-modes, and accurate predictions of temperature-induced structural transitions can thus be made, see eg. Ref.^5^. This approach is particularly accurate if extended to the so-called quasi-harmonic approximation, where harmonic phonons are calculated at thermally expanded volumes.

In many cases, however, the (quasi-)harmonic approximation is insufficient. For instance for systems that are unstable with respect to the equilibrium atomic positions of some high symmetry phase at 0 K, but where these positions are stabilized by anharmonic phonon-phonon interactions at elevated temperatures. In such anharmonic cases, re-normalized phonon schemes, including self-consistent phonons^6–9^, and effective harmonic models fitted from molecular dynamics (MD)^10–12^, have proven to be successful.

Corrections to the free energy of an (effective) harmonic model from higher order terms in the phonon expansion can also be derived^13^ and have been applied in a first-principles context^14^.

Beyond these “phonon-based” techniques, the concept of thermodynamic integration (TI) provides a set of formally exact methods to extract free energies, which are increasingly being used in connection with first-principles DFT-related techniques^15–17^. While machine learning force fields (MLFFs) are emerging as powerful tools for such studies , accurately capturing the delicate meV-scale energy balances in dynamically disordered systems without extensive specialized training sets remains challenging. Thus, we rely on fully first-principles AIMD to ensure our thermodynamic integration provides an unbiased, rigorously exact reference. Methods in this class are based on the fact that while the full free energy, or more precisely the entropy, is very difficult to extract directly from an MD simulation, certain derivatives of the free energy can be more easily extracted and the free energy itself can then be calculated by integration over several different simulations. This is founded on the basic result that free energy derivatives can be evaluated as thermal averages over energy derivatives, i.e. for an external parameter $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ coupled to the system it holds that^18^

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left( \frac{\partial F(V,T,\lambda )}{\partial \lambda } \right) = \overline{ \frac{\partial H(\{\textbf{q},\textbf{p}\};\lambda )}{\partial \lambda }}, \end{aligned}$$\end{document}

where H is he Hamiltonian of the sytem. The phase-space ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\textbf{q},\textbf{p}\}$$\end{document}$ ) ensemble average on the right-hand-side is practically evaluated as a time-average over an MD simulation under the assumption of ergodicity.

The most commonly used TI method is the so called Kirkwood coupling constant integration. Here, one adds to the free energy a parametrical dependence on an external parameter $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ , i.e., $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = F(V,T,\lambda )$$\end{document}$ and the free energy change upon changing $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ is obtained by integration of the right-hand-side of Eq. (1).

A less commonly employed TI technique involves obtaining free energy changes by integrating pressure over a volume change, or more generally, as will be the basis of the calculations in this work, integrating a stress related quantity over a deformation path connecting two (real or artificial) phases of interest^19–21^.

A certain set of solids, which we will refer to as dynamically disordered solids, necessitates the use of free energy methods beyond phonon based ones. These are systems which possess some type of time-averaged long-range order, characteristic of a crystalline solid, but where well-defined, unique, equilibrium positions cannot be statically assigned to each individual atom.

Fig. 1. Illustration of the deformation path (parametrized by $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ ) from the ordered orthorhombic (left) to the dynamically disordered cubic (right) phase with the supercells used in this work. Li and C atoms are colored green and brown, respectively. In the cubic phase a snapshot of the positions of C atoms are taken from an AIMD simulation, in order to illustrate the rotational disorder of the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimers, while the Li atoms are displayed in their high-symmetry positions. The conventional orthorhombic unit cell is illustrated by solid black lines. The Cartesian reference coordinate system, $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\textbf{x},\textbf{y},\textbf{z}\}$$\end{document}$ , used in our simulations is indicated by black vectors and letters, while the coordinate system, $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\mathbf {x'},\mathbf {y'},\mathbf {z'}\}$$\end{document}$ , aligned with the conventional orthorhombic cell is displayed in red. $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{z}$$\end{document}$ and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{z}'$$\end{document}$ are collinear and directed out of the figure. The spheres at the bottom schematically represent the angular distribution of the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimers over AIMD simulations at 600 K. Bright colors and dark colors correspond to high density and low density, respectively.

They involve certain fast ion-conducting solids, so called superionics, which posses a “liquid-like” diffusional behavior of certain ionic species^22^. Such systems are being increasingly intensively studied, with envisioned applications ranging from solid electrolytes in batteries and fuel cells to thermoelectric materials^23,24^.

Another clear example of dynamically disordered solids are those containing rotating molecular units, sometimes referred to as plastic or orientation-disordered solids. Such systems typically have a high temperature phase where the center of mass of each molecular unit is ordered on certain high symmetry positions, but where the individual atoms making up the molecule are positionally disordered. To cite an emerging application, we mention that some such systems are highly promising barocaloric materials for solid-state cooling applications^1,2,25–27^.

In this work, we investigate from first principles the temperature induced orthorhombic-to-cubic phase transformation of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Li}_2\hbox {C}_2$$\end{document}$ at ambient pressure. We choose $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Li}_2\hbox {C}_2$$\end{document}$ as it is a very “clean” example of a system with a temperature-induced phase-transformation between an ordered and a dynamically disordered phase with rotating molecules. The low temperature orthorhombic phase of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Li}_2\hbox {C}_2$$\end{document}$ is ordered with $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimers aligned along the orthorhombic b axis. Upon heating, this phase transforms into the cubic phase, where the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimers are dynamically disordered among shallow local minima corresponding to alignment along the [110] directions of the cubic structure^28^. Figure 1 shows the orthorhombic and cubic phases and the nature of the dimer alignment in these phases. We employ the stress-strain TI framework to investigate the nature of this transformation. We find that the stabilizing entropy associated with the onset of disorder among the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimers is seen through the behavior of the stresses on the deformation path between the two phases, and fair agreement with experiments for the critical temperature of the transition.

Methodology

\documentclass[12pt]{minimal}

			\usepackage{amsmath}
			\usepackage{wasysym} 
			\usepackage{amsfonts} 
			\usepackage{amssymb} 
			\usepackage{amsbsy}
			\usepackage{mathrsfs}
			\usepackage{upgreek}
			\setlength{\oddsidemargin}{-69pt}
			\begin{document}$$\hbox {Li}_2\hbox {C}_2$$\end{document} polymorphs and stress-strain thermodynamic integration

Figure 1 illustrates the supercells used to represent the low-temperature orthorhombic and high temperature cubic phases. The relation between the monoclinic and the conventional orthorhombic representation of the low-temperature phase is displayed. Viewed in the coordinate system $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{ \textbf{x}',\textbf{y}',\textbf{z}' \}$$\end{document}$ , aligned with the orthorhombic axes, the transformation is a compression along $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{y}'$$\end{document}$ and an expansion along $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{x}'$$\end{document}$ and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{z}'$$\end{document}$ . The supercell matrices used for the two phases in this work are given in the Appendix.

The stress-strain TI methodology starts from the differential of the Helmholtz free energy, F, for a rotationally invariant crystal^20,29^,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} dF(\{\textbf{X},\boldsymbol{\eta }\},T) = -SdT + V(\textbf{X}) \boldsymbol{\tau } \mathbin {\mathbf {:}} d\boldsymbol{\eta } . \end{aligned}$$\end{document}

where S and T are the entropy and temperature, $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{X}$$\end{document}$ is an arbitrary reference configuration, $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\eta }$$\end{document}$ is the Lagrangian strain tensor, defined w.r.t $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{X}$$\end{document}$ , and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V(\textbf{X})$$\end{document}$ and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\tau }$$\end{document}$ are the second Piola-Kirchoff stress tensor, respectively.

Integrating along a deformation path parametrized by $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ at constant T finally yields the free energy change as^20^

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Delta F(\lambda ) = \int _0^\lambda V(\lambda ') \boldsymbol{\sigma }(\lambda ') \boldsymbol{h}^{-T}(\lambda ') \mathbin {\mathbf {:}} \frac{\partial \boldsymbol{h}(\lambda ')}{\partial \lambda '} d \lambda ', \end{aligned}$$\end{document}

where $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\sigma }$$\end{document}$ is the Cauchy stress tensor, related to $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\tau }$$\end{document}$ as $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\tau } = \mathbin {J} \boldsymbol{\alpha }^{-1} \boldsymbol{\sigma } \boldsymbol{\alpha }^{-T}$$\end{document}$ . The deformation gradient $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\alpha }$$\end{document}$ is the derivative of each component of a material point in the configuration $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{x}$$\end{document}$ with respect to each component of the material point in the reference configuration $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{X}$$\end{document}$ , and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J = \det (\boldsymbol{\alpha })$$\end{document}$ . For $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{x}$$\end{document}$ = $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{x}(\boldsymbol{X})$$\end{document}$ , the elements of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\alpha }$$\end{document}$ are

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \alpha _{ij} = \frac{\partial x_i}{\partial X_j}. \end{aligned}$$\end{document}

We define the matrix $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{h}$$\end{document}$ to contain as its columns the supercell lattice vectors $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{a}$$\end{document}$ , $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{b}$$\end{document}$ and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{c}$$\end{document}$ used in the AIMD simulations,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \boldsymbol{h} = \begin{pmatrix} a_x & b_x & c_x \\ a_y & b_y & c_y \\ a_z & b_z & c_z \\ \end{pmatrix}. \end{aligned}$$\end{document}

For homogeneous deformations, which is what we treat here, $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\alpha }$$\end{document}$ relates the lattice vectors in the initial configuration $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{X}$$\end{document}$ and a configuration $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{x}$$\end{document}$ as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \boldsymbol{h}(\boldsymbol{x}) = \boldsymbol{\alpha } \boldsymbol{h}(\boldsymbol{X}). \end{aligned}$$\end{document}

$\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\eta }$$\end{document}$ is related to $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\alpha }$$\end{document}$ as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \boldsymbol{\eta } = \frac{1}{2} \left( \boldsymbol{\alpha }^T\boldsymbol{\alpha } - \boldsymbol{I} \right) , \end{aligned}$$\end{document}

where $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{I}$$\end{document}$ is the identity matrix.

It is important to notice that in thermodynamic integration as employed here, the Helmholtz free energy $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F(\lambda )$$\end{document}$ is evaluated along a strain path $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda \in [0,1]$$\end{document}$ between two structures, where $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ parametrizes a smooth deformation of the supercell matrix. Since we work in a finite supercell with periodic boundary conditions, the system has a finite number of degrees of freedom. The canonical partition function is then given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z(\lambda , \beta = \frac{1}{k_B T}) = \int e^{-\beta H(\lambda )}\, d\Gamma ,$$\end{document}

where $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H(\lambda )$$\end{document}$ is the Born–Oppenheimer potential energy surface, and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\Gamma$$\end{document}$ denotes the configurational phase-space measure. In this finite system setting, it is a well-established result from statistical mechanics that $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z(\lambda , \beta )$$\end{document}$ and hence the free energy $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F(\lambda ) = -\frac{1}{\beta } \ln Z(\lambda ,\beta )$$\end{document}$ are analytic functions of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ , provided that the Hamiltonian depends smoothly on $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ . This analyticity was rigorously demonstrated for quantum systems with bounded interactions by Bogoliubov^30^ and extended to physically relevant unbounded interactions (such as Coulomb potentials) by Maison^31^. While these results were originally derived for quantum systems, the analyticity property extends to classical systems through the correspondence principle, as the classical partition function can be viewed as the high-temperature limit of the quantum case. In our AIMD simulations, the only $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ -dependence enters through the smooth strain deformation of the potential energy surface: as $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ changes, the lattice vectors $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{h}(\lambda )$$\end{document}$ vary linearly. For fixed phase-space coordinates (e.g., fixed fractional atomic positions), this induces a smooth and differentiable $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ -dependence of the Hamiltonian through the cell matrix. Since the interatomic potentials are smooth functions of the resulting interatomic distances, the total potential energy surface $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V(\lambda )$$\end{document}$ is smooth by composition of smooth functions. Therefore, both $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F(\lambda )$$\end{document}$ and its derivative $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial F / \partial \lambda$$\end{document}$ (the generalized stress) are guaranteed to be continuous and differentiable functions of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ , even when the system undergoes dramatic microscopic changes such as the onset of rotational disorder of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dumbbells. While such transitions may involve rapid increases in both potential energy (due to less favorable configurations) and entropy (due to many more accessible states), these changes occur through continuous crossovers rather than true phase transitions in the finite supercell. The free energy $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = U - TS$$\end{document}$ remains smooth because the partition function $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z(\lambda ) = \int e^{-\beta H(\lambda )} d\Gamma$$\end{document}$ is a smooth integral over the finite configurational space, despite the system exploring qualitatively different dynamical regimes at different $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ values. Consequently, if any discontinuities or jumps in the stress–strain curve are observed in actual calculations, they cannot reflect true thermodynamic singularities, but instead arise from insufficient equilibration between distinct dynamical regimes. This systematic error can be diagnosed by performing reverse integration and checking for path independence. The transition region may require longer simulation times to properly sample both ordered and disordered configurational states, but the underlying free energy surface remains mathematically well-defined and integrable. Even if the precise $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ value at which dynamic disorder sets in is not fully resolved due to finite simulation times, the resulting error in the integrated free energy difference is proportional to the area under the stress error–strain curve in the poorly equilibrated region. Because this region is narrow in $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ , the associated contribution to the total free energy difference remains negligible compared to the overall free energy change between the ordered and disordered phases. Finally, because the supercell is finite, true phase coexistence is impossible, ensuring that the system remains in a single, well-defined thermodynamic state throughout the integration path.

We find it revealing to analyse the Cauchy stress in the coordinate system aligned with the conventional orthorhombic cell, $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{ \textbf{x}',\textbf{y}',\textbf{z}' \}$$\end{document}$ , see Fig. 1. The Cauchy stress $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\sigma }$$\end{document}$ transforms as $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\sigma }' = \textbf{Q}\boldsymbol{\sigma }\textbf{Q}^{T}$$\end{document}$ , where $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{Q}$$\end{document}$ is the transformation matrix from the coordinate system $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\textbf{x},\textbf{y},\textbf{z}\}$$\end{document}$ to $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{ \textbf{x}',\textbf{y}',\textbf{z}' \}$$\end{document}$ given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{Q} = \begin{pmatrix} \cos \theta & -\sin \theta & 0\\ \sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{pmatrix}, \end{aligned}$$\end{document}

where $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document}$ is defined in Fig. 1. Note that $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document}$ and hence $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{Q}$$\end{document}$ changes along the deformation path, i.e., they are $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ dependent.

Computational details

We carried out Kohn-Sham (KS) DFT calculations using the projector augmented wave (PAW)^32^ method and a plane wave basis set as implemented in the Vienna ab initio simulation package (VASP)^33–35^. Exchange and correlation effects were treated using the Perdew-Burke-Ernzerhof (PBE)^36^ form of the generalized gradient approximation (GGA). PAW potentials with Li(1s2s) and C(2s2p) states treated as valence were used. The plane-wave cutoff energy of the KS orbitals was 600 eV. We used a 128 atom supercell constructed as a 2 $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document}$ 2 $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document}$ 2 expansion of the conventional cubic anti-fluorite unit cell, for which k-point sampling was performed at the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma$$\end{document}$ -point. This k-point sampling is enough to converge the stress tensor within $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}$ 1 kbar, as judged by comparing to $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2 \times 2 \times 2$$\end{document}$ and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3 \times 3 \times 3$$\end{document}$ k-point grids on 10 snapshots from AIMD simulations. AIMD simulations in the NVT ensemble were performed using a short 0.25 fs timestep and a Nosé-Hoover thermostat with a Nosé mass corresponding to a time-constant of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}$ 86 timesteps. The presence of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_2$$\end{document}$ dimer stretching modes at around 60 THz, separated from the rest of the vibrational modes in the system, makes it necessary to carefully choose the timestep and the Nosé mass in order to avoid having the C and Li subsystems slowly decoupled and equilibrate at two separate temperatures over long timescales.

This was observed to happen using more standard parameters, e.g. a 1 fs timestep and the default VASP value for the Nosé mass. The combined simulation time at each point on the deformation path was typically at least $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}$ 360000 timesteps, corresponding to 90 ps of simulation time, and for select cases significantly longer to judge convergence. We typically skip at least $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}$ 40 ps as equilibration. More rationale on the simulation times follows below.

We also note that while the high-frequency $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ stretching mode possesses a non-negligible quantum zero-point energy (ZPE), the intramolecular covalent C-C bond remains completely intact with virtually identical stretching frequencies in both the orthorhombic and cubic phases. Consequently, the ZPE contributions effectively cancel out when calculating the free energy difference between the two phases, justifying our classical treatment of the nuclei for the phase stability calculation.

Results and discussion

Thermal expansion

We start by investigating the temperature dependence of the lattice constants in the orthorhombic and cubic phases of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Li}_{2}\hbox {C}_{2}$$\end{document}$ . We extract these from AIMD by iteratively changing the lattice constants until all elements of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\sigma }$$\end{document}$ are <2 kBar (more details on this are discussed below). In practice this is done by taking the average stress from an AIMD run in the monoclinic cell, transforming it to the orthorhombic cell and changing the a, b, and c lattice vectors according to the diagonal components of this transformed stress tensor in a steepest decent fashion. The procedure is then repeated until convergence.

The obtained values at 300 K are $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a =$$\end{document}$ 3.68 Å, $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b =$$\end{document}$ 4.83 Å, and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c =$$\end{document}$ 5.44 Å, which agree very well with the experimental values at room temperature^37^. Figure 2 shows the thermal expansion from T = 300 K, i.e., $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(a(T)-a(300 \text { K}))/a(300 \text { K})$$\end{document}$ and similar for b, c and the volume V. While the a and c axes expand with temperature, the b axis, which is collinear with the alignment of the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_2$$\end{document}$ dimers, shows a non-conventional decrease at high temperature. This is an effect of the temperature induced weakening of the alignment of the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_2$$\end{document}$ dimers. This peculiarity was pointed out experimentally for $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Li}_2\hbox {C}_2$$\end{document}$ ^37^.

Fig. 2. Thermal expansion of the a, b and c lattice parameters and volume of the orthorhombic structure of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Li}_2\hbox {C}_2$$\end{document}$ . The thermal expansion is evaluated as $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(a(T)-a(300 \text { K}))/a(300 \text { K})$$\end{document}$ for a and similar for b, c, and the volume. Inset shows the conventional orthorhombic unit cell.

Stress-strain thermodynamic integration

We obtain the free energy difference $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta F(\lambda )$$\end{document}$ on a deformation path between the ordered orthorhombic phase ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda = 0$$\end{document}$ ) and the disordered cubic phase ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda = 1$$\end{document}$ ) using Eq. (3). We use 10 points obtained by linear interpolation between $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{h}(\lambda = 0)$$\end{document}$ and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{h}(\lambda = 1)$$\end{document}$ . Since the path connects an ordered phase with a dynamically disordered one, the disorder is expected to set in somewhere along the path. The convergence of the stress with respect to simulation time thus becomes a critical issue. The approach we take to attempt to ensure proper convergence is the following: At each temperature and each point on the deformation path we perform two separate AIMD runs, one with initial (fractional) atomic positions and velocities taken from an AIMD snapshot of the orthorhombic phase, and one with a snapshot from the cubic phase. Convergence can then be gauged by how well these two sets of stress-strain TIs agree with each other, as a function of simulation time. Indeed, insufficient simulation time in the two sets should produce errors of opposite type. Roughly speaking, the simulations started from the orthorhombic snapshot should be biased towards producing errors associated with too much order, while those started from the cubic snapshot should be biased towards too much disorder. We perform simulations long enough so that the integrand in Eq. (3) differs between the two AIMD runs by less than $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}$ 1.6 meV/f.u at each point on the deformation path. The mean difference over the whole path is significantly smaller than this, around $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}$ 0.5 meV/f.u.. The stresses used to produce the final results are averaged over the two simulations, in order to increase the accuracy further. This two-initialization strategy inherently serves as a rigorous forward-reverse hysteresis check. The fact that the stresses converge to within $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}$ 0.5 meV/f.u. from completely different starting states provides empirical proof that the system has reached ergodic equilibrium at each point along $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ , and establishes a strict upper bound on the statistical error of our thermodynamic integration.

Figure 3 (a) shows the free energy differences between the two phases at 600 and 650 K. We see that at 600 K the orthorhombic phase is favored by 1 meV/f.u., while at 650 K the cubic phase is 4 meV/f.u. lower in free energy. A linear interpolation results in a phase transition temperature of 611 K. Comparing this temperature to an experimental transition temperature is not straightforward. It was noted^37^ that the first appearance of an XRD reflection corresponding to the high temperature phase upon heating appears at around 693 K, and differential scanning calorimetry (DSC) measurements^38^ puts the transition at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}$ 725 K. Our thermodynamic transition temperature is thus likely underestimated, although in fair agreement considering all approximations inherent to the simulations, including the exchange-correlation functional and the finite and defect-free supercells. We note that this also highlights a well known difficulty in modelling solid-solid phase transitions, namely that small shifts in the free energy difference between the phases, can lead to substantial shifts in the predicted transition temperature. As is clear from Fig. 3, we clearly deal with small free energy differences, and the methodology is thus certainly susceptible to this difficulty.

Fig. 3. Results of the stress-strain thermodynamic integration on the deformation path from the ordered orthorhombic phase at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ = 0 to the dynamically disordered cubic phase at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ = 1. The top panel shows the free energy change along the deformation path at 600 K (solid blue line) and 650 K (dashed red line). Middle and lower panel show the diagonal components of the Cauchy stress in the coordinate system $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{ \textbf{x}',\textbf{y}',\textbf{z}' \}$$\end{document}$ aligned with the conventional orthorhombic unit cell, at 600 and 650 K, respectively.

It has also noted^37^, however, that the two phases continue to co-exist over a broad temperature range. This is in agreement with the result from Fig. 3 (a) that the two phases appear to be free energy local minima both above and below the transition temperature.

The slight non-zero values of the stresses at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ = 0 and 1 are due to the fact that the configurations chosen as starting points, turn out, once very long simulations have been run, to be slightly offset from the equilibrium ones. That is, the path we choose through configuration space passes slightly to the side of the actual minima. The change in free-energy in adjusting for this is, however, negligible, in the “worst” case, which is $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ = 1 at 650 K, this free energy correction can be estimated to be below 0.5 meV/f.u.

In order to gauge how the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimer disorder sets in on the deformation path we calculate a dimer self-correlation function,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} C_{self}(t) = \langle \tilde{\textbf{S}}_i(t_0) \cdot \tilde{\textbf{S}}_i(t+t_0) \rangle _{N,t_0}, \end{aligned}$$\end{document}

where $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\textbf{S}}_i(t) = \textbf{S}_i(t) - \langle \textbf{S}_i(t)\rangle _N$$\end{document}$ , $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S}_i(t)$$\end{document}$ is the vector between the two C atoms making up $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimer i and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \cdots \rangle _{N,t_0}$$\end{document}$ denotes an average over both the N $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimers and all time-origins $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_0$$\end{document}$ .

In Fig. 4 we plot $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{self}(t)/C_{self}(0)$$\end{document}$ for select values of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ along the deformation path at 600 K. In the orthorhombic phase at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ = 0 we observe, as expected, a quick initial decay in the very short time limit followed by a virtually constant behaviour, consistent with the ordered state of the dimers. At the first point along the path ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ = 0.09) a similar initial decay is followed by a slow decay, indicating a still largely ordered state of the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimers. Moving further along the deformation path, the behaviour qualitatively changes into significantly faster decay rates. Around this point on the deformation path, we see a change in the behaviour of the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{yy}$$\end{document}$ stress component from the standard decrease upon compression, to an increase, which we can thus associate with the onset of proper disorder of the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimers, which effectively weakens the alignment of the dimers along the y direction.

Starting from $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ = 0.27, $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{self}(t)/C_{self}(0)$$\end{document}$ is reasonably well fitted to a stretched exponential $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e^{-(t/\tau )^\beta }$$\end{document}$ , with $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau \approx$$\end{document}$ 2.9 ps and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta \approx$$\end{document}$ 0.67 at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ = 0.27. For $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ approaching 1, i.e. the cubic structure, we can instead fit to a single exponential $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e^{-(t/\tau )}$$\end{document}$ . For the cubic phase we obtain a relaxation time of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau \approx$$\end{document}$ 0.69 ps.

Fig. 4. Self-correlation function $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{self}(t)/C_{self}(0)$$\end{document}$ for a select number of points along the deformation path between the orthorhombic and cubic phases at 600 K.

Figure 5 shows the changes in internal energy U along the deformation paths at 600 and 650 K. We see that the difference in U between the two phases amounts to $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}$ 82 and 71 meV/f.u. at 600 and 650 K, respectively. At the phase-transformation temperature this gives, by linear interpolation, a transition enthalpy of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}$ 80 meV/f.u. It is thus clear that the cubic phase is stabilised by a large entropy contribution originating from the disordered $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimers. It is interesting to note that U is increasing for $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ across 0, at 600 K the free energy minima at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ = 0 is at least 8 meV/f.u. higher in U compared to $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda \approx$$\end{document}$ -0.09, which implies that there is a quite significant entropic part responsible for producing the free energy minima even in the orthorhombic case. Moving towards lower values of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ corresponds to elongation of the orthorhombic b axis and contraction of a and c, which will act to suppress librational motion of the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimers, i.e., to increase their directional order. It thus follows that the position of cell-configuration that corresponds to the orthorhombic free energy minimum is a competition between the energetic and entropic tendencies to suppress and promote $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ librations, respectively.

Fig. 5. Internal energy change $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta U$$\end{document}$ , along the deformation path from the ordered orthorhombic phase at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ = 0 to the dynamically disordered cubic phase at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}$ = 1. Results at temperatures 600 and 650 K are shown with blue circles and red squares, respectively.

Since we at this point have access to both the free and internal energy differences between the ordered and disorder phases we can extract their entropy difference $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta S$$\end{document}$ . At 600 K and 650 K, T $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta S$$\end{document}$ corresponds to $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}$ 80 and 75 meV/f.u., respectively, i.e $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta S$$\end{document}$ corresponds to 0.1333 and 0.1153 meV/K/f.u., respectively.

Thus both the energy and the entropy differences between the two phases are larger at the lower temperature. This may look counter-intuitive, and lead to conceptualization that the transition is not happening due to a rapid increase of the entropy with temperature, but instead happens primarily due to an increase of the energy of the orthorhombic phase. Let us show this would be an incomplete statement, since the changes of energy and entropy with temperature are not independent.

If we consider the system at constant pressure and neglect the volume change with temperature, which is a reasonable approximation when comparing the states at 600 and 650 K, the thermodynamic equality following directly from a fundamental equation of thermodynamics reads

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left( \frac{\partial U}{\partial T}\right) _{V} = T \left( \frac{\partial S}{\partial T}\right) _{V}. \end{aligned}$$\end{document}

Therefore, for a small temperature change from $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {T}_{1}$$\end{document}$ to $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {T}_{2}$$\end{document}$ we can write in finite differences

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} U(T_{2}) - U(T_{1}) = T_{1} [S(T_{2}) - S(T_{1})]. \end{aligned}$$\end{document}

Equation (11) can be easily re-written for the Helmholtz free energy as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} F(T_{2}) = F(T_{1}) - (T_{2} - T_{1}) S(T_{2}). \end{aligned}$$\end{document}

If we define a difference between thermodynamic functions of the two phases, dynamically ordered (I) and disordered (II), as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Delta X^{II-I}(T_{i}) = X^{II}(T_{i}) - X^{I}(T_{i}), \end{aligned}$$\end{document}

where X can be F, U, or S, and i stands for either 1 or 2, then the equation for the change in the Helmholtz free energy difference between the two phases upon changing temperature from $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {T}_{1}$$\end{document}$ to $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {T}_{2}$$\end{document}$ is given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Delta F^{II-I}(T_{2}) = \Delta F^{II-I}(T_{1}) - (T_{2} - T_{1}) \Delta S^{II-I}(T_{2}). \end{aligned}$$\end{document}

Thus, if the free energy difference between the dynamically ordered and disordered phases is positive (as we have at 600 K, i.e. $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta F^{II-I}(T_{1}) > 0$$\end{document}$ ), then the only necessary condition for the transition into the dynamically disordered phase with increasing temperature is the positive $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta S^{II-I}(T_{2})$$\end{document}$ , and the size of this $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta S^{II-I}(T_{2})$$\end{document}$ defines how high should be the transition temperature. We notice that such a simple estimation is very accurate in the considered case. For $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {T}_{1}$$\end{document}$ = 600 and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {T}_{2}$$\end{document}$ = 650 K, the right-hand side of Eq. (14) gives – 3.8 meV, while its left-hand side calculated directly from TI gives – 4 meV.

It is also interesting to compare the entropy between the two phases with a very simple model where the entropy difference is taken to be purely associated with a set of static configurations of the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_2$$\end{document}$ dimers. In the high temperature limit, we then have $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta S_{conf} = k_{B}\log (6)$$\end{document}$ , which gives a contribution T $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta S_{conf} \approx$$\end{document}$ 93 meV/f.u. at T = 600 K. This is larger, but comparable, to the actual entropy contribution to the free energy difference between the two phases at 600 K of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}$ 80 meV/f.u.

The $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimers clearly have a preference of aligning along the [110] directions of the cubic cell, as seen in Fig. 1. However, it can also be seen from the figure that the density of dimer-alignments on paths between different [110] directions is clearly non-negligible. Furthermore, the self-relaxation time of the dimers at 600 K was estimated above to be 0.69 ps, which is comparable to the periods of atomic vibration in the system. These considerations indicate that the conceptual partitioning of the total entropy into vibrational and configurational parts, $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{tot} = S_{vib.} + S_{conf.}$$\end{document}$ , as is commonly done for solids, is questionable in $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Li}_2\hbox {C}_2$$\end{document}$ , and likely other similar systems. Thus estimating, as one could typically do, the free energy by averaging the energy and harmonic vibrational entropy over many stable configurations of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimers and then adding an estimate of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{conf}$$\end{document}$ on top of this, is not advisable for these types of systems. Instead, the free energy (or entropy) should preferably be calculated by a method which does not need to make this type of entropy partitioning, such as the stress-strain TI applied in this work.

Conclusion

To conclude we have studied the phase transformation from the low-temperature ordered orthorhombic phase to the high-temperature, dynamically disordered, cubic phase of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Li}_2\hbox {C}_2$$\end{document}$ using the stress-strain thermodynamic integration (TI) methodology. The stabilizing entropy of the cubic phase, associated with the rotational disorder of the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {C}_2$$\end{document}$ dimers, is captured through the behaviour of the stress on the deformation path that connects the two phases. We note good agreement with available experimental results, in particular with respect to the finding that both phases appear to be local free-energy minima above as well as below the phase transition temperature. Our findings indicate that the stress-strain TI methodology can accurately describe phase transformations from ordered phases to phases with dynamical disorder related to molecular units.

Bibliography8

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Escorihuela–Sayalero, C. et al. Prediction and understanding of barocaloric effects in orientationally disordered materials from molecular dynamics simulations. Npj Comput. Mater.10(1) , 13. 10.1038/s 41524-024-01199-5 (2024).
2Jani, J. M., Leary, M., Subic, A. & Gibson, M. A. A review of shape memory alloy research, applications and opportunities. Mater. Design 56, 1078–1113. 10.1016/j.matdes.2013.11.084 (2015).
3Knoop, F. et al. Tdep: temperature dependent effective potentials. J. Open Sourc. Softw.9(94), 6150. 10.21105/joss.06150 (2024).
4Wallace, D. C. Thermodynamics of Crystals (Dover Publications, inc., 1998).
5Landau, L. D. & Lifshitz, E. M. Statistical Physics (Elsevier, 1980).
6Li, B. et al. Colossal barocaloric effects in plastic crystals. Nature 567(7749), 506–510. 10.1038/s 41586-019-1042-5 (2019).10.1038/s 41586-019-1042-530918372 · doi ↗ · pubmed ↗
7Lloveras, P. et al. Colossal barocaloric effects near room temperature in plastic crystals of neopentylglycol. Nat. Commun.10, 11. 10.1038/s 41467-019-09730-9 (2019).10.1038/s 41467-019-09730-9PMC 647242331000715 · doi ↗ · pubmed ↗
8Wallace, D. C. Thermoelastic theory of stressed crystals and higher-order elastic constants. Solid State Phys.25, 301–404 (1970). https://www.sciencedirect.com/science/article/pii/S 0081194708600107.