Basic operations to optimize a Birch-Murnaghan EoS

from static energies and vibrational frequencies, at different cell volumes, computed at the ab initio level

This tutorial shows the basic usage of the BMx system to optimize the EoS of a mineral (pyrope, in the example) starting from a set of values of static energy, at different volumes ($V$), and frequencies of vibrational normal modes also computed at different values of the unit cell volume of the crystal.

The calculation is done within the framework of the Quasi-Harmonic Approximation (QHA).

A bit of theory

The Helmholtz free energy $F$, as a function of $V$ and $T$ (temperature), is computed following its definition in statistical termodynamics:

$$F=U-TS = -k_BT \log Z$$

where $U$ is the internal energy, $S$ the entropy, $Z$ the partition function and $k_B$ the Boltzmann's constant. In turn, $Z$ is defined as:

$$Z = \sum_a e^{-\epsilon_a/k_BT}$$

where $\epsilon_a$ is the energy of the $^{th}a$ state of the crystal and the sum is extended to all the possible states. In particular, at a fixed volume, the contribution of the $^{th}i$ normal mode with frequency $\nu_i$ to the vibrational Helmholtz free energy ($F^{vib}_{i}$) is given by

$$F^{vib}_{i} = -k_BT\sum_n e^{-(n+1/2)h\nu_i/k_BT}$$

where $n$ is the vibrational quantum number of the oscillator. The total vibrational free energy is:

$$F^{vib}=\sum_i F^{vib}_{i}$$

where the sum is extended to all the normal modes of the crystal. Note that:

$$F = F^{static} + F^{vib} = U^{static} + U^{vib} - TS^{vib}$$

where $U^{static}$ is the so named static energy of the crystal ($U^{static}\equiv F^{static}$) that depends on the volume of the cell only.

Pressure ($P$) at a given volume an temperature is evaluated from its definition in term of free energy:

$$P=-\left(\frac{\partial F}{\partial V}\right)_T = -\left(\frac{\partial U^{static}}{\partial V}\right)_T -\left(\frac{\partial F^{vib}}{\partial V}\right)_T = P^{static}(V) + P^{vib}(V,T)$$

The vibrational pressure $P^{vib}$ is here defined as the partial derivative of $F^{vib}$ with respect to $V$ (at $T$ constant). In turn, $F^{vib}$ depends upon $V$ through the variation of the $\nu_i$'s by $V$; in a perfectly harmonic crystal, the vibrational frequencies are by definition independent of volume, so that the vibrational pressure is exactly zero and the total pressure is given by the static contribution only. This is the harmonic model.

The harmonic model is usually not valid as the frequencies normally do depend of volume. The Quasi harmonic approximation (QHA) is thus generally used. In fact, to compute the vibrational pressure, the frequencies variation with the unit cell volume must be estimated, so that $F$ and its derivatives with $V$ can be computed at whatever temperature.

Static energies and frequencies (at different volumes) must be calculated with some quantum-mechanical program; the one that was used here, in the case of pyrope, is CRYSTAL17.

A possible strategy to get the parameters of a chosen equation of state (EoS) is to calculate the pressure of the crystal (first derivative of $F$ with respect to $V$) in a range of volumes, at a fixed temperature. The EoS is then fitted to the pressure/volume set thus obtained.

The strategy followed in the current project if to directly fit the volume-integrated EoS to the $F$ free energy/volume set. In fact, by defining with $P_{EoS}(V)$ the $P(V)$ EoS function (for instance, a Birch-Murnaghan equation), its volume integrated form follows by:

$$P_{EoS}(V)=-\left(\frac{\partial F_{EoS}}{\partial V}\right)_T \ \rightarrow \ {\rm d}F_{EoS}(V)=-P_{EoS}(V){\rm d}V \ \rightarrow \ F_{EoS}(V) = -\int_{V_0}^{V} P_{EoS}(V){\rm d}V + F_{EoS}(V_0)$$

where $F_{EoS}(V_0)$ is the free energy at the equilibrium volume (at the fixed temperature).

Practice

  1. Prepare a folder containing;
    • a file of primitive unit cell volumes (in $\mathrm A^3$) and the corresponding energy values ($U^{static}$ in a.u.; see the theoretical section above). In the example provided, such file is 'pyrope_static.dat': the first column lists volumes; the second one lists energies;
    • a file of volumes at which frequencies of vibrational modes are computed. In the example provided, such files is 'volume.dat';
    • a file of set of frequencies computed at each volume listed in the volume file. Such file is named 'pyrope_freq.dat' in the example. Each column of the file lists the frequencies of all the modes computed at a given volume. Note that the first column lists the degeneracy of the modes. Frequency values are in $\rm{cm}^{-1}$;
    • an input.txt file with instructions to run the program. Have a look at the input file provided in the example, with the comments it contains about the meaning of each keyword.
  1. For a quick start, in the master folder (the one with the program code) do provide a file quick_start.txt just containing the name of the data folder (pyrope in the example).

Start the calculation by running the bm3_thermal_2 Python code:

The parameters of a BM3 EoS fitted to the static data (no vibrational contributions considered) are printed.

