Main features of the program and more important functions implemented

The program is useful for the post-processing of ab initio results concerning static energies and vibrational frequencies as functions of the unit cell volume of a crystal. Within the limit of the quasi-harmonic approximation (QHA) and within the framework of statistical thermodynamics, the Helmholtz free energy ($F$) is computed as a function of the volume of the unit cell and of the temperature; by knowing $F$ and by means of classical thermodynamics relations, a number of properties are then evaluated. These include:

Elastic properties and thermal expansion

  1. Equation of state (bulk modulus, its pressure derivative and its temperature dependence, equilibrium volume) at the static condition and in a range of temperatures. Currently, the 3^rd and 4^st order Birch-Murnaghan P(V) equations are implemented (BM3 and BM4 respectively);
  2. Thermal expansion as a function of temperature and pressure.

Thermodynamics properties

  1. Specific heat at costant volume and at constant pressure;
  2. Entropy;
  3. Helmholtz free energy;
  4. Gibbs free energy;
  5. Phase transitions and mineral reactions: equilibrium curves in the P/T space. Currently the program does not consider solid solutions; only end-member mineral phases can be considered.

Frequencies analysis

  1. Analysis of the variations of vibrational frequencies with pressure, volume or temperature of the crystal.

Phonon dispersion

Concerning dispersion relations in the reciprocal space, and their impact of the calculated properties, unless super-cell computations of the frequencies are available, the modified Kieffer model (see this paper and the references therein) can be used to evaluate the contribution to the thermodynamics properties of the acoustic vibrational modes.

In case supercell calculations are available, the frequencies of the off-center vibrational modes, together with their volume dependence, can be taken into account. See the tutorial here.

LO-TO splitting

The impact of the LO-TO splitting on properties can also be evaluated. See the tutorial here.

Anharmonicity

The program anharm.py of the package can be used to take into account anharmonic effects on selected vibrational modes, provided that scans of the anharmonic potential as functions of the corresponding normal coordinates of the modes are available, at different values of the unit cell volume. See the relevant tutorial here.

The required data (scans of the anharmonic potential) can be obtained by using the CRYSTAL17 code by means of the SCANMODE keyword, after a standard frequency calculation has been performed.

Information concerning the required input files can be found in the Basic EoS Tutorial.

The examples in the documentation below are executed with reference to the case of pyrope. The program is launched and the input files are loaded. The input.txt file is shown as soon as the program starts (and a static BM3 EoS is computed).

Dealing with the frequencies

In the calculation of the equation of state it is possible to choose among several options; the first two choices are:

  1. the optimization of a volume-integrated EoS (V-BM3, or V-BM4) on the F(V) values computed at a given temperature
  2. the optimization of the EoS (BM3, or BM4) on the P(V) values computed at a given temperature.

Concerning the choice (1), there are three possibilities:

In the current example, a polynomial fit is requested:

FITVOL
POLY
725. 769. 16 3

The FITVOL activates the frequency fit; the POLY keyword below specifies a polynomial fit. The first two float numbers following the POLY keyword are the minimum and maximum volumes for the fit; the integer 16 in the number of points, in the volume range, at which the frequencies will be obtained (from the fit) in order to calculate the Helmholtz energy; the EoS will then be optimized on this set of 16 F(V) points.

The command for the EoS optimization is eos_temp(tt), where the argument tt is the temperature:

The other possible choice for the frequency fitting is the spline fit. This can be required in input.txt with:

FITVOL
SPLINE
725. 769. 16 3 1.

The last float number is the smoothness parameter. Information about the spline fit can be found here; look at the documentation concerning the s parameter for the meaning of smoothness.

The polynomial or spline fittings can be requested at any time, overriding the choice made in input.txt by the commands

set_poly(deg, npoint)
set_spline(deg, smooth, npoint)

default values for the parameters are:

As an example:

Note: The volume range originally requested in input.txt was $[725., 769. A^3]$; by requesting the spline fit, by the command above, the volume range changed. This behaviour is due to the fact that the volume ranges for the different type of fittings are allowed to be different: as no volume range was specified for the spline fitting, the default choice is to take the minimum and maximum volumes at which the frequencies have been computed, with an applied offset of 0.2 A^3. Such volumes can be visualized by the command info.show:

To change the volume range, the command set_volume_range can be issued; as we currently are in the spline fitting mode, the command will act on the volume range for the spline fit:

Concerning the command above, the optional parameter npoint (default=16) is available; the default value for prt is False: the volume range is changed without any printed notification.

The command eos_temp will now output the equation of state (BM3) by using a spline fit of the frequencies:

To deactivate any fitting option (POLY or SPLINE of whatever degree), the FITVOL keyword can be removed from the input file, or the command fit_off can be issued:

To see the effect on the equation of state:

A wrote above (option 2), the EoS can also be optimized by fitting a BM3 (or BM4) function on a P(V) set of data. In this case, pressures at each volume are computed as the partial derivative of F with respect to V, at constant temperature:

