# Efficient Calculation of the Negative Thermal Expansion in ZrW_{2}O_{8}

- Department of Physics, University of Washington, Seattle, WA, United States

We present a study of the origin of the negative thermal expansion (NTE) on ZrW_{2}O_{8} by combining an efficient approach for computing the dynamical matrix with the Lanczos algorithm for generating the phonon density of states in the quasi-harmonic approximation. The simulations show that the NTE arises primarily from the motion of the O-sublattice, and in particular, from the transverse motion of the O atoms in the W–O and W–O–Zr bonds. In the low frequency range these combine to keep the WO_{4} tetrahedra rigid and induce internal distortions in the ZrO_{6} octahedra. The force constants associated with these distortions become stronger with expansion, resulting in negative Grüneisen parameters and NTE from the low frequency modes that dominate the positive contributions from the high frequency modes. This leads us to propose an anharmonic, two-frequency Einstein model that quantitatively captures the NTE behavior.

## 1. Introduction

There has been considerable interest in recent years in developing materials with electronic and structural properties tuned toward specific applications. Likewise, there is broad interest in pre-screening potential materials using theoretical simulations in an effort to minimize their complexity and simplify their synthesis. Although great progress has been made in the prediction of many properties (Jain et al., 2013) such as band gaps, electronic properties, equilibrium structures, and spectra, thermal properties depend on additional calculations of vibrational properties and often remain difficult for first principles computational methods. This is typically due to the complexity of the materials under consideration. For example, complex ceramics like ZrW_{2}O_{8} with large unit cells have heretofore been too demanding for routine screening using first principles calculations of both electronic and vibrational properties with conventional electronic structure codes. On the other hand, efficient calculations of many vibrational properties are now possible. For example, we have previously shown (Poiarkova and Rehr, 2001; Krappe and Rossner, 2002; Vila et al., 2007) that an efficient Lanczos algorithm for the projected phonon density of states (PDOS) obtained from standard density functional theory (DFT) calculations (Lee and Gonze, 1995; Rignanese et al., 1996; Baroni et al., 2001) of the dynamical matrix (DM) and the quasi-harmonic approximation (Allen and De Wette, 1969; Boyer, 1979) can produce accurate Debye-Waller factors for EXAFS and x-ray crystallography.

Here we develop an extension of this approach for efficient calculations of the thermal properties of ZrW_{2}O_{8}, a quintessential example of negative thermal expansion (NTE) ceramics, in which the NTE is large over a broad range of temperatures.

The origin of the NTE in ZrW_{2}O_{8} is controversial, having been previously attributed to a variety of mechanisms. For instance, the Rigid Unit Mode (RUM) model (Hammonds et al., 1996; Pryde et al., 1996, 1997; Hancock et al., 2004; Tucker et al., 2005, 2007) suggests that the tetrahedra and octahedra that make up the structure are mostly unaffected by the active modes. In contrast, other works (Cao et al., 2002, 2003; Bridges et al., 2014) have suggested that these units are distorted while the Zr-W distance remains more or less unchanged. Yet another alternative (Mary et al., 1996; Gupta et al., 2013; Sanson, 2014) suggests that the Zr–O–W bonds bend, pulling the Zr and W atoms closer and thus causing the NTE. In this paper we take advantage of the local nature of the PDOS to study the origin of the NTE. This procedure allows us to partition the contributions to the Helmholtz free energy from different atomic sites and deduce the origin of NTE from the shifts in the local free energy minima. In a second study in this collection (Vila et al., under review) we apply these methods to study the effects of the NTE on the EXAFS properties of the system, such as the EXAFS mean square relative displacements for relevant single- and multiple-scattering paths (Cao et al., 2003), as well as the anisotropic behavior of the crystallographic mean square atomic displacements. The following sections present a brief summary of the Lanczos algorithm for the PDOS, show how we compute the dynamical matrix efficiently, give a detailed discussion of the volume and temperature effects on the force constants, potentials and PDOS, and present a simple two-frequency model that reproduces the observed NTE quite accurately.

## 2. Methods

### 2.1. Helmholtz Free Energy

The quasi-harmonic approximation (QHA) (Allen and De Wette, 1969; Boyer, 1979) has been widely used to study thermal effects on the structural properties of both crystalline (Erba et al., 2015) and molecular solids (Erba et al., 2016; Brandenburg et al., 2017). QHA simulations usually rely on the full diagonalization of the force constant matrix obtained with density functional perturbation theory (DFPT) or other analytic derivative approaches (Baroni et al., 2010; Erba, 2014). We have previously shown (Vila et al., 2007) that it can also be combined with the Lanczos algorithm in order to avoid the full diagonalization and efficiently compute many thermal properties of materials. Briefly, the lattice thermal expansion *a*(*T*) can be obtained by minimizing the Helmholtz free energy *A*(*a, T*) at a given *T*,

Within the quasi-harmonic approximation *A*(*a, T*) is given by

where *U*(*a*) is the internal energy of the electronic system in its ground state, and *F*(*a, T*) is the total vibrational free energy (VFE) of the system

in the Born-Oppenheimer approximation. For convenience we define *F*(*a, T*) = *F*_{0}(*a*) + *F*_{T}(*a, T*). Here *F*_{0}(*a*) and *F*_{T}(*a, T*) are the zero-point and temperature-dependent components, respectively. They can be calculated in terms of the density of vibrational modes of the lattice ρ(*a*, ω) at a given *a*,

where