Next, issue the command eos_temp by giving a temperature (n K) as argument. This provides the BM3 EoS parameters at the specific temperature.

Besides other points, the output provides information concerning the method followed to deal with the frequencies change with the volume: in the present case, spline fits of $\nu_i(V)$ are performed for each mode $i$, as specified by the SPLINE keyword under the FITVOL directive in the input.txt file.

To deactivate the fitting of the frequencies, so to compute the EoS by using just the set of frequencies provided in input, at the corresponding volumes, use the command fit_off

Note the standard deviations on each parameter, as determined by the fit, which are higher than those obtained before, with the spline fit of the frequencies.

Beside spline fittings of the frequencies, polynomial fits can be used by issuing the command set_poly which accepts the degree of the polynomial as argument. Look at the example below, where the effect of the degree of the polynomial on the computed EoS parameters is evaluated:

The volume range for the fitting of frequencies and EoS can be changed by the function set_volume_range. The function takes, as arguments, the minimum volume (v_min), the maximum volume (v_max) and the number of points in the range:

Note that the set_volume_range function has no effect with fit_off. To change the volume range if the frequencies are not fitted, just use the SET keyword in input.txt.

EoS optimization can be performed by keeping Kp fixed. Use of the set_fix function can be used for the purpose:

The EoS can also be evaluated by fitting the BM3 equation to the pressures computed for a set of volumes (at a given temperature). In turn, such pressures are calculated as the volume derivative of the $F(V,T)$ Helmholtz free energy:

$$P=-\left(\frac{\partial F}{\partial V}\right)_T$$

At each volume, the static energy contribution to $F$ comes from a V-BM3 fit to the static energy values $[U^{static}(V)]$, whereas all the vibrational contributions comes from the vibrational partition function.

Setup the calculation by specifying a volume range for the frequency fit (and optionally check the EoS evaluated by the methods outlined above):

Then compute the EoS by fitting a BM3 function to the $P(V)$ data; the command to use is bulk_dir (which wants the temperature as argument):

4^ order Birch-Murnaghan EoS

To fit the F(V) curve by a volume integrated 4^ order Birch-Murnaghan EoS, just issue the commands start_bm4() and eos_temp:

Let's see the effect of a change in the volume range on a BM4 fit:

As in the case of BM3, also the BM4 EoS can be computed fitting P(V) instead of F(V) data. The function bm4_dir is to be used for the purpose:

Some information concerning volume ranges and fitting conditions can be gathered through the method info.show():

Direct calculation of the bulk modulus

Bulk modulus can also be computed from its definition

$$K=-V\left(\frac{\partial P}{\partial V}\right)_T$$

This is done by the function bulk_modulus_p in two different ways:

  1. a V-EoS is fitted to the $F(V)$ data at a fixed temperature; pressures at any given volume are computed from the fitted EoS and $K$ is computed as the derivative of P with respect to V;
  1. pressures are computed as derivatives of the F function, at a fixed temperature, and $K$ is obtained as the derivative of P with respect to V.

At variance with the first case, no reference to any EoS function is present in the second route.

Reload the input and start the calculation:

The route (1) explained above is followed by specifying the keyword noeos=False as argument of the function bulk_modulus_p (that is the default). The function requires the temperature (300K in the example) and the pressure (0 GPa in the example) as mandatory arguments:

The second route is the choice when the keyword noeos=True:

Note on numerical derivatives with respect to V

As said above, in a function like bulk_modulus_p, pressures are evaluated through numerical derivatives of F with respect to V (if noeos=True) and, likewise, the same type of numerical derivatives are employed to compute K. The computation of numerical derivatives involves a volume range at which the given $f(V)$ is computed, then it is fit by some polynomials of $V$ and, finally, that polynomials is analytically derived. The file parame.py contains some convenient parameters to perform such a task. For instance, the power of the polynomials is set to 3 (degree_v value, imported in the program at runtime and stored in the pr.degree_v variable). The V-range (pr.delta_v variable) is also imported as a default value but it is usually recomputed by the program as a function of the equilibrium static volume (v0) of the crystal being investigated. To do that, the parameter pr.v_frac is used (default value: 0.0015); precisely: delta=v0\pr.v_frac*.

The class volume_delta_class (whose instance is vd) is provided to redefine delta, as well as v_frac and v0, or to set delta to the default pr.delta_v value. As an example:

To use the default V-range, simply invoke the vd.off method:

In this case, a significant difference (about 2GPa) between the two computed $K$ values is observed.

Bulk modulus as a function of temperature

To compute the bulk modulus in a range of temperatures, the function bulk_modulus_p_serie can be used; its many arguments can be seen by requesting the help:

Interface to EoSFit

The function eosfit_dir can be used to write a PVT file to be used as input to the EoSFit. Such function computes the pressure, at each volume and pressure, as the derivative of the $F$ function with respect to $V$ (at constant $T$).

The function eosdir can also be used as an interface to EoSFit: it computes the pressure by using a BM3 EoS that was determined by fitting a volume-integrated BM3 EoS to $F(V)$ points.

The set of temperatures at which P(V) data must be computed is specified in input.txt, under the kewyword TEMP.