Static equation of state

By using the CRYSTAL17 program, the computation of equation of state (EoS) can be easily performed by means of the EOS keyword in the geometry section of the input; the keyword is followed by the specification of the volume range over which the computation of static energy is to be performed. The program proceeds by optimizing the geometry (constant volume optimization) at a number of volumes within the range, gathering the resulting energies and by performing a volume-integrated EoS fit to the set of volume-energy values.

While CRYSTAL17 keeps the classification of the integrals fixed over all the volume range, usually there is a lack of precision concerning the calculation of the energy at volumes that are significantly different from the chosen reference volume (where the integral classification is performed). The issue is likely due to the variation of the density and orientation of the grid used for the numerical evaluation of the DFT functionals. The problem can be seen by looking at the numerical integration of the electron density that should return a number of electrons in the unit cell as close as possible to the actual value. Expecially if the explored volume range is large, the problem can lead to a failure in the estimation of the equation of state.

A way out to the problem is the computation of the EoS by using the pressures, at each volume, as computed by CRYSTAL17 through the analytic derivatives of the energy with respect to the lattice parameters: the computed stress tensor is used to get the value of the hydrostatic pressure. The procedure is quite accurate provided that strict tolerances are imposed for the optimization of the geometry (the same tolerances used to get a geometry suitable for a frequency calculation should be set). The CRYSTAL17 keyword EXTPRESS can be used to optimize a structure under an imposed external pressure.

An example is shown in the following. The input contains a reference to a static energy file (under the keyword STATIC) that must always be present, together with the reference to a static pressure file (under the keyword PSTATIC); the file name (p_static.dat in the example) is followed by the value of the energy (E0, in hartree) at the static equilibrium unit cell volume. While E0 is unuseful in the computation of the EoS, it is required to get the value of the static energy at any different volume, by using a volume-integrated BM3-EoS whose parameters (but E0) are estimated by fitting a BM3-EoS on the pressure-volume data.

The EoS parameters shown above are from the V-BM3 fit of the energy-volume data (contained in the static.dat file specified under the keyword STATIC). To compute the EoS by using the P/V data instead, the function p_static is to be invoked:

A difference of about 1 GPa between the K0's computed by means of the two methods can be seen.

In the second figure above a plot of the energies is reported: the dashed line is from the P/V BM3 fit, whereas the continuos line refers to the the original E/V V-BM3 fit; stars mark the actual energy values found in the static.dat file.

Note that, for values of the volume larger than about 391 A^3 (the star marker that is located furthest to the right), the energy from the STATIC case (continuous line) is extrapolated from the V-BM3 fit done on the actual points at lower volume. On the contrary, the pstatic energy curve in not extrapolated at any volume in the whole range (but it reflects the actual P/V values found in the p_static.dat file).

The discrepancy between the two curves is small: about 390 J/mol on average (over a total energy of about -7.5 10^9 J/mol) and increases with the volume, expecially in the region where the V-BM3 energies are extrapolated (as apparent from the latter figure above).

The reason for which the p_static procedure has been implemented is evident in the following run, where 2 further E/V point have been added (v_add and e_add optional arguments to the p_static function):

The 2 red markers in the second figure above are the added V/E values, that are not used (and cannot be used) for the V-BM3 fit. Such values are from the the very same EOS calculation performed by CRYSTAL17 that was originally extended up to a volume close to 406 A^3. For some reason (likely related to the DFT grid problem quoted above), those energies are totally out of scale with respect to those computed at lower volumes. Indeed the bulk modulus computed by CRYSTAL17 at the end of the EOS calculation was negative!. Apart from the EOS procedure, the computation of the energies at those high volumes led to the same (wrong) values even in FIXINDEX/GEOM calculations by keeping the integral classification fixed at the reference low volume (358 A^3).

Wanting to get static energy curves reliable enough for the calculation of thermodynamic properties ranging from high pressures (low volumes) and high temperatures (large volumes), the p_static procedure here shown is advisable.