## computer programs

*QUANTAS*: a Python software for the analysis of thermodynamics and elastic behavior of solids from *ab initio* quantum mechanical simulations and experimental data

^{a}Department of Biological, Geological and Environmental Sciences, University of Bologna, Piazza di Porta San Donato 1, Bologna, 40126, Italy^{*}Correspondence e-mail: gianfranco.ulian2@unibo.it, giovanni.valdre@unibo.it

Mineralogy, petrology and materials science are fundamental disciplines not only for the basic knowledge and classification of solid phases but also for their technological applications, which are becoming increasingly demanding and challenging. Characterization and design of materials are of utmost importance and usually need knowledge of the thermodynamics and mechanical stability of solids. Alongside well known experimental approaches, in recent years the advances in both quantum mechanical methods and computational power have placed theoretical investigations as a complementary useful and powerful tool in this kind of study. In order to aid both theoreticians and experimentalists, an open-source Python-based software, *QUANTAS*, has been developed. *QUANTAS* provides a fast, flexible, easy-to-use and extensible platform for calculating the thermodynamics and elastic behavior of crystalline solid phases, starting from both experimental and *ab initio* data.

Keywords: *QUANTAS* computer program; solid phases; thermodynamics; equations of state; second-order elastic constants; quantum mineralogy.

### 1. Introduction

An important issue for materials scientists (mineralogists, petrologists, solid-state chemists and physicists, and materials engineers) is to obtain good-quality data of the physical and chemical properties of a solid phase at different pressure and temperature conditions. This fundamental knowledge can drive further research on both the thermodynamic stability of the system under investigation, for example the construction of phase diagrams, and possible innovative specific technological applications.

Among several experimental methods, two approaches are commonly employed to study the elastic properties of a solid material: (i) hydrostatic compression and (ii) uniaxial/biaxial deformation. In the first case, the mechanical data on the solid can be obtained from the variation of unit-cell volume/axis lengths with pressure (*e.g.* by using diamond-anvil cells) by means of equation-of-state (EoS) fitting. Such experiments are conducted either at constant temperature or by varying it during diffraction experiments (Fukui *et al.*, 2003; Gatta *et al.*, 2013, 2015; Pawley *et al.*, 1995). The greatest challenge for this kind of investigation is to describe the structural dynamics on simultaneous variation of both pressure and temperature, which still requires very complex experimental setups. In the second case, the stress–strain relationship allows the determination of the fourth-rank elastic tensor of the material, whose components are also called `elastic moduli' and `second-order elastic constants' by physicists and engineers, respectively. The analysis of the elastic properties can reveal many details about both single-crystal and polycrystalline behaviors and is of utmost importance for designing materials with tailored mechanical properties. Several methodologies could be employed to obtain the elastic moduli, such as Brillouin spectroscopy (Jiang *et al.*, 2006) and shock wave experiments (Duffy *et al.*, 1991; Yoon & Katz, 1976).

In the past two decades, quantum mechanical (QM) methods have been increasingly applied in solid-state and materials science fields for both simple and complex systems. For instance, in mineralogy, this new scientific branch even started being referred to as `quantum mineralogy'. The use of quantum mechanical approaches is mainly for structural, electronic, vibrational and catalytical (surface) properties. The simulation of the thermomechanical behavior of solids is only recent, principally because of the high computational costs required by such calculations. The basis to obtain thermomechanical insights of crystals is given by the quasi-harmonic approximation (QHA), as described by Anderson (1995). The QHA offers a simple formulation valid for several solids presenting different bond types, such as pure covalent (Ottonello, Zuccolini & Civalleri, 2009), mixed covalent–ionic (Ottonello *et al.*, 2010; Ulian & Valdrè, 2015*a*), pure ionic (Erba, 2014) and mixed covalent–dispersive (Ulian & Valdrè, 2015*b*,*c*). This method, which is derived from the harmonic treatment of the of solids, introduces an explicit dependence of phonon frequencies on volume, thus overcoming the well known limitations of the harmonic approximation (HA) as reported by Baroni *et al.* (2010).

Several formulations describing the volume dependence of vibrational modes have been proposed: (i) Grüneisen's mode-γ parameters (Anderson, 1995; Ottonello, Civalleri *et al.*, 2009); (ii) a volumetric polynomial fit of the thermodynamic quantities calculated by the HA (Prencipe *et al.*, 2011); and (iii) a direct polynomial fit of the individual phonon frequencies with respect to volume (Erba, 2014; Erba *et al.*, 2015). In all of these approaches, it is required to calculate the phonon modes of the solid under investigation at different unit-cell volumes. The main difference between the thermodynamic and phonon interpolation schemes resides in the number of fitting procedures to be performed, which is usually lower for the direct polynomial fit of the individual phonon frequencies. Typically, QHA calculations are carried out considering only hydrostatic pressure variations on the However, it is known that phonon frequencies depend also on the changes in the cell parameters expressed as strains, as also reported and discussed in recent literature (Murri *et al.*, 2018; Ulian & Valdrè, 2018; Destefanis *et al.*, 2019). However, in the present work we refer to the most common approach to QHA.

Nowadays, there are various quantum chemistry codes that can calculate the quantities necessary to perform QHA and elastic analyses of solids. Phonon modes and elastic tensors can be obtained directly from the outputs of the quantum mechanical software thanks to automated algorithms implemented in the package and/or indirectly. For example, the elastic moduli require knowledge of the stress–strain relationship based on total energy calculations, which are performed for a systematic series of deformations of the ). In the simplest approach, it is necessary to (i) determine the number and type of deformations, (ii) calculate the total energy of the crystalline solid for each deformed cell, and (iii) numerically derive the elastic tensor components from the energy versus deformation curves. Both the *CRYSTAL* (Dovesi *et al.*, 2018) and the *VASP* (Kresse & Furthmüller, 1996; Kresse & Hafner, 1993) codes can provide the elastic moduli directly, because they implement routines that automatically perform the steps cited above. If the quantum mechanical code does not provide these routines, it is necessary to perform these operations by hand or rely on external codes/scripts that both generate specific input files for the QM software and analyze the output results, as performed for example by the *ElaStic* scripts (Golesorkhtabar *et al.*, 2013). The same applies for the calculation of the phonon modes required by QHA analysis and of the equation of state of the solid under analysis.

In addition, there is an increasing need of codes that can post-process data from both theoretical simulations and experimental measurements, to obtain other information such as the thermodynamic and thermoelastic properties, or the directional and averaged mechanical behavior. The available tools in the scientific community are usually very specific to just one or a few of the cited analyses, and a comprehensive and inclusive platform is still missing.

Starting from all these considerations, and trying to satisfy the needs of both experimentalists and theoreticians, we developed *QUANTAS* (acronym of `quantitative analysis of solids'), a software platform that can aid solid-state and materials scientists of different fields in obtaining a multiplicity of properties of a crystalline phase.