is the total PDOS, and ρ_{iα}(*a*, ω) is the PDOS projected onto the α = {*x, y, z*} coordinate of atom *i*. For simplicity we assume in this work that the system contracts isotropically with a cubic unit cell of lattice constant *a* and ignore equilibrium distortions.

### 2.2. Projected Phonon Density of States

Although highly simplified phenomenological approaches like the Einstein and Debye models based on empirical data can be quite useful, they are generally inadequate to treat complex materials (Dimakis and Bunker, 1998; Poiarkova and Rehr, 1999). As an alternative which overcomes many of these limitations Poiarkova and Rehr (1999, 2001) introduced a method for EXAFS calculations in which the unit-normalized PDOS, projected onto a simple bond or arbitrary multiple-scattering path, is calculated from the imaginary part of the lattice dynamical Green's function. Vila et al. (2007) introduced an extension of the method to compute the PDOS projected onto atomic displacements ρ_{iα}(*a*, ω),

Here |*iα* > is a Lanczos seed vector representing a normalized, mass-weighted displacement of atom *i* along the Cartesian direction α such that the position of the center of mass of the cell is not changed, and **D**(*a*) is the *dynamical matrix* (DM) of force constants computed at lattice constant *a*,

where *u*_{jlα} is the displacement of atom *j* in unit cell *l* along the direction α, and *M*_{j} is the mass of atom *j*. As described below, **D**(*a*(*T*)) can be interpolated from **D**(*a*), which makes the force constants depend parametrically on the temperature and thus includes the dominant effects of anharmonicity. Efficient calculations of the lattice dynamical Green's function can be accomplished using a continued fraction representation, with parameters obtained with the iterative Lanczos algorithm (Deuflhard and Hohmann, 1995). This approach yields an efficient many-pole representation for the PDOS; since the pole locations can be interpreted as Gaussian quadrature points (Haydock, 1980), the method is well suited for accurate spectral integrations, as needed for example in the calculation of the VFE. The first iteration of the Lanczos algorithm corresponds to a *correlated Einstein* model for the |*iα* > displacement such that ρ_{iα}(*a*, ω) = δ(ω−ω_{E}) with an Einstein frequency ${\omega}_{E}^{2}=\langle i\alpha \left|\text{D}(a)\right|i\alpha \rangle $. Poiarkova et al. showed that a second tier continued fraction gave about 10% errors for EXAFS Debye-Waller factors, while Krappe and Rossner (2002) later showed that six iterations are needed to achieve convergence to within 1%. We have found (Vila et al., 2007) that 16 or more iterations may be needed to achieve the same precision for ρ_{iα}(*a*, ω). Nevertheless, this algorithm avoids the explicit calculations of phonon-spectra and hence permits highly efficient calculations.

### 2.3. Dynamical Matrix

The main computational bottleneck to using this approach for many complex materials of industrial or technological interest then lies in obtaining sufficiently accurate DMs. Although empirical or model potential estimates of the interatomic force constants are sometimes available, they are not useful in screening a broad range of potential new materials. Moreover, their temperature dependence usually limits their accuracy and generality. Nevertheless, we have previously shown (Vila et al., 2007) that *ab initio* DMs obtained from density functional theory (DFT) calculations provide accurate force constants that can be used to compute both EXAFS Debye-Waller factors and thermal expansion. In that study we restricted our attention to simple periodic systems using the methodology implemented in ABINIT (Gonze et al., 2002), described in detail in Gonze and Lee (1997). This method computes ${\stackrel{~}{D}}_{j\alpha {j}^{\prime}\beta}(a,\text{q})$, the dynamical matrix in reciprocal space:

and then Fourier transforms it back into real space. This approach is efficient for systems with simple unit cells of just a few atoms and high symmetry. In this paper we extend this method by using our recently developed parallelization scheme (Vila et al., Submitted) for the efficient computation of the DM in more complex unit cells based on first principles DFT calculations. This method computes the real space DM in Equation (6) directly using finite differences in supercells, taking advantage of the real space efficiency of codes like VASP (Kresse and Hafner, 1993, 1994; Kresse and Furthmüller, 1996a,b; Kresse and Joubert, 1999). This real space approach is particularly well suited for large unit cells with low symmetry such as those in supercell simulations of supported nanoparticles. In such cases the use of analytic methods such as DFPT results in very demanding calculations both in computing time and memory. A similar approach has recently been described for molecular systems (Liu et al., 2017) that suggest that the high availability of commodity processors can be used to take advantage of the parallelization of the finite differences. Similarly, our approach also relies on the parallel computation of complete rows of the DM, directly from the analytical DFT forces in VASP as a function of the lattice constant:

where *F*_{jlα} is the force on atom *j* in unit cell *l* along the direction α. For the case of ZrW_{2}O_{8} treated here we use centered finite differences with displacements of 0.015 Å, which provide sufficiently accurate force constants. The most efficient row partition of the DM depends on the configuration of the computing system used. Details for the partition used here are presented in section 2.5.

### 2.4. Other Considerations