$$P=-\left(\frac{\partial F}{\partial V}\right)_{\!T}$$

To this end, some kind of fitting of the frequencies is mandatory. In fact, as we switched off fitting, by requesting the EoS P(V) by the bulk_dir command, we get:

where the argument "300" is the temperature. We have to switch on the frequency fitting, and then we can do the EoS optimization; for instance:

A number of commands are available to show the behaviour of frequencies with the volume and the quality of fits. For instance, the command check_poly_list will make the analysis on a specified set of modes:

Take a look at this tutorial to see other commands of the same type.

Equation of state

As already seen above, the two commands for optimizing the equation of state at a given temperature are:

the command static can be issued to get the static equation of state (3^rd order only): this is based on the static energy values only (no vibrational effects included):

A 4^th order Birch-Murnaghan EoS can also be optimized; see the relevant section in basic_eos_tutorial for details.

Fixing Kp

For BM3 optimization (with no fit of frequencies, as well as with all the possible fit types), the first derivative of the bulk modulus with respect to the pressure (Kp) can be fixed. This can be achieved by issuing the command set_fix(kp) whose argument is the fixed value given for Kp (default: 4). For instance:

Note that some commands can override the directive set_fix by providing the optional keyword argument fix: precisely, if fix=0., Kp is optimized:

The command reset_fix deactivates the set_fix directive:

Note however that in:

the keyword fix (not set equal to 0.) forces Kp to be fixed at the given value (of 4.5). Other commands behave in a similar way; however eos_temp always optimizes Kp unless it is fixed with set_fix.

Note Due to correlations among K0, Kp and V0, it is generally advised to keep Kp fixed when the bulk modulus is used to compute other properties depending on its variations with P, V or T. For instance:

In particular, the function bulk_serie is used to compute the K0 dependence upon T (in a given T range). The $K_0(T)$ values are then fitted by a polynomial of a given degree; in the example, the T range is the linear one and $dK_0/dT = -0.0371$, with Kp fixed, or $-0.0348$ if Kp is not fixed.

Thermal expansion

The command thermal_exp_p(t,p) can be used to compute the thermal expansion at a given temperature (t) and pressure (p). Note in the example that we a still under the set_fix directive issued above, as demonstrated by the fix_status command:

To inquire about the current fit status of the frequencies, the command fit_status must be given:

The thermal expansion at T = 300K and P = 0 GPa is:

In this calculation Kp (which is an intermediate result of the algorithm) can be optimized by using the optional keyword fix:

Thermal expansion (at constant pressure) in a temperature range and at a given pressure, can be computed by the function alpha_serie. In the example, $\alpha(T)$ is computed in the range $[20, 600K]$ (18 points) at a pressure of 0GPa. The fit optional keyword (default: True) requests a fit of the $\alpha(T)$ values by a polynomial function whose $T$ powers are specified in the input file:

ALPHA
0. -2  -1 -0.5

This means that the thermal expansion is fitted as:

$$\alpha(T) = k + aT^{-2} + bT^{-1} + cT^{-0.5}$$

where $k, a, b, c, $ are the optimized constants.

The array printed at the bottom of the cell contains the optimized parameters (following the order of the corresponding given powers). If the optional keyword prt is set to True (default), a list of the computed $\alpha$ values is provided.

Entropy and specific heat

The function entropy_p computes the entropy at given temperature and pressure; it also provides the specific heat at constant volume ($C_V$). In the example, the entropy is computed at T=300K and P=0GPa:

The specific heat at constant pressure ($C_P$) is provided by the function cp; the temperature and the pressure must be specified; the function also provides $C_V$, $\alpha$ and $K_0$:

A list of values of $C_P$ in a given temperature range can be obtained by the function cp_serie; in the example, we compute $C_P$ in the range $[40, 400K]$ (12 points) at P=0GPa; a fit is requested (optional keyword fit set to True, by default) with a polynomial whose powers are specified in the input file under the keyword CP (same structure as the ALPHA fit described above):

To have a printout of the coefficients of the fitting polynomial, set prt=False to suppress the printing of the $C_P$ values at each temperature in the interval (to suppress the plot, too, use the graph optional keyword):

Such values of the coefficient can be saved in a variable like:

and used to compute $C_P$ at any temperature by using the polynomial fitting function cp_fun:

In the last example, the Python round command has been used to approximate the printed value returned by cp_fun to 2 decimal digits.

LO-TO splitting

In the example above, the presence of the LO keyword in the input.txt file, followed by a file name (LO.txt in this case), requests the LO-TO splitting correction to the frequencies.

The LO.txt file contains two columns: the first one specifies the mode numbers (for those modes affected by LO-TO splitting only); the second column specifies the corresponding values of the split (in $cm^{-1}$).

The presence of the LO keyword in the input file produces corrected frequencies for the modes affected by LO-TO splitting. Precisely, the new frequencies $\nu_c$ are computed as follows

$$\nu_c = 2/3\cdot \nu_{TO} + 1/3\cdot \nu_{LO}$$