At present, *QUANTAS* allows the user to calculate

(i) the thermodynamic and thermoelastic properties of a material at selected pressure and temperature conditions from *ab initio* quantum mechanical results;

(ii) the equation of state of crystalline phases from both theoretical and experimental data;

(iii) elastic properties derived from the second-order elastic moduli, independently of the means used to obtain them.

The present work is organized as follows: Section 2 lists and briefly explains some of the extant codes to better compare their functions with those of *QUANTAS*. Section 3 reports and discusses the non-functional requirements of *QUANTAS*. Sections 4–6 show all the functionalities and implementations, and Sections S1–S4 in the supporting information present, for the sake of completeness, some test cases used to validate *QUANTAS*. Finally, possible future developments of the presented package are explained and discussed. This paper is not intended as a user manual but provides just some of the details that can be found in the online documentation (see Section 3). The focus is on the basic science, general concepts and engineering aspects of *QUANTAS*.

### 2. Extant software

In the following, we provide a brief presentation of other related codes, explaining the differences with the *QUANTAS* software.

#### 2.1. *Phonopy*

*Phonopy* (Togo & Tanaka, 2015) is an open-source code written in Python (with some functions in C) which aims to calculate phonon properties of crystalline solids. It has several pre-processing routines, *e.g.* creation of input files (unit-cell structures) for the calculation of unit-cell energy related to atomic displacements, and post-processing ones, such as determination and plotting of phonon bands in *k* paths, specifically developed for first-principles simulations. *Phonopy* has several interfaces for different *ab initio* codes and is widely employed by quantum chemical researchers. It allows also the calculation of harmonic thermodynamic properties and includes a quasi-harmonic approximation framework. It employs a thermodynamic interpolation scheme, using either the Murnaghan (1937) or the Vinet (Vinet *et al.*, 1987) equation of state to minimize the volume at 0 GPa. In *QUANTAS*, we provide four equation-of-state formulations and polynomials to minimize the volume at different temperature and pressure states, and both thermodynamic and phonon interpolation schemes (see below).

#### 2.2. *EosFit7*

The *EosFit* software (Angel *et al.*, 2014) is a copyleft suite of codes developed to calculate *P*–*V*–*T* EoSs and cell parameter variations with pressure *P* and temperature *T* from volume *V*, cell parameters and elasticity data (Milani *et al.*, 2017). The current version is based on the *CrysFML* Fortran library (Rodríguez-Carvajal & González-Platas, 2005, 2008) and includes both a console and a graphical user interface (Gonzalez-Platas *et al.*, 2016). It does not implement volume-integrated equation-of-state formulations that could be used by theoreticians to fit their total energy versus unit-cell volume curves.

#### 2.3. *ElAM*

*ElAM* (Marmier *et al.*, 2010), implemented in Fortran 90, is an open-source command-line software that provides analysis of the second-order elastic tensor using well known solid-state physics formulations as described by Nye (1957). It has both 2D and 3D plotting capabilities (the first in PostScript format, the second in virtual-reality modeling language format). However, the software seems discontinued since the end of 2009. It does not provide the calculation of directional and averaged seismic wave speeds.

#### 2.4. *ELATE*

A successor of *ElAM*, *ELATE* is an open-source online tool, entirely written in Python by Gaillac *et al.* (2016), which provides a detailed analysis of the second-order elastic tensor, together with both bi- and tri-dimensional plots. It is also integrated in the Materials Project, a large database of several properties of many solid phases. As for *ElAM*, no routine has been implemented to calculate the seismic wave speeds.

### 3. Non-functional requirements of *QUANTAS*

Many scientific software codes have very similar non-functional requirements, such as license type and programming language. For the development of *QUANTAS* we followed whenever possible the best practices for scientific code development that have been recently proposed by Wilson *et al.* (2014).

#### 3.1. Software license

*QUANTAS* is released as a free and open-source code under the New Berkeley Software Distribution (BSD) software license. This aims to provide the necessary conditions for the verification and validation of research data (Joppa *et al.*, 2013), while also allowing advanced users to modify the source code, collaborate on future implementations, and ensure long-term maintenance and sustainability of the software.

#### 3.2. Programming language

*QUANTAS* is developed in Python 3 (https://www.python.org), with the numerical sections implemented through the *NumPy* and *SciPy* scientific computing libraries (van der Walt *et al.*, 2011). Heavy computations are performed by a small portion of the software that was written in Cython (Behnel *et al.*, 2011), namely as C extensions for Python. The computation is in double precision for the thermodynamic/thermoelastic properties using the HA or QHA framework, whereas single precision is employed for the equation-of-state fitting and the analysis of the elastic tensor. The choice of a high-level language is motivated by its easier comprehensibility and maintenance, and also the need to provide cross-platform support. *QUANTAS* was developed keeping in mind the different environments where it could be employed, from desktop computers to servers. At the moment, the package is shipped as a Python library that can be installed by the user, but, in the near future, it will be available also from well known package repositories such as PyPI (https://pypi.org/).

#### 3.3. User interface

*QUANTAS* uses a command-line interface (CLI), which is invoked by using a Python entry point. In general, a calculation is started from the console shell (or command prompt on Windows) by typing

`quantas sub_command input_f``ile_name [options]`

The available sub-commands are related to the different routines implemented, whose names were chosen to be intuitive for users. For example, to run a quasi-harmonic approximation calculation, the sub-command that has to be called is `qha`.

In addition, a specific sub-command (`inpgen`) can be employed to generate input files for post-processing analyses of the (quasi-)harmonic approximation and second-order elastic constants:

`quantas inpgen [generator] input_f``ile [options]`

This routine is particularly useful given the complexity of extracting the phonon modes from quantum mechanical simulations and formatting them for a *QUANTAS* input file.

#### 3.4. Online services

The *QUANTAS* source code is hosted at https://github.com/gfulian/quantas, and detailed documentation on how to obtain (download), install and use the software is provided by the ReadTheDocs service (hosted on https://quantas.readthedocs.io). The latter is also home of some tutorials designed to guide users in creating the input files, running the different analyses and collecting the results. The web site and documentation were created using the *Sphinx* code (https://www.sphinx-doc.org/).

### 4. (Quasi-)Harmonic approximation

#### 4.1. Theory

In quantum mechanical simulations of periodic three-dimensional systems (crystals), the harmonic thermodynamics of any system can be obtained from its phonon properties, which can be calculated just in the central point of the first Γ point) or in several **k** points in the Several reports in the literature describe different approaches to calculating the Γ point (**k** = 0) vibrational modes and the phonon dispersion relations, whose detailed description is beyond the scope of the present work. The interested reader may refer to the fundamental work of Parlinski *et al.* (1997), Erba (2014) and Togo & Tanaka (2015).