As discussed in section 2.1, the temperature-dependence of the lattice constant *a*(*T*) is obtained by minimizing the free energy *A*(*a, T*) in Equation (2) with respect to *a* at a given temperature *T*. The temperature dependence can be obtained quite efficiently since it only requires the recomputation of the integral in Equation (3) given a PDOS ρ_{iα}(*a*, ω) at a certain lattice constant. The variation with lattice constant, however, requires the recomputation of ρ_{iα}(*a*, ω), and thus of **D**(*a*), which is the most computationally demanding step. In previous work (Vila et al., 2007) we solved this issue for simple materials by computing *A*(*a, T*) in a lattice constant grid that can be fitted or interpolated to find the minimum lattice constant as a function of temperature. To improve efficiency we took advantage of two properties of these simple systems: (i) we assumed that the shape of *A*(*a, T*) can be accurately described using a Morse potential form as a function of lattice constant; and (ii) we reduced the number of calculations of **D**(*a*) by taking advantage of the nearly linear behavior of **D**(*a*) as a function of lattice constant near equilibrium,

where **D**_{2} is the constant harmonic DM at *a*_{0}, the lattice constant of minimum internal energy, and **D**_{3} is the constant matrix of anharmonic constants. This parametrization assumes that anharmonic contributions of order higher than cubic are small, and allows us to compute the DM for only two lattice constants that can be used to determine **D**_{2} and **D**_{3}. This approach typically cuts the computational time by a factor of about 2/3 or more. Given the complexity of some of the systems in the present study and the importance of having reliable efficiency enhancements for computational screening of materials, here we have opted to verify that this approximation is still valid by computing the DMs over a more complete grid of lattice constants. We have also changed the form used to fit *A*(*a, T*) to a simple polynomial.

### 2.5. Computational Details

All structural optimizations and DM calculations were performed with VASP (Kresse and Hafner, 1993, 1994; Kresse and Furthmüller, 1996a,b; Kresse and Joubert, 1999) using PAW potentials (Kresse and Joubert, 1999) and the PBEsol exchange-correlation functional (Perdew et al., 2008). An 8 × 8 × 8 *k*-point grid was used in all calculations, which was sufficient to achieve converged results. To reduce the effect of Pulay stress, the planewave energy cutoff was set at 350 eV. The finite difference force calculations used a three point centered stencil with 0.015 Å displacements. To obtain the variation of the DM with respect to lattice constant, we start by optimizing the size of unit cell and the reduced positions of the atoms in it while preserving its symmetry. This produces an optimized unit cell expanded by 1% from the experimental value. Given that the reduced positions do not change with volume until the phase transition at about 430 K (Evans et al., 1999), we freeze them and expand/contract around the optimal lattice size. We compute the DM at expansions of 0.0, 0.5, 1.0, 1.5, and 2.0% with respect to the experimental lattice constant, which correspond to contractions of 0.5 and 1.0%, and expansions of 0.5 and 1.0% with respect to the optimal lattice contant. This procedure provides an adequate number of points for the interpolation of both the energies and components of the DM. For this material and taking into consideration the computing system being used, we partitioned the DM into 44 independent blocks, one per atom in the simulation cell, each using 512 processors for a total of 22,528. This procedure reduced the computing time from 33 h/core to 0.8 h/core. Equivalent simulations using DFPT in ABINIT would require on the order of 30–40 h/core, with a maximum number of processors of about 200–300.

## 3. Results and Discussion

### 3.1. Structure

Figure 1 shows the structure of the unit cell of ZrW_{2}O_{8} used here. The cell is composed of four ZrO_{6} octahedra linked to eight WO_{4} tetrahedra. The ZrO_{6} units are all equivalent, with three short and three long Zr–O bonds, each connecting to a WO_{4} unit (Figure 2A). The WO_{4} units have three long W–O bonds, connecting to the ZrO_{6} units, and a short one that is essentially free (Figure 2B). Unlike the ZrO_{6} units, there are two types of WO_{4} units depending on the direction of the free W–O bond: Four of them point toward a ZrO_{6} unit and have a truly free W–O bond, while the other four point to another WO_{4} unit and have a more restricted W–O bond. These two types of tetrahedra are paired (Figure 2B). Here we label the different types of atoms in the cell as follows: O_{SL} and O_{LS} are O atoms bridging WO_{4} and ZrO_{6} units, with short W–O and long Zr–O bonds, and vice versa. They correspond, respectively, to the O_{2} and O_{1} atoms in the more traditional notation (Mary et al., 1996). The O_{FF} and O_{FR} are O atoms in truly free and somewhat restricted W–O bonds, respectively. They correspond to the O_{4} and O_{3} atoms in the traditional notation (Mary et al., 1996). Finally, the W atoms bonded to O_{FF} and O_{FR} center are similarly labeled W_{FF} and W_{FR}, and correspond to the W_{1} and W_{2} atoms in the traditional notation (Mary et al., 1996), respectively. Table 1 presents a comparison between the experimental (Auray et al., 1995) and optimized structural parameters obtained in the PBEsol optimization. The changes in internal structure are minimal except for the overall expansion of the cell of 1% predicted by PBEsol.

**Figure 1**. Structure of the ZrW_{2}O_{8} cell. The cell is composed of four equivalent Zr atoms (magenta), eight W atoms, of type W_{FF} (blue) and W_{FR} (green), and 32 O atoms of type O_{LS} (yellow), O_{SL} (red), O_{FR} (orange) and O_{FF} (black). Examples of the different types of O and W atoms are highlighted, and the coordination tetrahedra and octahedra are also included for clarity. The W_{FF}, W_{FR}, O_{LS}, O_{SL}, O_{FR}, and O_{FF} atoms correspond, respectively, to the W_{1}, W_{2}, O_{1}, O_{2}, O_{3}, and O_{4} in the traditional notation (Mary et al., 1996).

