research papers
A selfconsistent approach to describe unitcellparameter and volume variations with pressure and temperature
^{a}IGG, CNR, Via G. Gradenigo, 6, Padova, Padova I35131, Italy, ^{b}Mainz Institute of Multiscale Modeling and Institute of Geosciences, JohannesGutenberg University of Mainz, J.J.BecherWeg 21, Mainz D55128, Germany, ^{c}Departamento de Física, Instituto Universitario de Estudios Avanzados en Física Atómica, Molecular y Fotónica (IUDEA), MALTA Consolider Team, Universidad de La Laguna, La Laguna, Tenerife E38204, Spain, and ^{d}Department of Earth and Environmental Sciences, University of Pavia, Via A. Ferrata 1, Pavia I27100, Italy
^{*}Correspondence email: rossjohnangel@gmail.com
A method for the selfconsistent description of the large variations of unitcell parameters of crystals with pressure and temperature is presented. It employs linearized versions of equations of state (EoSs) together with constraints to ensure internal consistency. The use of polynomial functions to describe the variation of the unitcell angles in monoclinic and triclinic crystals is compared with the method of deriving them from linearized EoSs for d spacings. The methods have been implemented in the CrysFML Fortran subroutine library. The unitcell parameters and the compressibility and tensors of crystals can be calculated from the linearized EoSs in an internally consistent manner in a new utility in the EosFit7c program, which is available as freeware at https://www.rossangel.net.
Keywords: equations of state; unitcell parameters; EosFit; pressure.
1. Introduction
The measurement of the response of the unitcell parameters and therefore the volume of crystals to hydrostatic pressure (P) and temperature (T) provides fundamental information about the nature and anisotropy of the bonding within the structures of crystals. The variation of the unitcell volume, molar volume and density of a material with pressure and temperature is described by its equation of state (EoS). The formulae that relate the volume or density of a material to the applied pressure are based on various assumptions. They include assumed interatomic potentials and structural geometries, or an assumed relationship between parameters and pressure (e.g. the Murnaghan and Tait EoSs), or an assumed relationship between the strain arising from compression and the free energy of the solid (e.g. the Birch–Murnaghan EoS). Thermal EoSs include various purely parametric forms and those such as the Mie–Grünesien–Debye EoS that involve assumptions about the phonon of the solid and its contribution to the and thus the thermal pressure. A full review of EoSs and these issues has been given by Anderson (1995).
An EoS for volume (and therefore density) describes isotropic properties. Such EoSs do not describe the anisotropy of the response of crystal structures to pressure or temperature, which is described by the secondrank tensors of thermal expansivity and compressibility. While the variations of individual elastic moduli of a crystal (a fourthrank tensor) with P and T are related to the Helmholtz free energy (e.g. Stixrude & LithgowBertelloni, 2005), the expression involves derivatives of the individual elastic moduli that are often not available, while the equations cannot easily be reduced to describe the anisotropic and compressibility under hydrostatic pressure.
Therefore, the variations of individual unitcell parameters of crystals as P and T are changed have been described by using linearized equations of state in which the individual cell parameters are cubed and then fitted using the volume EoS (e.g. Angel, 2000). When applied to unitcell parameters, these have been called `axial EoS'. Alternatively, the equations for or compressibility can be converted to equivalent linear forms to fit unitcell data (e.g. Kroll et al., 2014; Murshed et al., 2015). These methods can be extended to any direction within the In most cases they provide an accurate description of the variation of the individual unitcell parameters with P and/or T, and yield linear compressibilities and moduli that agree with those determined by direct measurements of the elastic tensors. The advantage of these approaches is that the underlying theory of volume EoSs is well developed and their limitations in parameter and P and T space are well understood (e.g. Holzapfel, 2001; Anderson, 1995; Angel et al., 2019), so linearized EoSs can be extrapolated beyond the range of data with more confidence and justification than simple polynomial functions of the unitcell parameters with P and T.
However, the volume and the unitcell parameters of a crystal are not independent of one another, so the use of independent equations for the variation of the unitcell parameters and the volume introduces additional false P and T of a and its volume with independent EoSs is not physically consistent.
into the description of the behaviour of the crystal. Furthermore, the relationship between the volume and the unitcell parameters imposes specific constraints that relate the bulk modulus, , and its derivatives to the axial linear moduli, , and their derivatives. As we prove below, except for cubic crystals, as pressure is increased the different axial moduli increase at different rates, and the symmetry constraints relating the bulk modulus and its derivatives to the linear moduli and their derivatives are violated if conventional volume and axial EoSs are used to describe the variation of the moduli with pressure. The same problem occurs with the coefficients of the volume and the unitcell parameters. Therefore, describing the anisotropic changes withThere is an additional but distinct difficulty that arises when the crystal has monoclinic or triclinic symmetry. The use of an axial EoS does not address how to describe the variation of unitcell angles with P and T, or composition. This can be achieved by using the strain tensor and its derivative tensors of compressibility and thermal expansivity (Nye, 1957; Ohashi & Burnham, 1973). These define the instantaneous variation of all of the unitcell parameters of a crystal, including the unitcell angles. The linearized EoS defines the components of these property tensors that correspond to the axes of the but there is no independent underlying theory that defines how the other tensor components vary with P and T in monoclinic and triclinic crystals. This means that a strain tensor analysis cannot be used to extrapolate the behaviour of the crystal to P and T ranges beyond that of the data. Furthermore, when two unitcell determinations are made at both different P and different T, the strain is defined but the partitioning of the strain into T and Pinduced components is not, so the data cannot be interpreted in terms of compression and And when the unitcell angles of monoclinic and triclinic crystals change significantly from one measurement to another, the strain tensors do not provide a unique description of the unitcellparameter variation or of the directions of the principal axes of the strain, which include the directions of greatest and least strain (e.g. Paufler & Weber, 1999; Knight, 2010; Langreiter & Kahlenberg, 2015). This is important when trying to relate the directions of greatest or least compressibility or to the crystal architecture and thus the bonding within the crystal structure.
All of these issues have recently become more important in mineralogy as hostinclusion piezobarometry has been developed. This technique uses the stress and the strain measured in an inclusion crystal trapped inside a host mineral to determine the P and T under which the inclusion was trapped. In its simplest isotropic approximation only the volume equations of state of the minerals are required (Angel, Mazzucchelli et al., 2014). More recent developments of the method for anisotropic phases (Alvaro et al., 2020; Mazzucchelli et al., 2019; Gonzalez et al., 2021) require an accurate, precise and internally consistent description of how the unitcell parameters of both the host and inclusion phases change with pressure and temperature. If the equations for the unitcell parameters of a phase do not result in exactly the same volume as that calculated from its volume EoS, then the small discrepancies can propagate into significant errors in the calculation of entrapment conditions.
In this paper we present a simple phenomenological approach to describe consistently the variation of the unitcell parameters of crystals with P and T. It is based on using axial EoSs to describe the cellparameter variations, but with constraints to ensure that full internal consistency is maintained between the predicted unitcell parameters and volume. We then use the method of Paufler & Weber (1999) to calculate the compressibility and thermal expansivity tensors directly from these linearized EoSs at any P and T. This allows the elastic behaviour of any direction in the to be described in a consistent manner and the principal axes of the tensors, which include the directions of greatest and least strain, to be unambiguously defined. We also document how this method is implemented in a new version of the eos module in the CrysFML software library and within the established EosFit7c program for EoS calculations (Angel, Alvaro & GonzalezPlatas, 2014). A new utility in the EosFit7c program allows the user to perform all of the data analysis and EoS calculations described in this paper. The program is freely available in compiled form for Windows, Linux and macOS operating systems at https://www.rossangel.net, together with example data sets and complete documentation. The CrysFML subroutine library (RodriguezCarvajal & GonzalezPlatas, 2003) is open source and is available at https://code.ill.fr/scientificsoftware/crysfml. The architecture of the EosFit7c program also allows it to be called from other software such as MATLAB (The MathWorks Inc., Natick, MA, USA) to perform EoS calculations without the need to crosscompile software directly with the CrysFML library.
2. Constrained equations of state
2.1. The theoretical basis
The fundamental constraint on the unitcell parameters of a crystal is that they must obey the relationship
which for ease of notation we will write as
The A in (2) represents the entire squareroot expression in (1), which involves only the unitcell angles. The derivative of (2) relates changes in the volume to the changes induced in the unitcell parameters:
As an example, it follows that the volume compressibility is also related to the compressibilities of the individual cell parameters [e.g. ] as
An equivalent expression exists for thermal expansivity, . In crystal systems of orthorhombic and higher symmetry, the unitcell angles are fixed by symmetry and the volume compressibility is just the sum of the linear compressibilities of the unitcell axes:
Further differentiation of (4) gives the relationship between the first pressure derivatives of the compressibilities of crystals:
The volume EoS is usually parameterized and expressed in terms of the isothermal Reuss bulk modulus, K = . The linear compressibilities are the inverse of the corresponding isothermal Reuss linear moduli, e.g. , so from (4) the relationship between these moduli must always be
The relationship between the pressure derivative of any individual compressibility and the pressure derivative of its corresponding modulus is
so that the relationship (6) between the first pressure derivatives of the compressibilities can be expressed in terms of moduli and their pressure derivatives as
Thus, in order to be completely consistent in the description of the properties of a but also those expressed in equations (7) and (9), at all pressures.
the linearized EoSs of the unitcell axes and the volume EoS must obey not only the relationship given by (1)Examination of these equations shows that this is only possible for three special cases: first, when the compressibilities, moduli and their derivatives are identical for all directions in the crystal, which is only true in cubic crystals; second, when the pressure derivatives of the moduli are all zero, so that all of the moduli remain equal to their roompressure values at all pressures; and third, when the pressure derivatives () of the compressibilities are independent of pressure, which allows the constraints expressed in equations (4) and (6) to be met at all pressures. However, the last two cases are physically unrealistic, as nonzero pressure derivatives describe the stiffening of the structure under compression, and the rate of this stiffening changes with pressure, which is an intrinsic property of all commonly used EoSs.
As a simple example, the Murnaghan EoS assumes a linear variation of K (or M for the cell axes) with pressure, corresponding to constant values of K′ and M′. However, equation (8) shows that constant K′ and M′ do not correspond to constant values of the pressure derivatives of the compressibilities . Therefore, the Murnaghan EoS can only meet the criterion set by (4) at a single pressure. Even if all of the axes have the same value of M′, equation (9) shows that K′ for the volume calculated from the unitcellparameter EoS with constant M′ will vary with pressure, because the values of M increase with increasing pressure. But the corresponding Murnaghan EoS for the volume will have, by definition, constant K′. The consequences of this are shown in Fig. 1, with a simple example of the volume variation of a soft tetragonal crystal predicted from the axial EoSs for the unitcell parameters and the EoS for the volume. The parameters V_{0}, K_{0}^{} and of the volume EoS at the reference conditions (Table 1) are exactly those required by the constraints given in equations (1)–(9) above. However, although all of the elastic properties of the cell and volume obey these constraints at the reference conditions, they do not obey them at any other pressure. For example, the axial moduli at 1 GPa will be , so M_{a}^{} = 45 GPa and M_{c}^{} = 40 GPa, from which one obtains [equation (7)] a bulk modulus of K = 14.4 GPa, whereas the bulk modulus predicted by the Murnaghan EoS for the volume (Table 1) would be 15.625 GPa. The Birch–Murnaghan thirdorder EoS produces similar discrepancies in the bulk modulus values. Using the same roompressure parameters (Table 1), it yields for a pressure of 1 GPa values of M_{a}^{} = 44.20 GPa and M_{c}^{} = 38.09 GPa, which imply a bulk modulus of K = 13.99 GPa, whereas that from the volume EoS is 15.32 GPa. Therefore, the volume and bulk modulus predicted from the volume EoS show increasing divergence from the values predicted by the axial EoS as pressure is increased (Fig. 1); the same behaviour is exhibited by all other EoS formulations.

The rate at which these differences in volumes and bulk moduli increase with pressure depends on the bulk modulus, the elastic anisotropy of the crystal [equations (4) and (7)] and in particular on the anisotropy of the pressure derivatives of the axial compressibilities [equations (6) and (9)]. Therefore, the differences in volume calculations depend on the amount of compression (fractional volume decrease) applied to the crystal, not on the pressure itself. The difference in predicted volumes for the example in Fig. 1 reaches 1% at approximately 15% volume compression, which occurs at a pressure of about 2 GPa for this example, corresponding to ∼K_{0}/5. Thus, a stiffer material with the same degree of anisotropy will only show similar differences in volumes at higher pressures. For example, the volume differences for TiS_{2} are only slightly larger than the estimated experimental uncertainties up to pressures of 9 GPa (Fig. 2), because it has a K_{0}^{} three times larger than the example in Fig. 1. Therefore, whether or not the effects of anisotropy on the calculation of cell parameters and volumes is significant depends not only on the elastic properties of the crystal but also on the precision of the experimental data, the pressure and compression range being considered, and the precision and internal consistency required in calculations.
Clearly, the unitcell volume, the unitcell parameters and their derivatives [equations (1)–(9)] are not independent quantities, and therefore the description of their behaviour with P and/or T by independent EoSs introduces false additional We have also demonstrated here that, if the variations of all of the unitcell parameters and volume are described by realistic EoSs, then the constraints on the relationship between the elastic properties of the volume and the unitcell parameters cannot be met exactly at all pressures except in special cases. Therefore, it is not physically consistent to describe the anisotropic evolution of a with P and T with completely independent EoSs. The following examples, grouped by show how this problem can be overcome and a fully selfconsistent description of the variation of the parameters of a of a crystal can be obtained. Most of the examples illustrate the problem in the context of isothermal compression, for which the effects are largest. However, exactly the same methodology can be applied to or, in general, an EoS describing how the unitcell parameters change with P and T, as illustrated below with the example of quartz.
2.2. Cubic crystal system
For cubic crystals the general relationship (1) between the unitcell edge a and the volume is simplified to V = a^{3}, and it follows from the constraint equations given above that the elastic properties of all of the unitcell edges have the same relationship to those of the volume. Thus,
This means that the variation of the a or an EoS for the volume with the parameter relations given in (10), which will then hold for all temperatures and pressures, unlike the example shown in Fig. 1. In addition, because the and compressibility are secondrank tensor properties, these are identical for all directions in a cubic crystal.
of a cubic crystal can be described equally well with either a linearized EoS for the unitcell parameter2.3. Uniaxial crystal systems
In these crystal systems the values of the unitcell angles are fixed by the symmetry so the angle factor A [equations (1) and (2)] is 1.0 in the tetragonal system and A = sin(120°) for the trigonal and hexagonal systems in their conventional settings. Being constant, the angle factor A does not contribute to the volume derivatives as defined by equation (3). Therefore, for the uniaxial crystal systems the variation in the unitcell parameters can be described by fitting EoSs to any two of a, c and V, setting the EoS of b to be equal to that of the a axis, and calculating the properties of the third symmetrically independent parameter from the other two. Fig. 2 shows the unitcell variation of TiS_{2} with pressure (Allan et al., 1998) treated in this way; this is used as an example because TiS_{2} has a layered structure with very strong anisotropy in both the axial moduli and their pressure derivatives, and it therefore provides a challenging test. The unitcell volume calculated from the axial Birch–Murnaghan EoSs for a and c is indistinguishable on the scale of the figure from the volume EoS fitted directly to the volume data. In detail, the average misfit to the data quantified as (P_{obs} − P_{calc}) is 0.024 GPa from the unitcell EoS and 0.015 GPa from the direct volume EoS. From room pressure to 5 GPa, in the middle of the data set where the elastic parameters are best constrained, the difference in K and K′ between the two descriptions of the volume is less than 1.5 e.s.d., but then increases to the level of 2 e.s.d. at 8 GPa. Thus, for this very anisotropic example, the fit of the axial EoS represents the volume properties almost as well as the direct fit to the volume data.
The mineral zircon, ZrSiO_{4}, has tetragonal symmetry and is also strongly elastically anisotropic with the c axis being almost twice as stiff as the a and b axes (e.g. Özkan et al., 1974; Ehlers et al., 2022). This makes it challenging to determine the linear modulus of the c axis from just measurements of the unitcell parameters under pressure. Fig. 3 shows that the behaviour of the c axis of zircon can be well represented by calculating it as c = V/a^{2}, with the values of V and a obtained from their corresponding EoSs fitted to the experimental data (Ehlers et al., 2022). The calculated caxis variation gives a value of M_{c0} that is within 1 e.s.d. of the value determined by a combined fit to elasticity and compressional data (Ehlers et al., 2022), and also reproduces the slight softening ( < 0) from 0 to ∼4 GPa that is apparent in the data (Fig. 3).
A further stringent test of the method is provided by the extreme values of properties in the neighbourhood of continuous displacive phase transitions (e.g. Carpenter, Salje & GraemeBarber, 1998; Carpenter & Salje, 1998). In quartz, the α–β at ∼848 K and room pressure (e.g. Carpenter, Salje, GraemeBarber, Wruck et al., 1998) is accompanied by almost infinite and compressibilities in the lowtemperature α phase, the latter corresponding to the softening of the isothermal bulk modulus and linear moduli to zero. Fig. 4 shows that the behaviour of the c axis calculated as c = V/a^{2}sin(120°) from the volume and aaxis EoSs is reproduced just as well as the aaxis properties are predicted by the EoS fitted directly to the aaxis data. In particular, this calculation for the c axis captures the reduction of the linear modulus M_{c} as the transition is approached in the α phase, the rapid recovery of this modulus in the β phase, and the negative of the c axis in the β phase.
2.4. Orthorhombic crystal system
For orthorhombic crystal systems the unitcell angles are constrained to be 90°, the factor A = 1.0, and the volume of the unitcell is simply the product V = abc. Therefore, the anisotropy of the elastic properties of the of an orthorhombic crystal can be described completely by specifying the EoSs of three of these quantities and calculating the fourth from them. As for crystal systems of higher symmetry, this also allows the compressibility of any direction in the crystal to be calculated at any P and T from the EoS. As an example, Fig. 5 shows the pressure variation of the length of the [111] lattice vector and the d spacing of the (221) planes in pure forsterite, Mg_{2}SiO_{4}, calculated using the axial EoSs for the unitcell axes.
2.5. Monoclinic crystal system
We use the conventional setting of the monoclinic b axis being unique to describe our methodology for monoclinic crystals. Equivalent expressions for the relationships between cell parameters follow if a different cell setting with a different unique axis is chosen. For b unique, the angle factor in equations (1) and (2) is , and the unitcell parameters and volume are related by . It follows that the can be completely described by four of these five parameters. Selfconsistent calculations can therefore be achieved in different ways. We find that for precise data with β angles greater than 100° both the angle and the lengths of all directions can be well reproduced (Fig. 6) by using EoSs fitted to the volume and the cell parameters, and calculating . Note that because of the trigonometric identity one has to specify whether the β angle is obtuse or acute (greater or less than 90°).
with theThe uncertainty in the β angle (in degrees) calculated in this way is , where is the combined uncertainty in V/abc. The factor of in the denominator ensures that the uncertainty in β angle approaches infinity as β approaches 90°; the value calculated from the EoS becomes less reliable and more sensitive to uncertainties in V/abc to the extent that the value of β can be completely dominated by fitting errors and noise in the EoS and the underlying data, together with numerical precision in the computer code. The predicted behaviour of both the β angle and consequently some directions in the crystal can become completely wrong and unphysical as shown in Fig. 7. An additional problem is that uncertainties in the EoS parameters can sometimes result in the ratio V/abc calculated from their EoS being greater than 1, for which no value of β can be defined.
This is especially a problem in soft materials such as molecular crystals, in which the rapid increase in moduli with pressure and strong anisotropy in M′ means that the volume variation with pressure includes a significant component from the variation in the β angle [equation (3)], which in turn means that the EoSs do not fit the data as well as for stiffer materials over the same pressure range (compare Figs. 6 and 7). Because , one can also describe the by an EoS for this d spacing, or an EoS for , plus EoSs for three of a, b, c and V. However, when the β angle approaches 90°, the evolution of the d spacings with P and T approaches that of the corresponding unitcell axes a and c. The value of the β angle obtained from the two EoSs as then suffers the same unreliability just described.
Because the problem is in the recovery of the β angle from the calculated cell parameters, an alternative approach is to describe the variation of the angle directly in terms of P and T. There is no constraint from theories of elasticity or EoSs that define how unitcell angles should change with T or P, so a phenomenological approach can be adopted. For example, a polynomial of the angle in P and T or a polynomial of a trigonometric function or functions of the angle can be used. The only requirement is that the function chosen provides an invertible relationship that defines a unique angle for a given P (or T) and a unique P for a given angle and T. Then, given that the β angle is defined, EoSs for three of the four parameters a, b, c and V are required to complete the description of the anisotropic properties of the While this loses the advantages of using only the EoS, especially with respect to extrapolation to conditions beyond the range of the data, Figs. 6 and 7 show that, if the variation of the monoclinic angle is accurately represented by a polynomial, the method provides a good description of the variation of the lengths of all directions within the crystal while ensuring complete internal consistency between all calculated cell dimensions at a given P and T.
2.6. Triclinic crystal system
In the triclinic system there are six cell parameters and the volume, of which six are required to define the geometry of the lattice and a, b and c and also for the three d spacings 100, 010 and 001, then the is completely defined. However, this is not the case, because the relationships between the unitcell angles and the d spacings take the form
Thus it appears that if EoSs are defined forThis defines the unitcell angle α in terms of V, d_{100} and two unitcell parameters, but the unitcell angles are required to calculate the volume through equation (1). If one wants to retain the advantages of using only EoSs to describe the it is therefore necessary specify all seven EoSs for the three cell axes, three principal d spacings and volume. This introduces a spurious degree of freedom because the cell parameters and volume are treated independently and therefore the selfconsistency between the EoSs is not ensured. This approach also suffers from the same problem of uncertainties in the unitcell angles if they approach 90° (Fig. 8), which then leads to the prediction of unphysical behaviour of certain directions within the crystal, such as the [111] lattice vector shown in Fig. 8.
The alternative approach is to use a separate polynomial for each of the unitcell angles, in which case the unitcell is uniquely and consistently determined in combination with EoSs for three of the four parameters a, b, c and V. Fig. 8 shows that this provides a more accurate description of the variation of lengths of different lattice vectors and d spacings within the of triclinic crystals. The volume calculated from the six cell parameters is indistinguishable from the EoS fitted to the volume, at the scale of the measurement uncertainties in the data.
2.7. Implementation in EosFit
These procedures have been implemented in the EoS module of the CrysFML Fortran subroutine library (RodriguezCarvajal & GonzalezPlatas, 2003), including the description of unitcell angles by polynomials. Calculations can be performed with the CELL utility of the EosFit7c program. EoSs for volumes and cell parameters and, for triclinic crystals, the EoSs of the principal d spacings can be loaded from EosFit .eos files or be input directly by the user. Polynomials to describe the variation of the unitcell angles in monoclinic and triclinic crystals can also be stored in, and loaded from, .eos files. The current version supports simple polynomials up to third order in P and T, thus allowing a maximum of ten coefficients of all orders up to P^{3} and T^{3}, together with three crossterms PT, PT^{2} and P^{2}T. The program architecture will allow extension to higherorder polynomials or the implementation of other functions to describe angle variation, if these are required.
Apart from the constraints imposed by the symmetry of the a = b in uniaxial systems when the unique axis is specified as c, the choice of which unitcell parameter to calculate from the EoSs of the others is a matter of user choice. Therefore, the program only checks and restricts the choice of EoSs to those consistent with the Once a set of EoSs have been loaded, the program checks also for internal consistency of the EoSs, for example that the pressure scales and the units for cell parameters and volume are consistent between all of the EoSs. Calculations are prevented until the EoSs necessary for the are loaded. The available commands include the calculation and output of the full unitcell parameters at a P and T chosen by the user, and the output to a file listing the cell parameters over a range of P or T. The full properties (length, modulus, thermal expansivity etc.) of any individual direction can be calculated, as well as lines in P–T space of constant length (i.e. the linear equivalents of isochors of volume) and lines of constant ratio between a pair of directions. Pressure can be calculated from the length of any direction or d spacing in the crystal. The results of all of these calculations can be written into text files so that the calculated values can be imported into plotting software.
for example thatWith the approach described here there are four different types of directions (or volume) whose properties are calculated in different ways. The properties of unitcell axes (or the volume) for which EoSs are loaded, or which are symmetrically equivalent to a cell axis for which an EoS has been loaded, are calculated directly from their EoSs. All of the properties of other unitcell axes or the volume (i.e. those without a loaded EoS) are calculated directly from the loaded EoSs from the other axes through the relationships given in equations (1)–(9). For other, general, directions in the the cell parameters are first calculated from the loaded EoSs and the is constructed, from which the length of the chosen vector is calculated. Properties such as the linear modulus or thermal expansivity are calculated, numerically from a spline to a series of length calculations over a small range of P or T, as appropriate. This numerical approach, equivalent to that used by Paufler & Weber (1999) and Langreiter & Kahlenberg (2015), means that the calculated properties of general directions have lower precision than those calculated for unitcell axes for which EoSs are available, and the precision decreases as the order of the volume derivative increases; calculated values of the linear modulus M are more reliable than the derivatives M′ or .
3. Tensor properties
The directions in a ). In crystals of orthorhombic or higher symmetries the eigenvectors are constrained to be fixed in direction to be parallel to the axes of the conventional and only their magnitude is a property of the The magnitudes of the principal axes, and the corresponding tensor components, are therefore given by the compressibility and thermal expansivity of the corresponding unitcell axes, which are in turn defined by their axial EoSs.
of maximum and minimum or compression are given by the eigenvectors of the corresponding secondrank tensors. The eigenvectors are also called the principal axes of the tensor (Nye, 1957In monoclinic crystals one principal axis is constrained to lie parallel to the diad axis. The other two principal axes lie in the plane perpendicular to the diad axis, but their directions in this plane are not constrained by symmetry and they can therefore rotate around the unique axis with a change in P or T. In triclinic crystals there are no constraints on the directions of the principal axes of these property tensors, meaning that they are completely free to rotate. In these cases, the full tensors must be calculated. The traditional method was to first calculate the strains between a pair of determinations of the of the crystal at different P or T (e.g. Ohashi & Burnham, 1973; Schlenker et al., 1978) and then divide the strains by the P or T increment between the measurements to obtain the compressibility or tensors. This approach has several widely recognized disadvantages [discussed by Jessen & Küppers (1991), Paufler & Weber (1999) and Knight (2010)]. First, it returns a property tensor that is an average over the finite P or T range between the two data points. Second, the exact values of the tensor components are dependent on the initial and final conditions; thermal expansivity tensors calculated across different temperature intervals with the same midpoint temperature will have different values of their components. Third, in cases where there are significant changes in the unitcell angles between the two unit cells, the directions of the eigenvectors of the strain and property tensors are defined relative to the reference cell under the assumption that the Cartesian axial system used to describe the tensors is fixed and does not rotate during the P and T change. As a consequence, the directions of the eigenvectors relative to the cell at the midpoint conditions are undefined, and a different choice of reference cell will lead to principal axes with slightly different directions relative to the underlying crystal structure.
To avoid these problems, Paufler & Weber (1999) developed the alternative approach of defining the tensor directly in terms of the derivatives of the unitcell parameters with respect to temperature. They used polynomials to determine the unitcellparameter changes with temperature, from which the temperature derivatives are derived analytically and the tensor is unambiguously defined for any temperature within the range of the data. The values of the tensor components and the directions of the eigenvectors are not averages but are defined at the same temperature as one another (Paufler & Weber, 1999; Langreiter & Kahlenberg, 2015). The same approach can be applied to determine the compressibility tensor from a series of unitcell determinations with pressure (Knight, 2010). The method does not rely on the algebraic form of the equations used to describe the cellparameter variations, so the and compressibility tensors can be calculated from the temperature and pressure derivatives of the linearized EoSs of crystals.
3.1. Implementation in EosFit
In the EoS module of the CrysFML library (RodriguezCarvajal & GonzalezPlatas, 2003) we implement the principles of the method of Paufler & Weber (1999). In crystal systems of higher than monoclinic symmetry, the only nonzero tensor components for the conventional choices of orientation of the Cartesian axes are the diagonal components. These are calculated directly from the set of selfconsistent linearized EoSs fitted to the unitcell parameters, with the components of the compressibility tensor being the inverse of the corresponding linear moduli. In monoclinic and triclinic crystal systems, the values of the tensor components also depend on the choice of orientation of the Cartesian axes relative to the crystallographic axes, and on the derivatives of real and reciprocal angles. If the realspace angles are described by polynomial functions, their derivatives are taken directly from the functions. Otherwise their derivatives are calculated by splines over a small range of T and P around the point at which the tensor is required. The derivatives of the reciprocalspace angles are always calculated by splines. The equations of Paufler & Weber (1999) provide the equations for one orientation of the Cartesian axes, those of Tribaudino et al. (2011) for a second, and those for two further commonly used orientations of the Cartesian axes were derived from these and coded into the CrysFML library. All three possible orientations of monoclinic unit cells are also supported.
The CELL utility of EosFit7c allows the user to calculate the property tensors from a set of selfconsistent linearized EoSs. Facilities are provided to list the tensor components over a range of P or T, along with the eigenvalues of the properties and the directions of the principal axes (eigenvectors) relative to the Cartesian axes and the real and reciprocal unitcell axes. A calculation of the nearest lowindex plane normal and lattice vector is also provided to aid the user to relate the directions of the principal axes to the structure of the crystal. If the original unitcell data are available, the results can be compared with the property tensors calculated by the finite difference between consecutive cell determinations by the STRAIN utility of EoSFit7c. Note that the STRAIN utility supports several different strain definitions (Zotov, 1990), which should become equivalent in the infinitesimal limit represented by the derivative approach used here.
4. Conclusions
We have proved algebraically in equations (1)–(9) that the volume variation of a crystal with P and T predicted from its unitcell parameters described by conventional linearized (axial) EoSs is never exactly the same as that from a volume EoS, unless the crystal is cubic. Using independent EoSs for the symmetrically distinct unitcell parameters and the volume is therefore always physically inconsistent and introduces false additional This applies to all types of EoS. The physical inconsistency arises from the anisotropy in the pressure derivatives of the axial compressibilities. Whether or not this is a significant issue for a given crystal depends on the degree of elastic anisotropy, the range of compression of interest, and the precision of experimental data or the precision and internal consistency required in calculations. Discrepancies are greatest for highly anisotropic soft materials such as many metal–organic frameworks and molecular crystals. If internal consistency is important, we have described how the full variation of the unitcell parameters, volumes and elastic properties of a crystal can be described by suitable combinations of the axial and volume EoSs for each and we have provided the tools for such calculations in the 2021 release of the EosFit7c program. The consistency of such calculations has been demonstrated with examples from all crystal systems. The example of zircon demonstrates that the method is also useful for determining the elastic properties of very stiff directions in a crystal that are too stiff to determine precisely by direct experimental measurement.
Acknowledgements
We thank Tiziana BoffaBallaran (Bayreuth), Volker Kahlenberg (Innsbruck), Ronald Miletich (Wien) and Fabrizio Nestola (Padova) for extensive discussions over the years about how best to describe unitcellparameter variations, and Mara Murri (Milano), Daria Pasqual (Padova) and an anonymous reviewer for careful reviews of the manuscript.
Funding information
This project was funded from the European Research Council under the European Union's Horizon 2020 research and innovation programme grant agreement 714936 TRUE DEPTHS to MA. This work has been partially supported by Agencia Estatal de Investigación under the national programme of sciences and technological materials (PID2019106383GBC44/AEI/10.13059/501100011033) to JGP. MM was supported by a research fellowship from the Alexander von Humboldt foundation. MA acknowledges support from Italian MIUR Progetto di Ricerca di Interesse Nazionale (PRIN) No. 2017ZE49E7 to M. Scambelluri.
References
Allan, D. R., Kelsey, A. A., Clark, S. J., Angel, R. J. & Ackland, G. J. (1998). Phys. Rev. B, 57, 5106–5110. Web of Science CrossRef CAS Google Scholar
Alvaro, M., Mazzucchelli, M. L., Angel, R. J., Murri, M., Campomenosi, N., Scambelluri, M., Nestola, F., Korsakov, A. V., Tomilenko, A. A., Marone, F. & Morana, M. (2020). Geology, 48, 24–28. Web of Science CrossRef Google Scholar
Anderson, O. L. (1995). Equations of State of Solids for Geophysics and Ceramic Science. Oxford University Press. Google Scholar
Angel, R. J. (2000). Rev. Mineral. Geochem. 41, 35–59. Web of Science CrossRef Google Scholar
Angel, R. J., Alvaro, M. & GonzalezPlatas, J. (2014). Z. Kristallogr. 229, 405–419. Web of Science CrossRef CAS Google Scholar
Angel, R. J., Alvaro, M., Miletich, R. & Nestola, F. (2017). Contrib. Mineral. Petrol. 172, 29. Web of Science CrossRef Google Scholar
Angel, R. J., Mazzucchelli, M. L., Alvaro, M., Nimis, P. & Nestola, F. (2014). Am. Mineral. 99, 2146–2149. Web of Science CrossRef Google Scholar
Angel, R. J., Miozzi, F. & Alvaro, M. (2019). Minerals, 9, 562. Web of Science CrossRef Google Scholar
Angel, R. J., Mosenfelder, J. L. & Shaw, C. S. J. (2001). Phys. Earth Planet. Inter. 124, 71–79. Web of Science CrossRef CAS Google Scholar
Carpenter, M. A. & Salje, E. K. H. (1998). Eur. J. Mineral. 10, 693–812. CrossRef CAS Google Scholar
Carpenter, M. A., Salje, E. K. H. & GraemeBarber, A. (1998). Eur. J. Mineral. 10, 621–691. CrossRef CAS Google Scholar
Carpenter, M. A., Salje, E. K. H., GraemeBarber, A., Wruck, B., Dove, M. T. & Knight, K. S. (1998). Am. Mineral. 83, 2–22. CrossRef Google Scholar
ConesaEgea, J., GallardoMartínez, J., Delgado, S., Martínez, J. I., GonzalezPlatas, J., FernándezMoreira, V., RodríguezMendoza, U. R., Ocón, P., Zamora, F. & AmoOchoa, P. (2017). Small, 13, 1700965. Google Scholar
Ehlers, A. M., Zaffiro, G., Angel, R. J., BoffaBallaran, T., Carpenter, M. A., Alvaro, M. & Ross, N. L. (2022). Am. Mineral. https://doi.org/10.2138/am20217731. Google Scholar
Gonzalez, J. P., Mazzucchelli, M. L., Angel, R. J. & Alvaro, M. (2021). J. Geophys. Res. Solid Earth, 126, e2021j, b022080. Google Scholar
Holzapfel, W. B. (2001). Z. Kristallogr. 216, 473–488. Web of Science CrossRef CAS Google Scholar
Jessen, S. M. & Küppers, H. (1991). J. Appl. Cryst. 24, 239–242. CrossRef Web of Science IUCr Journals Google Scholar
Knight, K. (2010). Phys. Chem. Miner. 37, 529–533. Web of Science CrossRef CAS Google Scholar
Kroll, H., Kirfel, A. & Heinemann, R. (2014). Eur. J. Mineral. 26, 607–621. Web of Science CrossRef CAS Google Scholar
Lakshtanov, D. L., Sinogeikin, S. V. & Bass, J. D. (2007). Phys. Chem. Miner. 34, 11–22. Web of Science CrossRef CAS Google Scholar
Langreiter, T. & Kahlenberg, V. (2015). Crystals, 5, 143–153. Web of Science CrossRef Google Scholar
Murshed, M. M., Mendive, C. B., Curti, M., Šehović, M., Friedrich, A., Fischer, M. & Gesing, T. M. (2015). J. Solid State Chem. 229, 87–96. Google Scholar
Mazzucchelli, M. L., Reali, A., Morganti, S., Angel, R. J. & Alvaro, M. (2019). Lithos, 350–351, 105218. Web of Science CrossRef Google Scholar
Nye, J. F. (1957). Physical Properties of Crystals. Oxford University Press. Google Scholar
Ohashi, Y. & Burnham, C. W. (1973). Am. Mineral. 58, 843–849. CAS Google Scholar
Özkan, H., Cartz, L. & Jamieson, J. C. (1974). J. Appl. Phys. 45, 556–562. Google Scholar
Paufler, P. & Weber, T. (1999). Eur. J. Mineral. 11, 721–730. CrossRef CAS Google Scholar
Poe, B., Romano, C., Nestola, F. & Smyth, J. (2010). Phys. Earth Planet. Inter. 181, 103–111. Web of Science CrossRef CAS Google Scholar
RodriguezCarvajal, J. & GonzalezPlatas, J. (2003). IUCr Comput. Commission Newsl. 1, 50–58. Google Scholar
Schlenker, J. L., Gibbs, G. V. & Boisen, M. B. (1978). Acta Cryst. A34, 52–54. CrossRef CAS IUCr Journals Web of Science Google Scholar
Stixrude, L. & LithgowBertelloni, C. (2005). Geophys. J. Int. 162, 610–632. Web of Science CrossRef Google Scholar
Tribaudino, M., Bruno, M., Nestola, F., Pasqual, D. & Angel, R. J. (2011). Am. Mineral. 96, 992–1002. Web of Science CrossRef CAS Google Scholar
Zotov, N. (1990). Acta Cryst. A46, 627–628. CrossRef Web of Science IUCr Journals Google Scholar
This is an openaccess article distributed under the terms of the Creative Commons Attribution (CCBY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.