Let us assume that, given a crystal **k** points have been sampled. It is known that 3*N* oscillators (phonons), with *N* the number of atoms in the cell, are associated with each considered **k** point. In turn, each phonon is associated with an energy level given by the harmonic expression , with *m* an integer number and *ν _{i}*(

**k**) the frequency of the harmonic oscillator. It is then possible to express the vibrational canonical of the system as

where *k*_{B} is the Statistical thermodynamics defines the *S*(*T*), the thermal internal energy *U*_{th}(*T*) and the isochoric *C*_{V}(*T*) of a system as

By substituting equation (1) into equations (2)–(4), it is possible to write the harmonic expression of the thermodynamic properties as reported by Prencipe *et al.* (2011) as

While the harmonic approach has been successfully adopted in predicting vibrational (spectroscopic) and thermodynamic properties of several systems (Belmonte, 2017; Erba, 2014; Prencipe *et al.*, 2004; Ulian *et al.*, 2021; Ulian & Valdrè, 2019), it suffers from several well known limitations. In fact, many important properties of crystalline materials are wrongly described within this framework; for instance, the elastic constants and the bulk modulus do not depend on temperature, the constant-pressure and constant-volume heat capacities are equal, and the is null. There are several methods that can add the contribution of volume (pressure) to the thermodynamics of a solid system, among which one of the most powerful is the QHA (Anderson, 1995). The QHA includes an explicit dependence of the vibrational phonons on the crystal volume, *i.e.* , in the harmonic description of the Helmholtz free energy:

The Helmholtz free energy *F*^{QHA}(*V*, *T*) is the sum of the static energy of the system *U*_{0}(*V*) at *T* = 0 K and the vibrational (thermal) contribution . *U*_{0}(*V*) is obtained by any quantum mechanical code by geometry optimization of the at selected (and constrained) volumes. The thermal contribution in the QHA term is defined as

where is the zero-point vibrational energy. From equation (9), it is possible to calculate the equilibrium volume at selected temperatures by minimizing the term with respect to volume. The volumetric coefficient α_{V}(*T*) at selected pressure can be expressed as

It is possible to describe the isothermal bulk modulus (*K _{T}*) of the crystal from the energy second derivative of equation (8) at fixed temperature as

and also the adiabatic bulk modulus (*K _{S}*), which is a preferred way to report the elastic behavior of the solid when comparing the theoretical results with experimental data obtained from some techniques (

*e.g.*elastic wave analysis), as

The great advantage of the QHA approach is the combination of pressure and temperature effects, as the pressure is calculated by

and knowing the temperature we can calculate all the other properties at selected *P*–*T* conditions.

From these formulations, it is possible to calculate the isobaric and (11) and from *V*(*T*):

Finally, other thermodynamic properties could be calculated from the equations above, such as

and Gibbs free energy.The quasi-harmonic treatment, as reported above, needs an adequate knowledge of the phonon dispersion relations, namely how the phonon modes vary with the **k** point in the first This is particularly relevant to properly describe the three acoustic phonon bands of any under consideration, because they represent the main contribution to thermodynamic properties at low temperatures (Belmonte, 2017), but they are always null at the Γ point (central zone, **k** = 0). For simple systems, *e.g.* containing from two to 20 atoms, phonon dispersion relations can be calculated with approaches as described by Parlinski *et al.* (1997). However, such methods could become too computationally demanding for larger unit cells, and hence more approximate algorithms must be employed instead. From this perspective, it is possible to calculate the acoustic thermodynamic properties from sine wave dispersion relations as described by Kieffer (1979*a*,*b*,*c*):

where *X* = *hν*/(*k*_{B}*T*), *R* is the *Z* is the number of unit formulae per and the integrals are evaluated up to the acoustic phonon boundary, the latter deriving, for example, from second-order elastic moduli. Albeit approximate, this approach can provide the acoustic contributions to thermodynamic and thermoelastic properties with adequate accuracy for a system with more than 20–30 atoms in the (Prencipe *et al.*, 2011; Ulian & Valdrè, 2015*b*; Ulian *et al.*, 2021).

#### 4.2. Implementation

*QUANTAS* implements both HA and QHA methods for the analysis of thermodynamics, as summarized in Fig. 1. A scheme of the workflow specifically related to the quasi-harmonic approximation is provided in Fig. 2. Formulations based on both (i) the complete phonon dispersion relations and (ii) the approximated sine wave dispersion relations are included, which can be selected according to the available input data.

While the harmonic approximation routines are quite straightforward, *i.e.* they involve summations of quantities calculated over each phonon frequency for each unit-cell volume (see above), the QHA approach needs special care, as it deals with both fitting and minimization procedures (see Fig. 2). There are two critical points in the quasi-harmonic treatment:

(i) how the phonon/thermodynamics dependence on volume is described;

(ii) how the minimization of the Helmholtz free energy [] curves is performed.

##### 4.2.1. Phonon frequency dependence on volume

For the phonon frequency dependence on volume, two possible strategies are implemented: (i) exploiting an explicit dependence of the phonon modes on unit-cell volume and (ii) providing an implicit dependence using the harmonic approximation thermodynamics. In both cases, linear, quadratic or cubic polynomial functions can be employed. A functional description of the phonon versus unit-cell volume curves, *ν _{i}*(

**k**,

*V*), should be the preferable option, as it allows the thermodynamics of the system to be calculated at each desired

*P*–

*T*–

*V*condition from statistical thermodynamics.