**Figure 2**. Local structure of ZrW_{2}O_{8} around the Zr octahedra **(A)** and W tetrahedra **(B)**. The color key is the same as in Figure 1. Examples of the different types of O and W atoms are highlighted. The bottom figure shows the pairing of the two types of W atoms. The W_{FF}, W_{FR}, O_{LS}, O_{SL}, O_{FR}, and O_{FF} atoms correspond, respectively, to the W_{1}, W_{2}, O_{1}, O_{2}, O_{3}, and O_{4} in the traditional notation (Mary et al., 1996).

**Table 1**. Experimental and theoretical structural parameters for the ZrW_{2}O_{8} cell used in this work.

### 3.2. Bond Force Constants

Figure 3 shows the lattice constant dependence of the mean force constants for the parallel stretch and perpendicular scissoring displacements for the different bond types in ZrW_{2}O_{8} (The numerical data for this and all other plots in this work can be found in Supplementary Data Sheet 1). The parallel force constants were obtained by rotating the 6 × 6 ${(jl,{j}^{\prime}{l}^{\prime})}_{\alpha ,\beta}$ blocks of the DM formed by *jl* and *j*′*l*′ near-neighbor atoms into a local coordinate system along the bond. The scissors force constants were obtained similarly by rotating the 3 × 3 (*jl, jl*)_{α, β} block for O atom *jl* into coordinates parallel and perpendicular to the W–O bond. For the W_{FF}–O_{FF} and W_{FR}–O_{FR} bonds, motion of the O atoms orthogonal to the bonds is nearly equivalent in all directions. For the W_{FF}–O_{LS} and W_{FR}–O_{SL} bonds, we average over the in-plane and out-of-plane directions defined by the W–O–Zr plane. These directions are also nearly equivalent, with anisotropies of at most 20%.

**Figure 3**. Mean force constants for the parallel OW and OZr stretch **(Top)** and perpendicular OW scissoring **(Bottom)** displacements in ZrW_{2}O_{8}. The perpendicular force constants are averaged over two directions in the plane orthogonal to the W–O bond for the W_{FF}–O_{FF} and W_{FR}–O_{FR} bonds, and over the in-plane and out-of-plane directions for the W_{FF}–O_{LS} and W_{FR}–O_{SL} bonds. See text for further details.

We find that of all the elements of the DM, only those associated with the O atoms exhibit the positive anharmonicity essential to observe NTE. Indeed, only the deviations perpendicular to the bonds have positive slope (Figure 3, bottom). The force constants associated with the perpendicular O distortions are also much weaker than the parallel ones, and are therefore activated at much lower temperature. The strength of the parallel bond distortions follows closely the trends in the bond distances shown in Table 1, with the shorter W_{FF}–O_{FF} and W_{FR}–O_{FR} bonds having the largest force constants, followed by the slightly longer W_{FF}–O_{LS} and W_{FR}–O_{SL} bonds. As discussed in section 2.4, the components of the dynamical matrix should vary mostly linearly with cell expansion in simple systems. Given the large NTE observed for ZrW_{2}O_{8}, which implies the existence of significant anharmonicity, it is worth examining the accuracy of the linear approximation. Figure 3 shows that the linearity of the bond-parallel force constants is excellent for small displacements, while the transverse ones show a slight curvature. The mean absolute errors from a linear approximation are ± 3 N/m for the strong parallel displacements and only ± 0.4 N/m for the transverse ones. Consequently, the linear approximation can be used to estimate the the lattice dependence of the dynamical matrix to high accuracy.

### 3.3. Helmholtz and Vibrational Free Energies

Figure 4 (top) shows the variation with lattice constant and temperature of the internal energy *U*(*a*), total vibrational free energy *F*(*a, T*), and Helmholtz free energy *A*(*a, T*), together with the optimal lattice constant as a function of temperature obtained by fitting *A*(*a, T*) to a 5th-order polynomial and minimizing numerically. The NTE clearly arises from the lower relative VFE at smaller lattice constants. To help understand the origin of this behavior, the bottom of Figure 4 also shows a decomposition of the VFE at low and high temperature, obtained by partitioning the sum over atoms in Equation (3) into its O, W, and Zr components. This decomposition reveals that the W and Zr contributions to the VFE are virtually constant with variations in lattice constant; consequently the variation of total VFE closely follows the variation of only the O component. This further substantiates our finding that the perpendicular O distortions are the primary origin of NTE in ZrW_{2}O_{8}.

**Figure 4. (Top)** Internal energy *U*(*a*), total vibrational free energy *F*(*a, T*) and Helmholtz free energy *A*(*a, T*) as a function of lattice constant and temperature. The potentials are interpolated using a 5th-order polynomial. The black dots indicate the optimal lattice constant as a function of temperature. **(Bottom)** Total *F*(*a, T*) and contributions from the O, W, and Zr sites.

### 3.4. Negative Thermal Expansion

The relative thermal expansion obtained by minimizing *A*(*a, T*) with respect to lattice constant is shown in Figure 5. The agreement with experiment is clearly very good except at very low temperatures, where the theoretical curve shows the expected flattening due to zero-point effects, which does not appear to be pronounced in the experiment (Evans et al., 1996). Figure 5 also shows the NTE that would be obtained by only considering the components of the VFE arising from the O, W, and Zr lattices. As expected from the previous discussion, the NTE arises exclusively from the O component. The Zr and W contributions have a minimal effect on the system, in contrast with the hypothesis (Mary et al., 1996; Gupta et al., 2013; Sanson, 2014) that the bending of the Zr–O–W bonds causes the NTE by pulling the Zr and W atoms together.