where $\nu_{LO}=\nu_{TO}+\Delta\nu$ and $\Delta\nu$ is the LO-TO split value. The TO values are found in the frequencies file specified by the keyword FREQ.

The corrected frequencies $\nu_c$ are internally stored in the class lo, along with a copy of the original TO frequencies, and used in all the subsequent calculations instead of the TO frequencies.

To deactivate the LO-TO correction, either the LO keyword can be deleted from the input file, or the method lo.off can be issued (lo.on re-activates the correction):

To see the effect of the correction on the EoS parameters in this case:

which is to be compared with:

The bulk modulus is not significantly affected by such LO-TO correction; however thermodynamic properties can be changed:

to be compared with:

The Kieffer's model for the acoustic branches

The presence of the KIEFFER keyword in the input.txt file activates the evaluation of the contribution to the $F$ Helmholtz free energy coming from the acoustic phonon branches; such contribution is neglected in a standard calculation of frequencies at the center of the Brillouin zone. The keyword is followed by three numbers representing the frequencies of the acoustic modes at the Brillouin zone boundary (in cm^-1), that can be estimated from the knowledge of the elastic tensor of the crystal.

Apart from being activated by the keyword KIEFFER, the correction can be applied by the methods kieffer.on (to set up the calculation) and kieffer.freq to provide the frequencies (if not provided in the input file, or to change them). The method kieffer.plot can be used to plot the $F_{acoustic}(T)$ function:

At any temperature, the contribution to $F$ from the acoustic branches can be recovered with the method kieffer.get_value; for instance:

Note that the Kieffer's model does not take into account a possible variation of the frequencies of the acoustic modes with the unit cell volume; therefore it has no impact on all those properties depending uniquely upon the derivative of $F$ with respect to $V$; e.g.:

On the other hand, other properties can be significantly affected:

to be compared with:

Thermodynamics and mineral equilibria

Functions exist to compute and to plot the Gibbs free energy $G=F+PV$ as function of $T$ or $P$; for instance:

In the above axample, the Gibbs pree energy of pyrope has been computed at 10 points, in the $[200, 400K]$ temperature range, at a fixed pressure of 0GPa.

The keywords v0 and g0, among the optional arguments of the function above, specify values from other sources of the molar volume and of the Gibbs free energy at the standard conditions. Such values are retrieved from a database and stored as attributes of the variable py of the mineral class. The py.info method prints out all the stored information concerning pyrope (in the present case, from the Holland & Powell database):

To compute the $G$ energy in the range $[0, 10GPa]$ of pressure (10 points), at T=300K, use the function gibbs_serie_p:

Below, the equilibrium curve in the P/T space, for the reaction enstatite + corundum <--> pyrope is computed.

The Gibbs free energy of each mineral phase involved is evaluated by integrating

$$dG = V{\rm d}P - S{\rm d}T$$

The contribution to $G$ from the integration over $T$ is computed at the reference pressure of 1 bar, by knowing the specific heat $C_P$ as a function of $T$ (at $P=1$ bar). The contribution to $G$ from the integration over $P$ is computed by knowing $V$ as a function of $P$ at the fixed temperature $T$; $V(T,P)$ is derived from the knowledge of the thermal expansion (as a function of $T$) and the equation of state, for which

$$K=K_0 + Kp\cdot P + b\cdot(T-T_{ref})$$

where $K_0$ is the bulk modulus at $P=1$ bar and $T=T_{ref}=298.15$ K, and $b$ is a constant.

Values for all these quantities are stored in the database defined within the program (Holland and Powell database, and printed by the methods py.info, ens.info and cor.info).

By issuing the command equilib, specifying a range of temperatures (300, 500K), and a number of $T$ points in the range (8), the equilibrium pressure (in GPa) of the reaction, at each temperature is computed:

Note the syntax of the command equilib: two lists must be provided (enclosed in square brackets) and, precisely, the prod list of the products of the reaction, and the rea list of the reactants; the mineral species in each list are identified by their names, as stored in the internal database of the program, and are followed by the corresponding stoichiometric coefficients they have in the reaction. In this example, we have 3/2 ens + cor <--> py

Quantum-mechanical based properties for bulk modulus (and its pressure and temperature derivatives), thermal expansion, specific heat and entropy can be uploaded. The upload_mineral command can be used to make appropriate fits of all those quantities in a given temperature range (300-500K in the example), that are subsequently uploaded in the py object.

Here, the volume range for the frequency fitting is expanded to cover an higher temperature range at low pressure; the conditions under which the calculation is done are checked by the method info.show(), and then the upload_mineral command is issued:

Note the mqm keyword appearing as argument of the upload_mineral command: it specifies that ab initio thermodynamics data are referred to pyrope (data for enstatite and corundum are from the Holland & Powell database).

The new data for pyrope are:

The new equilibrium curve for the reaction, in the P/T space, is then:

which is in very good agreement with that derived by using the Holland & Powell data for pyrope.