However, to obtain the *ν _{i}*(

**k**,

*V*) curves a correct description of the phonon frequency continuity over the explored volumes is required, taking into account possible crossing of the phonon bands by performing scalar products of the normal mode eigenvectors. In contrast, the employment of a functional (polynomial) form of the thermodynamic quantities over the volume (at constant

*T*) is more straightforward than the former method (since vibrational frequency continuity is not necessary), but there are more fitting procedures involved (for each considered temperature, there are six polynomial fits). Our tests, and the results reported by Erba and co-workers (Erba, 2014; Erba

*et al.*, 2015), suggest that the first method using the phonon continuity over volume provides results in very good agreement with experimental data, depending on the quantum mechanical approach employed to obtain the electronic energy and the vibrational frequencies. However, in

*QUANTAS*both formulations were included, because some users may not be able to obtain (or use) the phonon mode eigenvectors to check their continuity over volume.

##### 4.2.2. Minimization of Helmholtz free energy

Two strategies were implemented in *QUANTAS* for the minimization of the Helmholtz free energy at each temperature value: (i) an equation-of-state fitting procedure of by means of volume-integrated formulations or (ii) minimizing an *n*th-order polynomial that fits as a function of volume. This procedure leads to the equilibrium volume at selected pressure and temperature conditions. Generally, the first method is the preferred experimental way to obtain a phenomenological description of the system under hydrostatic compression, and its main advantage in the *QUANTAS* framework is that it yields with just one fit both the equilibrium volume *V*(*T*) and the bulk modulus *K _{T}*(

*T*) at a selected temperature. The equation-of-state formulations implemented in

*QUANTAS*are the Murnaghan (1937), third-order Birch–Murnaghan (Birch, 1947), Vinet (Vinet

*et al.*, 1987) and Poirier–Tarantola, also called `natural strain' (Poirier & Tarantola, 1998), all of them in their volume-integrated form (Fu & Ho, 1983; Hebbache & Zemzemi, 2004). However, despite being more experimentally friendly, this approach may result in slight numerical noise at high temperature (higher than about 1500 K) in the thermomechanical data, whereas the adoption of a numerical (polynomial fitting) approach can provide more stable results. With polynomial fitting, the isothermal bulk modulus

*K*(

_{T}*T*) can be computed according to equation (11).

The equilibrium volume at each desired pressure is obtained, either from EoSs or from polynomial functions, by minimizing

where *P*_{0} is the target pressure. The *V*(*P*, *T*) values are then used to calculate the thermodynamic and thermoelastic properties at any selected *P**–T* conditions.

#### 4.3. Input data for (Q)HA calculations

In order to use the routine proposed in the present work, some data from quantum mechanical simulations have to be obtained. *QUANTAS* is developed to be `code blind', meaning that any QM code that uses its specific algorithms and measurement units could be used to produce the required input information.

The following input data from QM simulations are required:

(1) Starting bulk geometry (with volume *V*_{e}) of the system under analysis, fully optimized for both lattice parameters and atomic coordinates.

(2) Several bulk structures relaxed in both compression and expansion regimes with respect to the equilibrium volume *V*_{e}. A suitable approach is described in very recent literature (Erba *et al.*, 2015), which considers a number of volumes *N _{V}* between

*V*

_{e}−

*sV*

_{e}% and

*V*

_{e}+ 2

*V*

_{e}

*s*% (including

*V*

_{e}), with

*s*the step of the compression/expansion. Reliable values of

*N*are 4, 7 and 13, and the step

*s*can be set in the range 2–5.

(3) A complete set of vibrational frequencies (or phonon dispersion relations) for each geometrically relaxed unit-cell volume.

In the input file, the unit-cell volume and the lattice static energy (0 K, no thermal contributions) are provided as arrays of length *m* (the number of compressed/expanded volumes), and the phonon modes are reported as a matrix of shape *k* × *m* × *p*, where *k* is the number of sampled **k** points in the Brillouin space and *p* is the number of phonon modes (3*n*, with *n* the number of atoms in the unit cell). For example, the user could employ the Γ-point vibrational frequencies calculated for a containing 40 atoms, studied in five compression states: in this case, *k =* 1, *m =* 5 and *p* = 120, whereas *k* > 1 when phonon dispersion relations were simulated. Compared with the use of only Γ-point vibrational frequencies, the accuracy of the quasi-harmonic approximation results increases when phonon dispersion relations (*k* > 1) are employed as input data.

### 5. Equation of state

In both experimental and theoretical settings, the volumetric behavior of a solid phase with pressure can be described in a functional form called the `equation of state'. The volume of the solid is related to its ; Angel, 2000), and here only the relevant information is discussed.

which could be obtained from high-pressure diffraction experiments. The equation of state is a parametrized function, containing from two to four parameters that are adjusted to fit the experimental data. There are several places in the literature where a detailed description of the theory behind the EoS formulations is provided (Anderson, 1995At the moment, only five isothermal equation-of-state formulations are coded in *QUANTAS*:

(2) Birch–Murnaghan (Birch, 1947):

(3) Vinet (Vinet *et al.*, 1987):

(4) Modified Tait (Freund & Ingalls, 1989):

(5) Natural strain (Poirier & Tarantola, 1998):

Here, *V*_{0T} is the unit-cell volume, *K*_{0T} is the bulk modulus, is the first derivative of the bulk modulus with respect to pressure and is the second derivative of *K*_{0T}. The subscripts `0' and `*T*' mean that each parameter is obtained at reference pressure zero and temperature of interest *T*, respectively.

The modified Tait (T) and Murnaghan (M) EoSs are `invertible' formulations, as it is possible to express the unit-cell volume as a function of pressure by inverting the equation. In addition, the modified Tait equation of state can be reduced to the Murnaghan one by imposing = 0.

The Birch–Murnaghan (BM) and natural strain (NS) equations of state are `finite-strain EoSs', which were formulated considering that the energy of the compressed solid can be expressed as a Taylor series in the *f* (that is Eulerian strain, *f*_{E}, for the BM EoS and natural strain, *f*_{N}, for the NS EoS). Both of them are fourth-order expansions, but they can be truncated to third- and second-order expressions by using assumed values for the and parameters, respectively, as described by Angel *et al.* (2014).

The Vinet (V) equation of state was derived from molecular mechanics models for very high compression regimes and is a third-order EoS.