**Figure 5**. Total, O, W, and Zr components, and results from our two-pole model (see text) of relative thermal expansion in ZrW_{2}O_{8} as a function of temperature. The experimental results are taken from Evans et al. (1996).

### 3.5. Total and Projected Phonon Densities of States

Figure 6 shows a comparison between the theoretical total density of states at 0 K and experiment (Ernst et al., 1998) obtained from inelastic neutron scattering at 300 K. The PDOS can be roughly divided into three regions: high frequency range (above 20 THz), associated with bond stretching modes, middle range (between 5 and 13 THz), corresponding mostly to bond bending modes, and low frequency range (below 5 THz), which are thought to correspond to librational and translational modes (Chaplot, 2005). The overall agreement is excellent; the small discrepancy at the lowest frequencies around 1 THz, where the Lanczos approach yields a single peak, is likely due to the size of the 44-atom unit cell used in the calculations.

**Figure 6**. Comparison of the theoretical (0 K) and experimental (Ernst et al., 1998) (300 K) total phonon density of states (DOS) in ZrW_{2}O_{8}. The experimental result has been vertically shifted for clarity.

To understand the physical origin of the NTE more completely, we consider the behavior of the projected phonon density of states, or PDOS. Figure 7 shows the PDOS projected onto the O, W, and Zr sites. The high frequency PDOS only has significant weight on the O sites. The W and Zr sites have most of their weight in the low and medium regions, respectively. The variation of the PDOS with lattice constant depends greatly on the projection site as shown by the Grüneisen parameters γ = *d*lnω/*d*ln*V* of the different sites and frequency ranges: The middle range is essentially independent of lattice constant with a γ = − 0.3 ± 0.3, except for the Zr sites which show a very small decrease of the frequencies with expansion and γ = 2.0. This positive γ results in the very small tendency to normal behavior seen in the Zr component of the thermal expansion shown in Figure 5. The high frequency range, found only on the O sites, also contributes to positive thermal expansion with a γ = 1.9 ± 0.2. These high frequencies result from the strong parallel force constants discussed in section 3.2, as demonstrated in Figure 8, which shows that the DOS projected onto the longitudinal (i.e., along the O–W or Zr–O–W direction) motion of the O atoms only has weight in the high frequency range. These stretch modes become weaker as the system expands, and would result in net positive thermal expansion if not countered by the low frequency region. In fact, if the high frequency range is neglected, the relative NTE at 0 K is overestimated by about 30%. The low frequency region exhibits a large negative γ = − 11 ± 2 for the O sites, and about γ = − 3 ± 1 for the W sites, in agreement with our finding that the W sites have a small contribution to the NTE, while the O sites account for most of it. For the overall low frequency range we estimate the same γ of − 8 ± 4, as in previous theoretical estimates, (Gupta et al., 2013) which is in reasonable agreement with experimental estimates ranging from − 7 to − 20 (Hancock et al., 2004). Figure 8 also shows the projection onto the transverse motion of the O atoms, which is localized in the medium and low frequency ranges, as expected for bending and librational modes, respectively.

**Figure 7**. Variation of the projected phonon density of states (PDOS) as a function of lattice constant in ZrW_{2}O_{8}. The partial PDOS projected onto O , W, and Zr sites are shown in the top, middle and lower panels, respectively. The curves for different lattice constants have been vertically shifted for clarity.

**Figure 8**. Variation of the O-sites phonon density of states (PDOS) as a function of lattice constant in ZrW_{2}O_{8}, projected onto longitudinal (along the O–W or Zr–O–W direction), and transverse in-plane (ip) and out of plane (op) displacements. The curves for different lattice constants have been vertically shifted for clarity.

To explore the role of the rigid WO_{4} units on the NTE, we have also projected the PDOS onto rotational linear combinations of the transverse in-plane and out-of-plane distortions described above, with axes along the different W–O bonds (Figures 9, 10). These projections eliminate all the internal O–W–O bending contributions seen in the 5–12 THz range of Figure 8, with only features associated with O–Zr–O bending remaining in that range. All these PDOS features show the expected negative Grüneisen parameters γ, and point to the need for a mixed model that is neither purely rigid (Hammonds et al., 1996; Pryde et al., 1996, 1997; Hancock et al., 2004; Tucker et al., 2005, 2007) nor fully distorted (Cao et al., 2002, 2003; Bridges et al., 2014) but one in which the NTE arises from the rigid motion of WO_{4} units and the internal distortion of the ZrO_{6} ones.

**Figure 9**. Variation of the ZrW_{2}O_{8} phonon density of states (PDOS) projected onto rotational coordinates around the W_{FF}–O_{FF} and W_{FF}–O_{LS} axes, as a function of the lattice constant. The curves for different lattice constants have been vertically shifted for clarity.

**Figure 10**. Variation of the ZrW_{2}O_{8} phonon density of states (PDOS) projected onto rotational coordinates around the W_{FR}–O_{FR} and W_{FR}–O_{SL} axes, as a function of lattice constant. The curves for different lattice constants have been vertically shifted for clarity.

### 3.6. Simple NTE Model

Based on the results discussed in the previous sections, we are led to formulate a simple anharmonic two-frequency Einstein model to describe the NTE. This model is based on approximations to the different contributions to *A*(*a, T*) and *F*(*a, T*) in Equations (2) and (3). Near equilibrium the internal energy *U*(*a*) can be approximated accurately with a harmonic potential $U(a)=(1/2){K}_{U}{(\Delta a)}^{2}$, and Δ*a* ≡ *a* − *a*_{0}. The VFE *F*(*a, T*) can be simplified by analyzing the different contributions arising from the PDOS of Figure 7: We have found that the total O-PDOS can be reduced to a very good approximation to just two constant weight poles representing the active regions of the O-PDOS: one for the low frequency (*L*) region below 5.4 THz with a negative Grüneisen coefficient, and a second high frequency (*H*) pole for the modes above 20 THz with a positive Grüneisen coefficient. The intermediate region of the O, Zr, and W PDOS do not contribute to the NTE since they are nearly constant with Δ*a*. The contribution from the *H* pole can be further simplified by noticing that the *H* modes only contribute to *F*_{0}(*a*) since at the temperatures relevant here the *H* modes are inactive. We can also take advantage of the linear behavior of the force constants as a function of Δ*a* as discussed in section 2.4: For a given pole ν with effective force constant *k*_{ν}, we can write *k*_{ν} = *k*_{2ν} + 6*k*_{3ν}Δ*a* (Frenkel and Rehr, 1993; Vila et al., 2007), where *k*_{2ν} and *k*_{3ν} are, respectively, the effective quasi-harmonic and anharmonic force constants. Thus, the frequency ω_{ν} associated with this pole can be approximated as

where we define ${\omega}_{\nu}^{0}=\sqrt{{k}_{2\nu}/{\mu}_{\nu}}$, ${\omega}_{\nu}^{\prime}=3{k}_{2\nu}/{\omega}_{\nu}^{0}{k}_{3\nu}$, and μ_{ν} is the reduced mass associated with this pole. The linear behavior of the pole frequencies around *a*_{0} can be seen in the PDOS of Figure 7, from which the ${\omega}_{\nu}^{0}$ and ${\omega}_{\nu}^{\prime}$ parameters can be obtained. The poles weights *w*_{L} and *w*_{H} are taken at Δ*a* = 0. With these approximations the Helmholtz free energy from Equation (2) can be written as:

where ${F}_{0}^{\nu}(a)$ and ${F}_{T}^{\nu}(a,T)$ are, respectively, the zero-point and temperature components of *F*(*a, T*) for ν = {*L, H*}. We can now find *a*(*T*) by applying the minimum condition in Equation (1) to the above *A*(*a, T*):

The last term in this expression can be further simplified by noticing in Figure 4 (bottom) the nearly linear behavior of *F*(*a, T*) around *a*_{0}, and expanding ${F}_{T}^{L}(a,T)$ to first order in *a*. After this, we can calculate *a*(*T*) simply as:

using the internal energy minimum and constant *a*_{0} = 9.246 Å and *K*_{U} = 85.82 eV/Å^{2}, respectively, weights *w*_{L} = 21.1 and *w*_{H} = 30.0, harmonic frequency ${\omega}_{L}^{0}=20.31$ THz, and effective anharmonic constants ${\omega}_{L}^{\prime}=81.90$ THz/Å and ${\omega}_{H}^{\prime}=-104.39$ THz/Å. Our results are presented in Figure 5; they clearly show that this simplified two-pole Einstein model yields a quantitative description of the observed NTE. Of these contributions the *L* pole dominates but overestimates the NTE, while the *H* pole provides a 30% correction.

## 4. Conclusions

We have performed calculations of the thermal expansion of ZrW_{2}O_{8} using an efficient approach for obtaining the variation of the dynamical matrix as a function of the lattice constant. The anomalous NTE arises almost exclusively from the transverse contributions of the O-atoms, as demonstrated by a partition of the vibrational free energy over the O, W, and Zr sites. This behavior results from the interplay between the normal positive Grüneisen parameter and hence PTE of the O-projected high frequency modes, and the anomalous negative-Grüneisen parameter driving the NTE from the low frequency modes. The low-frequency modes overestimate the effect of NTE, while the high-frequency modes provide a 30% correction. These low frequency modes are associated with distortions at the O sites transverse to the W–O or W–O–Zr bonds. More precisely, the rotational linear combination of these transverse distortions results in nearly rigid rotations of the WO_{4} tetrahedra and internal distortions of the ZrO_{6} octahedra. The mixed character of the low frequency modes is consistent with the variety of models previously proposed to account for the NTE. We have simplified this complex behavior by proposing a two-frequency anharmonic Einstein model of the O-PDOS which captures both the dominant negative- and weak positive-contributions to thermal expansion.

## Author Contributions

FV co-developed the theory, was the main developer of the specialized software used in the analysis, and performed most of the simulations. SH contributed to the simulations and helped in verifying the soundness of the software. JR co-developed the theory presented here.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

This work was supported by U.S. Department of Energy, Office of Science CSGBD Grant No. DE-FG02-03ER15476, with computational support from NERSC, a DOE Office of Science User Facility, under Contract No. DE-AC02-05CH11231. We thank Frank Bridges and Joshua Kas for many helpful discussions.

## Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2018.00296/full#supplementary-material

**Supplementary Data Sheet 1**. Numerical data for all plots in this work.

## References

Allen, R. E., and De Wette, F. W. (1969). Calculation of dynamical surface properties of noble-gas crystals. I. The quasiharmonic approximation. *Phys. Rev.* 179, 873–886. doi: 10.1103/PhysRev.179.873

Auray, M., Quarton, M., and Leblanc, M. (1995). Zirconium tungstate. *Acta Crystallogr. Sect. C* 51, 2210–2213. doi: 10.1107/S0108270195001545

Baroni, S., de Gironcoli, S., Dal Corso, A., and Giannozzi, P. (2001). Phonons and related crystal properties from density-functional perturbation theory. *Rev. Mod. Phys.* 73, 515–562. doi: 10.1103/RevModPhys.73.515

Baroni, S., Giannozzi, P., and Isaev, E. (2010). Density-functional perturbation theory for quasi-harmonic calculations. *Rev. Mineral. Geochem.* 71:39. doi: 10.2138/rmg.2010.71.3

Boyer, L. L. (1979). Calculation of thermal expansion, compressiblity, an melting in alkali halides: NaCl and KCl. *Phys. Rev. Lett.* 42, 584–587. doi: 10.1103/PhysRevLett.42.584

Brandenburg, J. G., Potticary, J., Sparkes, H. A., Price, S. L., and Hall, S. R. (2017). Thermal expansion of carbamazepine: systematic crystallographic measurements challenge quantum chemical calculations. *J. Phys. Chem. Lett.* 8, 4319–4324. doi: 10.1021/acs.jpclett.7b01944

Bridges, F., Keiber, T., Juhas, P., Billinge, S. J. L., Sutton, L., Wilde, J., et al. (2014). Local vibrations and negative thermal expansion in ZrW_{2}O_{8}. *Phys. Rev. Lett.* 112:045505. doi: 10.1103/PhysRevLett.112.045505

Cao, D., Bridges, F., Kowach, G. R., and Ramirez, A. P. (2002). Frustrated soft modes and negative thermal expansion in ZrW_{2}O_{8}. *Phys. Rev. Lett.* 89:215902. doi: 10.1103/PhysRevLett.89.215902

Cao, D., Bridges, F., Kowach, G. R., and Ramirez, A. P. (2003). Correlated atomic motions in the negative thermal expansion material ZrW_{2}O_{8}: a local structure study. *Phys. Rev. B* 68:014303. doi: 10.1103/PhysRevB.68.014303

Chaplot, S. L. (2005). Negative thermal expansion in ZrW_{2}O_{8}: do we give up the concept of normal mode? *Curr. Sci.* 88, 347–349.

Dimakis, N., and Bunker, G. (1998). Ab initio single- and multiple-scattering EXAFS Debye-Waller factors: Raman and infrared data *Phys. Rev. B* 58:2467. doi: 10.1103/PhysRevB.58.2467

Erba, A. (2014). On combining temperature and pressure effects on structural properties of crystals with standard ab initio techniques. *J. Chem. Phys.* 141:124115. doi: 10.1063/1.4896228

Erba, A., Maul, J., and Civalleri, B. (2016). Thermal properties of molecular crystals through dispersion-corrected quasi-harmonic ab initio calculations: the case of urea. *Chem. Commun.* 52, 1820–1823. doi: 10.1039/C5CC08982D

Erba, A., Shahrokhi, M., Moradian, R., and Dovesi, R. (2015). On how differently the quasi-harmonic approximation works for two isostructural crystals: thermal properties of periclase and lime. *J. Chem. Phys.* 142:044114. doi: 10.1063/1.4906422

Ernst, G., Broholm, C., Kowach, G., and Ramirez, A. (1998). Phonon density of states and negative thermal expansion in ZrW_{2}O_{8}. *Nature* 396, 147–149. doi: 10.1038/24115

Evans, J. S. O., David, W. I. F., and Sleight, A. W. (1999). Structural investigation of the negative-thermal-expansion material ZrW_{2}O_{8}. *Acta Crystallogr. Sect. B* 55, 333–340. doi: 10.1107/S0108768198016966

Evans, J. S. O., Mary, T. A., Vogt, T., Subramanian, M. A., and Sleight, A. W. (1996). Negative thermal expansion in ZrW_{2}O_{8} and HFW_{2}O_{8}. *Chem. Mater.* 8, 2809–2823. doi: 10.1021/cm9602959

Frenkel, A. I., and Rehr, J. J. (1993). Thermal expansion and x-ray-absorption fine-structure cumulants. *Phys. Rev. B* 48, 585–588. doi: 10.1103/PhysRevB.48.585

Gonze, X., Beuken, J.-M., Caracas, R., Detraux, F., Fuchs, M., Rignanese, G.-M., et al. (2002). First-principles computation of material properties: the ABINIT software project *Comput. Mater. Sci.* 25:478. doi: 10.1016/S0927-0256(02)00325-7

Gonze, X., and Lee, C. (1997). Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory. *Phys. Rev. B* 55:10355. doi: 10.1103/PhysRevB.55.10355

Gupta, M. K., Mittal, R., and Chaplot, S. L. (2013). Negative thermal expansion in cubic ZrW_{2}O_{8}: role of phonons in the entire Brilloiun zone from ab initio calculations. *Phys. Rev. B* 88:014303. doi: 10.1103/PhysRevB.88.014303

Hammonds, K. D., Dove, M. T., Giddy, A. P., Heine, V., and Winkler, B. (1996). Rigid-unit phonon modes and structural phase transitions in framework silicates. *Amer. Mineral.* 81, 1057–1079. doi: 10.2138/am-1996-9-1003

Hancock, J. N., Turpen, C., Schlesinger, Z., Kowach, G. R., and Ramirez, A. P. (2004). Unusual low-energy phonon dynamics in the negative thermal expansion compound ZrW_{2}O_{8}. *Phys. Rev. Lett.* 93:225501. doi: 10.1103/PhysRevLett.93.225501

Haydock, R. (1980). “The recursive solution of the schrodinger equation,” in *Solid State Physics*, Vol. 35, eds H. Ehrenreich, F. Seitz, and D. D. Turnbull (New York, NY: Academic), 215.

Jain, A., Ong, S. P., Hautier, G., Chen, W., Richards, W. D., Dacek, S., et al. (2013). The materials project: a materials genome approach to accelerating materials innovation. *APL Mater.* 1:011002. doi: 10.1063/1.4812323

Krappe, H. J., and Rossner, H. H. (2002). Bayes-Turchin approach to x-ray absorption fine structure data analysis. *Phys. Rev. B* 66:184303. doi: 10.1103/PhysRevB.66.184303

Kresse, G., and Furthmüller, J. (1996a). Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. *Comput. Mater. Sci.* 6, 15–50. doi: 10.1016/0927-0256(96)00008-0

Kresse, G., and Furthmüller, J. (1996b). Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. *Phys. Rev. B* 54, 11169–11186. doi: 10.1103/PhysRevB.54.11169

Kresse, G., and Hafner, J. (1993). Ab initio. *Phys. Rev. B* 47, 558–561. doi: 10.1103/PhysRevB.47.558

Kresse, G., and Hafner, J. (1994). Ab initio. *Phys. Rev. B* 49, 14251–14269. doi: 10.1103/PhysRevB.49.14251

Kresse, G., and Joubert, J. (1999). From ultrasoft pseudopotentials to the projector augmented-wave method. *Phys. Rev. B* 59:1758. doi: 10.1103/PhysRevB.59.1758

Lee, C., and Gonze, X. (1995). Ab initio calculation of the thermodynamic properties and atomic temperature factors of *SIO*_{2} α−quartz and stishovite. *Phys. Rev. B* 51, 8610–8613. doi: 10.1103/PhysRevB.51.8610

Liu, K. Y., Jie, L., and Herbert, J. M. (2017). Accuracy of finite? Difference harmonic frequencies in density functional theory. *J. Comput. Chem.* 38, 1678–1684. doi: 10.1002/jcc.24811

Mary, T. A., Evans, J. S. O., Vogt, T., and Sleight, A. W. (1996). Negative thermal expansion from 0.3 to 1050 Kelvin in ZrW_{2}O_{8}. *Science* 272, 90–92. doi: 10.1126/science.272.5258.90

Perdew, J., Ruzsinszky, A., Csonka, G., Vydrov, O., Scuseria, G., Constantin, L., et al. (2008). Restoring the density-gradient expansion for exchange in solids and surfaces. *Phys. Rev. Lett.* 100:136406. doi: 10.1103/PhysRevLett.100.136406

Poiarkova, A., and Rehr, J. J. (1999). Multiple-scattering x-ray-absorption fine-structure Debye-Waller factor calculations. *Phys. Rev. B* 59:948. doi: 10.1103/PhysRevB.59.948

Poiarkova, A., and Rehr, J. J. (2001). Recursion method for multiple-scattering XAFS Debye-Waller factors *J. Synchrotron Radiat.* 8:313. doi: 10.1107/S0909049599001685

Pryde, A. K., Hammonds, K. D., Dove, M. T., Heine, V., Gale, J. D., and Warren, M. C. (1997). Rigid unit modes and the negative thermal expansion in ZrW_{2}O_{8}. *Phase Transit.* 61, 141–153. doi: 10.1080/01411599708223734

Pryde, A. K. A., Hammonds, K. D., Dove, M. T., Heine, V., Gale, J. D., and Warren, M. C. (1996). *J. Phys.* 8:10973.

Rignanese, G., Michenaud, J., and Gonze, X. (1996). Ab initio study of the volume dependence of dynamical and thermodynamical properties of silicon. *Phys. Rev. B* 53, 4488–4497. doi: 10.1103/PhysRevB.53.4488

Sanson, A. (2014). Toward an understanding of the local origin of negative thermal expansion in ZrW_{2}O_{8}: limits and inconsistencies of the tent and rigid unit mode models. *Chem. Mater.* 26, 3716–3720. doi: 10.1021/cm501107w

Tucker, M. G., Goodwin, A. L., Dove, M. T., Keen, D. A., Wells, S. A., and Evans, J. S. O. (2005). Negative thermal expansion in ZrW_{2}O_{8}: mechanisms, rigid unit modes, and neutron total scattering. *Phys. Rev. Lett.* 95:255501. doi: 10.1103/PhysRevLett.95.255501

Tucker, M. G., Keen, D. A., Evans, J. S. O., and Dove, M. T. (2007). Local structure in ZrW_{2}O_{8} from neutron total scattering. *J. Phys. Condens. Matt.* 19:335215. doi: 10.1088/0953-8984/19/33/335215

Keywords: zirconium tungstate, NTE, DFT, quasi-harmonic approximation, phonon density of states

Citation: Vila FD, Hayashi ST and Rehr JJ (2018) Efficient Calculation of the Negative Thermal Expansion in ZrW_{2}O_{8}. *Front. Chem*. 6:296. doi: 10.3389/fchem.2018.00296

Received: 27 April 2018; Accepted: 26 June 2018;

Published: 30 July 2018.

Edited by:

Andrea Sanson, Università Degli Studi di Padova, ItalyReviewed by:

Simone Sanna, Justus Liebig Universität Gießen, GermanyAlessandro Erba, Università degli Studi di Torino, Italy

Copyright © 2018 Vila, Hayashi and Rehr. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Fernando D. Vila, fdv@uw.edu