According to the literature (Angel, 2000) and to the reported formulations, the fitting strategy considers the volume as the independent variable and the pressure as the dependent one, because the experimental uncertainties on *V* are generally much lower than those on *P*. Then, a least-squares method is used to fit the data, employing the errors on the variables as weights during the procedure. For example, this is the approach adopted by the well known *EosFit* software (Angel, 2000; Angel *et al.*, 2014). The goodness of fit is given by the residual variance (weighted chi squared, χ^{2}), which is equal to unity if the EoS model perfectly matches the weighted experimental data. In contrast, if χ^{2} > 1 it means that the equation of state correctly represents only a portion of the data for several possible reasons, *e.g.* some compression states were not adequately obtained, the errors of the values were underestimated or the model is not accurate enough to describe all of the data set. For example, it is discouraged to use the Murnaghan EoS for unit-cell compressions higher than 10%. A value of χ^{2} < 1 does not represent a better fit and may also express an overfitting of the data, namely the equation-of-state model contains more parameters than the number of *P*–*V* points.

A visual assessment of the goodness of fit can also be obtained from the strain-normalized pressure plots (*f*–*F*), which are usually clearer than the standard *P*–*V* plots.

#### 5.1. Implementation

While the physics behind the different formulations is well known within the community, the fitting procedure requires some explanation. Standard least-square procedures consider uncertainties only on the dependent variable, while the independent one is considered free from errors. In the EoS fitting, the dependent variable is pressure because the errors associated with pressure measurements are usually much larger than those of the unit-cell/lattice parameters. However, including the uncertainties on the unit-cell volume *V* would increase the accuracy of the fitted EoS parameters. Differently from *EosFit7* (Angel *et al.*, 2014), which employs the effective variance method described by Orear (1982), *QUANTAS* adopts the orthogonal distance regression (ODR) statistical approach (Boggs *et al.*, 1987), which allows the inclusion of the uncertainties of both dependent and independent variables. The ODR algorithm is coded in the *ORDPACK* library (Boggs *et al.*, 1989; Zwolak *et al.*, 2007), which is included in the *SciPy* Python package used by *QUANTAS*. If the user desires to weight the fit only for pressure, the software automatically switches to the ordinary least-squares procedure.

The fitting procedure is completely interactive, and users can have exact and complete control over it. As in *EosFit7* (Angel *et al.*, 2014), users may set weights and fix parameters and modify them to find the best EoS parameters for their data.

For this operation, the input file must contain the isothermal data of the material

at different pressure states. These quantities can derive from either experiments or theoretical simulations. The input data are organized in a table-like format, indicating which quantity is present in each column. Uncertainties on both pressure and unit-cell volume/lattice parameters can be included.### 6. Elastic properties from second-order elastic constants

In the theory of linear elasticity, the stress tensor can be expressed in terms of strain by

*C _{ijkl}* are the components of the fourth-rank modulus tensor, whose coordinates depend on the choice of the axes. Equation (23) can be inverted, leading to

where *S _{ijkl}* are the components of the compliance tensor, the inverse modulus tensor. A fourth-rank tensor has 81 components, but a maximum of 21 independent values, for triclinic crystals. It is also common to express the stiffness and compliance tensor components using the Voigt (engineering) matrix notation (see Nye, 1957), which is often adopted because of its simplicity. Note that the Voigt notation is not a tensor, but a matrix representation of it.

According to Nye (1957), several mechanical properties related to the polycrystalline (isotropic) elastic behavior of a material can be calculated from the elastic tensor by means of the Voigt, Reuss and Hill averages:

where *E*_{VRH} is Young's modulus, and *K*_{R}, *K*_{V}, *K*_{VRH}, μ_{R}, μ_{V} and μ_{VRH} are the Voigt, Reuss and Hill values of the bulk and shear moduli, respectively.

The mean shear, , and longitudinal, , wave speeds of a polycrystal with no

of the grains depend on the coupling between grains and can range from the Reuss limit (with free grain boundaries) to the Voigt limit (with locked grain boundaries). Most randomly oriented polycrystals have shear and Young's moduli close to, but not identical to, the VRH averages, for which the following approximation of the mean wave speeds is valid:where ρ is the density of the crystal.

The calculation of the six eigenvalues of the second-order elastic tensor allows us to define the mechanical stability of the solid as described by Mouhat & Coudert (2014) for the different crystal systems: if any of the eigenvalues are negative, the system is unstable.

It is also possible to derive single-crystal elastic properties from the stiffness matrix, and an excellent treatment of their calculation was recently proposed by Marmier *et al.* (2010). In brief, this requires the transformation of the stiffness tensor using the following rule for generic fourth-rank tensors ** T**:

which employs the Einstein summation rule with the terms *r _{αi}* being the direction cosines. In a Cartesian reference system, it is possible to represent a direction corresponding to an elastically significant distortion as a point on a unit sphere (unit vector

**a**), using two angles, θ(0, π) and φ(0, 2π):

A single vector is sufficient to calculate the spatial variation of some properties, such as Young's modulus or the linear compressibility, but a second vector **b**, perpendicular to **a**, is required to obtain the and Poisson's ratio. The second vector is characterized by a third angle, χ(0, 2π), and by the coordinates

Then, the coordinates of **a** and **b** represent the first two columns of the rotation matrix, which allows the calculation of all the components in the subvectorial space defined by directions 1 and 2:

The spatial dependence of the elastic modulus *E* and linear compressibility β are defined as

and the spatial variations of the μ and Poisson's ratio υ are given by the following formulae:

Directional variations of seismic (acoustic) wave speeds can be obtained by solving the Christoffel equation, as reported by Musgrave (1970).

#### 6.1. Implementation

The analysis of the second-order elastic constant matrix is performed by a forked version of the *ELATE* script (Gaillac *et al.*, 2016), which was modified in our program to work offline and adapted to the *QUANTAS* framework. This modified version includes also a routine to find the reference and provides analysis of the seismic wave speeds by solving the Christoffel equations as reported by Musgrave (1970). Regarding the plotting capabilities, the original framework was changed to the *Matplotlib* library (Hunter, 2007), providing polar plots of the spatial variations of the elastic properties on the *xy*, *xz* and *yz* planes.

The elastic moduli can be provided to *QUANTAS* as an input file containing the stiffness matrix (*i.e.* the elastic tensor in Voigt's notation) in either full, upper triangular or lower triangular form. The code performs a conversion from the matrix form to the tensorial form for the analysis of the elastic properties. This conversion preserves the Cartesian reference frame that was employed to obtain the elastic tensor in Voigt's notation. If the crystal density is supplied, analysis of the seismic wave speeds is also performed. The analysis is automatic, and the user may ask to produce polar plots of the results.

### 7. Future development

For the future development steps of *QUANTAS*, we are working on new functionalities, for example, the capability to construct the phase diagrams for both a single substance (*e.g.* silica compounds) and solid solutions, by using the calculated Gibbs free energy values. This would resemble the CALPHAD method (Lukas *et al.*, 2007), but applied to the theoretical/experimental materials science fields. Concerning the EoS capabilities, it is planned to add thermal and *P*–*T* and *P*–*T*–*V* equations of state to aid the analysis of high-pressure and high-temperature data.

In addition, a graphical user interface is currently under development, with the aim of easing the use of the software and providing a plug-in platform that could be easily implemented with new functionalities.

Finally, since *QUANTAS* is an open-source project, we encourage users to contribute to our code, extending its functionalities with new modules.

### 8. Conclusions

There is a continuous and growing interest in the development and characterization of different materials at both experimental and theoretical levels, with important basic, technological and industrial applications. In this context, a detailed knowledge of the thermodynamic and elastic stability of solid phases is of utmost importance. In addition, thanks to increasing computational power, density functional theory is rising as a competing tool to drive innovation in materials science, physics, chemistry and other disciplines. Indeed, this framework offers an unprecedented balance between accuracy and speed.

The scope of *QUANTAS* is to provide a tool to reach a better understanding of both experimental and theoretical results, and to extend the knowledge on the *P*–*T* behavior of synthetic materials and minerals. *QUANTAS* is a fast easy-to-use software that can support researchers in speeding up their calculations, improving data quality over a wide range of cases. The software structure is simple and written in the Python programming language, which is truly cross platform.

As developers, we also encourage and support the integration of *QUANTAS* into other software that relies on quantum mechanical simulations of solid phases. The source code is flexible and reusable with little modification by other developers, a feature that is also facilitated by the distribution of the code under the Simplified BSD license.

Last, but not least, *QUANTAS* can be considered as a useful teaching software aid as well as a practical and powerful research tool. Students of various levels of study may find it helpful in understanding how thermodynamics and mechanical properties of materials behave and evolve.

### 9. Related literature

The following additional literature is cited in the supporting information: Anderson *et al.* (1991, 1992, 1993), Angel *et al.* (1997), Chopelas (1990), Comodi *et al.* (2002), Cynn *et al.* (1995), Gatta *et al.* (2006, 2014), Isaak *et al.* (1989), Kantor *et al.* (2012), Kresse & Joubert (1999), Lotti *et al.* (2017), Monkhorst & Pack (1976), Pascale *et al.* (2004), Peintinger *et al.* (2013), Perdew *et al.* (1996), Perger *et al.* (2009) and Tosoni *et al.* (2005).

### Supporting information

Explained and discussed test cases for the software. DOI: https://doi.org/10.1107/S1600576722000085/vb5027sup1.pdf

Archive containing the software and the related documentation (part 1). DOI: https://doi.org/10.1107/S1600576722000085/vb5027sup2.bin

Archive containing the software and the related documentation (part 2). DOI: https://doi.org/10.1107/S1600576722000085/vb5027sup3.bin

Archive containing the software and the related documentation (part 3). DOI: https://doi.org/10.1107/S1600576722000085/vb5027sup4.bin

Archive containing the software and the related documentation (part 4). DOI: https://doi.org/10.1107/S1600576722000085/vb5027sup5.bin

Archive containing the software and the related documentation (part 5). DOI: https://doi.org/10.1107/S1600576722000085/vb5027sup6.bin

### Acknowledgements

The authors thank the beta testers of *QUANTAS* for their feedback on the program. Open access funding provided by Universita degli Studi di Bologna within the CRUI-CARE Agreement.

### References

Anderson, O. L. (1995). *Equations of State of Solids for Geophysics and Ceramic Science.* New York: Oxford University Press. Google Scholar

Anderson, O. L., Isaak, D. L. & Oda, H. (1991). *J. Geophys. Res.* **96**, 18037–18046. CrossRef CAS Web of Science Google Scholar

Anderson, O. L., Oda, H., Chopelas, A. & Isaak, D. G. (1993). *Phys. Chem. Miner.* **19**, 369–380. CrossRef CAS Google Scholar

Anderson, O. L., Oda, H. & Isaak, D. (1992). *Geophys. Res. Lett.* **19**, 1987–1990. CrossRef Web of Science Google Scholar

Angel, R. J. (2000). *Rev. Mineral. Geochem.* **41**, 35–59. Web of Science CrossRef Google Scholar

Angel, R. J., Allan, D. R., Miletich, R. & Finger, L. W. (1997). *J. Appl. Cryst.* **30**, 461–466. CrossRef CAS Web of Science IUCr Journals Google Scholar

Angel, R. J., Alvaro, M. & Gonzalez-Platas, J. (2014). *Z. Kristallogr.* **229**, 405–419. Web of Science CrossRef CAS Google Scholar

Baroni, S., Giannozzi, P. & Isaev, E. (2010). *Rev. Mineral. Geochem.* **71**, 39–57. Web of Science CrossRef CAS Google Scholar

Behnel, S., Bradshaw, R., Citro, C., Dalcin, L., Seljebotn, D. S. & Smith, K. (2011). *Comput. Sci. Eng.* **13**, 31–39. Web of Science CrossRef Google Scholar

Belmonte, D. (2017). *Minerals* **7**, 183. Web of Science CrossRef Google Scholar

Birch, F. (1947). *Phys. Rev.* **71**, 809–824. CrossRef CAS Web of Science Google Scholar

Boggs, P. T., Byrd, R. H. & Schnabel, R. B. (1987). *SIAM J. Sci. Stat. Comput.* **8**, 1052–1078. CrossRef Web of Science Google Scholar

Boggs, P. T., Donaldson, J. R., Byrd, R. H. & Schnabel, R. B. (1989). *ACM Trans. Math. Softw.* **15**, 348–364. CrossRef Web of Science Google Scholar

Chopelas, A. (1990). *Phys. Chem. Miner.* **17**, 142–148. CAS Google Scholar

Comodi, P., Gatta, G. D., Zanazzi, P. F., Levy, D. & Crichton, W. (2002). *Phys. Chem. Miner.* **29**, 538–544. Web of Science CrossRef CAS Google Scholar

Cynn, H., Anderson, O. L., Isaak, D. G. & Nicol, M. (1995). *J. Phys. Chem.* **99**, 7813–7818. CrossRef CAS Web of Science Google Scholar

Destefanis, M., Ravoux, C., Cossard, A. & Erba, A. (2019). *Minerals* **9**, 16. Web of Science CrossRef Google Scholar

Dovesi, R., Erba, A., Orlando, R., Zicovich-Wilson, C. M., Civalleri, B., Maschio, L., Rérat, M., Casassa, S., Baima, J., Salustro, S. & Kirtman, B. (2018). *WIREs Comput. Mol. Sci.* **8**, e1360. Google Scholar

Duffy, T. S., Ahrens, T. J. & Lange, M. A. (1991). *J. Geophys. Res.* **96**, 14319–14330. CrossRef CAS Web of Science Google Scholar

Erba, A. (2014). *J. Chem. Phys.* **141**, 124115. Web of Science CrossRef PubMed Google Scholar

Erba, A., Maul, J., Demichelis, R. & Dovesi, R. (2015). *Phys. Chem. Chem. Phys.* **17**, 11670–11677. Web of Science CrossRef CAS PubMed Google Scholar

Freund, J. & Ingalls, R. (1989). *J. Phys. Chem. Solids*, **50**, 263–268. CrossRef CAS Web of Science Google Scholar

Fu, C. L. & Ho, K. M. (1983). *Phys. Rev. B*, **28**, 5480–5486. CrossRef CAS Web of Science Google Scholar

Fukui, H., Ohtaka, O., Suzuki, T. & Funakoshi, K. (2003). *Phys. Chem. Miner.* **30**, 511–516. Web of Science CrossRef CAS Google Scholar

Gaillac, R., Pullumbi, P. & Coudert, F. X. (2016). *J. Phys. Condens. Matter*, **28**, 275201. Web of Science CrossRef PubMed Google Scholar

Gatta, G. D., Lotti, P., Merlini, M., Liermann, H.-P., Lausi, A., Valdrè, G. & Pavese, A. (2015). *Phys. Chem. Miner.* **42**, 309–318. Web of Science CrossRef CAS Google Scholar

Gatta, G. D., Merlini, M., Valdrè, G., Liermann, H.-P., Nénert, G., Rothkirch, A., Kahlenberg, V. & Pavese, A. (2013). *Phys. Chem. Miner.* **40**, 145–156. Web of Science CrossRef ICSD CAS Google Scholar

Gatta, G. D., Morgenroth, W., Dera, P., Petitgirard, S. & Liermann, H. P. (2014). *Phys. Chem. Miner.* **41**, 569–577. Web of Science CrossRef ICSD CAS Google Scholar

Gatta, G. D., Nestola, F. & Ballaran, T. B. (2006). *Phys. Chem. Miner.* **33**, 235–242. Web of Science CrossRef ICSD CAS Google Scholar

Golesorkhtabar, R., Pavone, P., Spitaler, J., Puschnig, P. & Draxl, C. (2013). *Comput. Phys. Commun.* **184**, 1861–1873. Web of Science CrossRef CAS Google Scholar

Gonzalez-Platas, J., Alvaro, M., Nestola, F. & Angel, R. (2016). *J. Appl. Cryst.* **49**, 1377–1382. Web of Science CrossRef CAS IUCr Journals Google Scholar

Hebbache, M. & Zemzemi, M. (2004). *Phys. Rev. B*, **70**, 224107. Web of Science CrossRef Google Scholar

Hunter, J. D. (2007). *Comput. Sci. Eng.* **9**, 90–95. Web of Science CrossRef Google Scholar

Isaak, D. G., Anderson, O. L. & Goto, T. (1989). *Phys. Chem. Miner.* **16**, 704–713. CrossRef CAS Google Scholar

Jiang, F. M., Speziale, S. & Duffy, T. S. (2006). *Am. Mineral.* **91**, 1893–1900. Web of Science CrossRef CAS Google Scholar

Joppa, L. N., McInerny, G., Harper, R., Salido, L., Takeda, K., O'Hara, K., Gavaghan, D. & Emmott, S. (2013). *Science*, **340**, 814–815. Web of Science CrossRef CAS PubMed Google Scholar

Kantor, A., Kantor, I., Merlini, M., Glazyrin, K., Prescher, C., Hanfland, M. & Dubrovinsky, L. (2012). *Am. Mineral.* **97**, 1764–1770. Web of Science CrossRef ICSD CAS Google Scholar

Kieffer, S. W. (1979*a*). *Rev. Geophys.* **17**, 1–19. CrossRef CAS Web of Science Google Scholar

Kieffer, S. W. (1979*b*). *Rev. Geophys.* **17**, 20–34. CrossRef Web of Science Google Scholar

Kieffer, S. W. (1979*c*). *Rev. Geophys.* **17**, 35–59. CrossRef CAS Web of Science Google Scholar

Kresse, G. & Furthmüller, J. (1996). *Comput. Mater. Sci.* **6**, 15–50. CrossRef CAS Web of Science Google Scholar

Kresse, G. & Hafner, J. (1993). *Phys. Rev. B*, **48**, 13115–13118. CrossRef CAS Web of Science Google Scholar

Kresse, G. & Joubert, D. (1999). *Phys. Rev. B*, **59**, 1758–1775. Web of Science CrossRef CAS Google Scholar

Lotti, P., Gatta, G. D., Comboni, D., Guastella, G., Merlini, M., Guastoni, A. & Liermann, H. P. (2017). *J. Am. Ceram. Soc.* **100**, 2209–2220. Web of Science CrossRef CAS Google Scholar

Lukas, H. L., Fries, S. G. & Sundman, B. (2007). *Computational Thermodynamics: the CALPHAD Method.* Cambridge University Press. Google Scholar

Marmier, A., Lethbridge, Z. A. D., Walton, R. I., Smith, C. W., Parker, S. C. & Evans, K. E. (2010). *Comput. Phys. Commun.* **181**, 2102–2115. Web of Science CrossRef CAS Google Scholar

Milani, S., Angel, R. J., Scandolo, L., Mazzucchelli, M. L., Ballaran, T. B., Klemme, S., Domeneghetti, M. C., Miletich, R., Scheidl, K. S., Derzsi, M., Tokár, K., Prencipe, M., Alvaro, M. & Nestola, F. (2017). *Am. Mineral.* **102**, 851–859. Web of Science CrossRef Google Scholar

Monkhorst, H. J. & Pack, J. D. (1976). *Phys. Rev. B*, **13**, 5188–5192. CrossRef Web of Science Google Scholar

Mouhat, F. & Coudert, F. X. (2014). *Phys. Rev. B*, **90**, 224104. Web of Science CrossRef Google Scholar

Murnaghan, F. D. (1937). *Am. J. Math.* **59**, 235–260. CrossRef Google Scholar

Murri, M., Mazzucchelli, M. L., Campomenosi, N., Korsakov, A. V., Prencipe, M., Mihailova, B. D., Scambelluri, M., Angel, R. J. & Alvaro, M. (2018). *Am. Mineral.* **103**, 1869–1872. Google Scholar

Musgrave, M. J. P. (1970). *Crystal Acoustics: Introduction to the Study of Elastic Waves and Vibrations in Crystals.* San Francisco: Holden-Day. Google Scholar

Nye, J. F. (1957). *Physical Properties of Crystals.* Oxford University Press. Google Scholar

Orear, J. (1982). *Am. J. Phys.* **50**, 912–916. CrossRef Web of Science Google Scholar

Ottonello, G., Civalleri, B., Ganguly, J., Perger, W. F., Belmonte, D. & Vetuschi Zuccolini, M. (2010). *Am. Mineral.* **95**, 563–573. Web of Science CrossRef CAS Google Scholar

Ottonello, G., Civalleri, B., Ganguly, J., Vetuschi Zuccolini, M. & Noel, Y. (2009). *Phys. Chem. Miner.* **36**, 87–106. Web of Science CrossRef CAS Google Scholar

Ottonello, G., Vetuschi Zuccolini, M. & Civalleri, B. (2009). *Calphad*, **33**, 457–468. Web of Science CrossRef CAS Google Scholar

Parlinski, K., Li, Z. Q. & Kawazoe, Y. (1997). *Phys. Rev. Lett.* **78**, 4063–4066. CrossRef CAS Web of Science Google Scholar

Pascale, F., Zicovich-Wilson, C. M., López Gejo, F., Civalleri, B., Orlando, R. & Dovesi, R. (2004). *J. Comput. Chem.* **25**, 888–897. Web of Science CrossRef PubMed CAS Google Scholar

Pawley, A. R., Redfern, S. A. T. & Wood, B. J. (1995). *Contrib. Mineral. Petrol.* **122**, 301–307. CrossRef CAS Web of Science Google Scholar

Peintinger, M. F., Oliveira, D. V. & Bredow, T. (2013). *J. Comput. Chem.* **34**, 451–459. Web of Science CrossRef CAS PubMed Google Scholar

Perdew, J. P., Burke, K. & Ernzerhof, M. (1996). *Phys. Rev. Lett.* **77**, 3865–3868. CrossRef PubMed CAS Web of Science Google Scholar

Perger, W. F. (2010). *Int. J. Quantum Chem.* **110**, 1916–1922. Web of Science CrossRef CAS Google Scholar

Perger, W. F., Criswell, J., Civalleri, B. & Dovesi, R. (2009). *Comput. Phys. Commun.* **180**, 1753–1759. Web of Science CrossRef CAS Google Scholar

Poirier, J. P. & Tarantola, A. (1998). *Phys. Earth Planet. Inter.* **109**, 1–8. Web of Science CrossRef Google Scholar

Prencipe, M., Pascale, F., Zicovich-Wilson, C. M., Saunders, V. R., Orlando, R. & Dovesi, R. (2004). *Phys. Chem. Miner.* **31**, 559–564. Web of Science CrossRef ICSD CAS Google Scholar

Prencipe, M., Scanavino, I., Nestola, F., Merlini, M., Civalleri, B., Bruno, M. & Dovesi, R. (2011). *Phys. Chem. Miner.* **38**, 223–239. Web of Science CrossRef CAS Google Scholar

Rodríguez-Carvajal, J. & González-Platas, J. (2005). *Acta Cryst.* A**61**, C22. Google Scholar

Rodriguez-Carvajal, J. & Gonzalez-Platas, J. (2008). *Acta Cryst.* A**64**, C46–C46. Google Scholar

Togo, A. & Tanaka, I. (2015). *Scr. Mater.* **108**, 1–5. Web of Science CrossRef CAS Google Scholar

Tosoni, S., Pascale, F., Ugliengo, P., Orlando, R., Saunders, V. R. & Dovesi, R. (2005). *Mol. Phys.* **103**, 2549–2558. CAS Google Scholar

Ulian, G., Moro, D. & Valdrè, G. (2021). *Am. Mineral.* **106**, 1928–1939. Web of Science CrossRef Google Scholar

Ulian, G. & Valdrè, G. (2015*a*). *Am. Mineral.* **100**, 935–944. Web of Science CrossRef Google Scholar

Ulian, G. & Valdrè, G. (2015*b*). *Phys. Chem. Miner.* **42**, 151–162. Web of Science CrossRef CAS Google Scholar

Ulian, G. & Valdrè, G. (2015*c*). *Phys. Chem. Miner.* **42**, 609–627. Web of Science CrossRef CAS Google Scholar

Ulian, G. & Valdrè, G. (2018). *J. Mech. Behav. Biomed. Mater.* **77**, 683–692. Web of Science CrossRef CAS PubMed Google Scholar

Ulian, G. & Valdrè, G. (2019). *Acta Cryst.* B**75**, 1042–1059. Web of Science CrossRef IUCr Journals Google Scholar

Vinet, P., Ferrante, J., Rose, J. H. & Smith, J. R. (1987). *J. Geophys. Res.* **92**, 9319–9325. CrossRef CAS Web of Science Google Scholar

Walt, S. van der, Colbert, S. C. & Varoquaux, G. (2011). *Comput. Sci. Eng.* **13**, 22–30. Google Scholar

Wilson, G., Aruliah, D. A., Brown, C. T., Chue Hong, N. P., Davis, M., Guy, R. T., Haddock, S. H. D., Huff, K. D., Mitchell, I. M., Plumbley, M. D., Waugh, B., White, E. P. & Wilson, P. (2014). *PLoS Biol.* **12**, e1001745. Web of Science CrossRef PubMed Google Scholar

Yoon, H. S. & Lawrence Katz, J. (1976). *J. Biomech.* **9**, 459–464. CrossRef CAS PubMed Web of Science Google Scholar

Zwolak, J. W., Boggs, P. T. & Watson, L. T. (2007). *Acm T Math. Softw.* **33**, 27. Web of Science CrossRef Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.