Input and Output

Input Description

Gibbs2 reads its input from a single file, typically with extension .ing. This file contains keywords parsed by the program that indicate the details of the system under study and the tasks to run. The input file is free-format and case-insensitive. Comments are preceded by a # symbol and they can be used anywhere in the input. Long lines can be continued using the backslash (\) character.

A simple input file corresponding to a pseudopotential/plane-waves LDA calculation in MgO is:

MM 40.3044
NAT 2
PHASE mgo
8883583665837   -73.5171659350000
0358791612784   -73.5360133400000
1833999559730   -73.5508544000000
3309207506677   -73.5624133500000
4784415453624   -73.5712754950000
6259623400570   -73.5779108750000
7734831347517   -73.5826963000000
9210039294464   -73.5859397750000
0685247241411   -73.5878968050000
2160455188357   -73.5887770750000
3635663135304   -73.5887555350000
5110871082251   -73.5879775400000
6586079029198   -73.5865593900000
8061286976144   -73.5846083150000
9536494923091   -73.5822069800000
1011702870038   -73.5794278200000
2486910816984   -73.5763274850000
3962118763931   -73.5729658900000
5437326710878   -73.5693791200000
0000000000000   -73.5662512150000
ENDPHASE

The MM keyword indicates the molar mass of the primitive cell in atomic mass units (amu). In the case of MgO, Mg weighs 24.3050 amu and O weighs 15.9994, and there is one atom of each in the primitive cell so Mg+O is 40.3044. NAT is the total number of atoms in the primitive cell.

The last keyword, PHASE is the most important part of a gibbs2 input. Every time a PHASE keyword appears, a new phase is defined in the input. The label immediately following PHASE is the name that phase will have in gibbs2. Additional options for the PHASE can be input following the label, although in this case we opted for simplicity.

The PHASE keyword can either read its data (mostly the $E(V)$ curve but also other information) directly from the input file, as above, or from an external file using the FILE option:

MM 40.3044
NAT 2
PHASE mgo FILE mgo.dat

where mgo.dat contains the numerical rows from the example above. If the energy-volume data is given in the main input file, then the data for a given phase must end with the ENDPHASE keyword. Each row between PHASE and ENDPHASE corresponds to a data point. In this case, the data is simply the unit cell volume (in bohr^3, first column) and the energy per unit cell (in Hartree, second column).

Gibbs2 can work with one or more phases, but they all must correspond to the same system (i.e. the same composition). If more than one phase is given, then gibbs2 will calculate the relative stability of each phase as a function of temperature and pressure. For instance, if we were interested in the B1 -> B2 transition in MgO, the input would be:

PHASE mgo:B1
 ...
ENDPHASE
PHASE mgo:B2
 ...
ENDPHASE

However, phases need not correspond to different atomic arrangements of a solid. For instance, if we were interested in how two different thermal models apply to the same set of data, we could do:

PHASE mgo:debye TMODEL debye
 ...
endphase
PHASE mgo:debeins TMODEL debye_einstein
 ...
endphase

where the difference between the two “phases” is that they use different thermal models (Debye and Debye-Einstein). In the output, each phase is identified by an integer, corresponding to the order in which that phase appears in the input. In the input above, mgo:debye corresponds to phase 1 and mgo:debeins is phase 2.

Going back to our simple example given by the input above, the input can be saved to file (e.g. mgo.ing) and then run using:

gibbs2 mgo.ing mgo.outg

The output of the calculation is written to mgo.outg. In addition, gibbs2 generates a number of files with the same root as the input file (mgo):

mgo_all_p.gnu
mgo_all_t.gnu
mgo.efit
mgo_efit.aux
mgo_efit.gnu
mgo.eos
mgo.eos_static

We now examine each of these files in turn.

Output description

The standard output file generated by gibbs2 (mgo.outg in the example above) contains a lot of useful information and a list of the auxiliary files generated. It is structured in blocks denoted by a header that starts with an asterisk. The first block (Input) contains information about the system and the input variables. It is mostly self-explanatory:

* Input
  Title:
  Output file (lu= 2): stdout
  Units: output is in atomic units, except where noted.
  Number of atoms per formula:   2
  Molecular mass (amu):       40.30440000
  ...

Next, the output shows the pressure range on which the equation will be calculated. The maximum pressure is determined by the slope at the first point of the static energy-volume curve.

* Pressure range examined
  Min_phases{p_max} (GPa):      118.530
  Pressure range (GPa):        0.000 ->      500.000
  Number of p points:    100

After that, information about each phase in the input is given, including their static equilibrium properties, statistical measures of the quality of the fit, and the static equation of state. The latter is similar to the contents of the .eos_static file (see below):

* Phase information after initial setup
+ Phase  1 (mgo)
  Number of volume points:   20
  p(V) input data? F
  Pressure range (GPa):      -27.120 ->      118.530
  Number of interpolated fields :    0
  Input units:
    Volume : bohr^3
  ...
  Static equilibrium volume (bohr^3):       121.1510427522
  Static equilibrium energy (Ha):       -73.5888704284
  Static equilibrium energy (kJ/mol):   -193207.5498742671
  Static bulk modulus (GPa):      171.876312
  ...
  Temperature model: Debye, Td from static B(V).
  All data points are ACTIVE for dynamic calculations

  Fit to static E(V) data:
# Copy in file : mgo.eos_static
#    p(GPa)          E (Ha)     V(bohr^3)      V/V0   p_fit(GPa)   B(GPa)         Bp      Bpp(GPa-1)
      0.0000  -7.358887030E+01   121.1510  1.0000000     0.0000   171.8631    4.1045942  -2.3291E-02
     ...
# Polynomial fit to strain: 
# Degree :  9
#      p(x) = a_0 + a_1 * f(V) + ... + a_n * f(V)^n
# a_00 =  -7.358877764647E+01
   ...
# a_09 =  -4.588683891761E+03
# V_scal (bohr^3) =   1.192160455188E+02
# p_scal (GPa) =       0.000000000000

The static energy-volume plots are calculated next, and then the Debye temperatures are computed, which are required to apply the Debye model in the Slater formulation. This is the default temperature model, where the Debye temperatures are calculated from the static bulk moduli $B(V)$ and the Poisson ratio is assumed to be volume-independent an equal to 1/4. The computed Debye temperatures are given in the output:

* Computed Debye temperatures from static data
+ Phase mgo
# ThetaD at static eq. volume:     836.48
# V(bohr^3)  Tdebye(K)   Tdebye_slater(K)
    81.8884    1563.12    1563.12
...
   160.0000     435.80     435.80

The smallest Debye temperature is usually found at the most expanded volume. This minimal temperature is used to set the default temperature range for the calculation, which is given next in the output:

* Temperature range examined
  Min_{DebyeT} (K):      435.799
  Temperature range (K):        0.000 ->      653.699
  Number of T points:    100

Once the static calculation is over, gibbs2 loops over temperatures and pressures and calculates the thermodynamic properties for at each (p,T) pair. The results are written to the eos file:

* Calculated temperature effects
  Writing file : mgo.eos

and some plotting scripts are generated to ease the interpretation of the .eos file.

Static energy plot (efit file)

The files mgo.efit, mgo_efit.aux, and mgo_efit.gnu can be used to create a plot of the static $E(V)$ data together with the fitted equation of state. To generate the plot, do:

gnuplot mgo_efit.gnu

For this to work you need the gnuplot program as well epstopdf and pdfcrop. (If the last two are not available, gnuplot will still generate an .eps file with the plot.)

Static equation of state (eos_static file)

The mgo.eos and mgo.eos_static contain the calculated equation of state and thermodynamic properties as a function of pressure and temperature. mgo.eos_static gives the static properties (i.e. without temperature) derived from fitting to the input $E(V)$ data. Each row in this file corresponds to a different pressure. The properties given in mgo.eos_static are:

Column	Symbol	Equation	Description	Unit
1	`p`	$p$	Pressure	GPa
2	`E`	$E_{\rm sta}$	Static energy	Hartree
3	`H`	$H = E_{\rm sta} + pV$	Enthalpy	Hartree
4	`V`	$V(p)$	Volume	bohr^3
5	`V/V0`	$V(p)/V(0)$	Volume divided by the equilibrium (zero-pressure) volume	–
6	`p_fit`	$p_{\rm fit}$	Pressure calculated from the EOS fit to the data (for testing purposes)	GPa
7	`B`	$B = -V\frac{\partial p}{\partial V}$	Bulk modulus	GPa
8	`Bp`	$B' = \frac{dB}{dp}$	First pressure derivative of the bulk modulus	–
9	`Bpp`	$B'' = \frac{d^2B}{dp^2}$	Second pressure derivative of the bulk modulus	1/GPa

The list of pressures at which the properties are calculated can be modified using the PRESSURE keyword. By default, gibbs2 calculates 100 pressure points from zero up to the maximum pressure allowed by the static energies (or 500 GPa). If several phases are given in the input, properties are calculated for each phase in turn. In the .eos_static file, each phase is indicated by a # Phase line.

Thermal equation of state (eos file)

The file mgo.eos contains the thermodynamic properties as a function of temperature and pressure. Each row in this file represents a different pair of pressure and temperature values. The columns are:

Column	Symbol	Equation	Description	Unit
1	`p`	$p$	Pressure	GPa
2	`T`	$T$	Temperature	K
3	`V`	$V(T,p)$	Volume	bohr^3
4	`Estatic`	$E_{\rm sta}$	Static energy	Hartree
5	`G`	$G = E_{\rm sta} + pV + F_{\rm vib} + F_{\rm el} + ...$	Gibbs free energy	kJ/mol
6	`Gerr`	$G_{\rm err}$	Difference in Gibbs free energy between the value calculated at the equilibrium value calculated at the equilibrium volume directly or interpolated from the EOS fit, per `nat` atoms. For testing purposes only.	kJ/mol
7	`p_sta`	$p_{\rm sta} = -\frac{d E_{\rm sta}}{dV}$	Static pressure	GPa
8	`p_th`	$p_{\rm th} = -\frac{\partial (F-E_{\rm sta})}{\partial V} = p - p_{\rm sta}$	Thermal pressure	GPa
9	`B`	$B_T = -V\left(\frac{\partial p}{\partial V}\right)_T$	Isothermal bulk modulus	GPa
10	`U-Esta`	$U_{\rm vib} + U_{\rm el} + ...$	Non-static contribution to the internal energy	kJ/mol
11	`Cv`	$C_V = \left(\frac{\partial U}{\partial T}\right)_V$	Constant-volume heat capacity	J/K/mol
12	`F-Esta`	$F_{\rm vib} + F_{\rm el} + ...$	Non-static contribution to the Helmholtz free energy	kJ/mol
13	`S`	$S$	Entropy	J/K/mol
14	`ThetaD`	$\Theta_D$	Debye temperature calculated from the static data	K
15	`gamma`	$\gamma = \frac{\alpha B_T V}{C_V}$	Average Grüneisen parameter	–
16	`alpha`	$\alpha = \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_p$	Volumetric thermal expansion coefficient	$10^{-5}$ /K
17	`dp/dT`	$\left(\frac{dp}{dT}\right)_V = \alpha B_T$	Constant-volume temperature-derivative of pressure	GPa/K
18	`Bs`	$B_S = -V\left(\frac{\partial p}{\partial V}\right)_S$	Adiabatic bulk modulus	GPa
19	`Cp`	$C_p = \left(\frac{\partial H}{\partial T}\right)_p$	Constant-pressure heat capacity	J/K/mol
20	`B_Tp`	$B_T' = \left(\frac{dB_T}{dp}\right)_T$	First pressure derivative of the isothermal bulk modulus	–
21	`B_Tpp`	$B_T'' = \left(\frac{d^2B_T}{dp^2}\right)_T$	Second pressure derivative of the isothermal bulk modulus	1/GPa
22	`Fvib`	$F_{\rm vib}$	Vibrational contribution to the Helmholtz free energy	kJ/mol
23	`Fel`	$F_{\rm el}$	Electronic contribution to the Helmholtz free energy	kJ/mol
24	`Uvib`	$U_{\rm vib}$	Vibrational contribution to the internal energy	kJ/mol
25	`Uel`	$U_{\rm el}$	Electronic contribution to the internal energy	kJ/mol
26	`Svib`	$S_{\rm vib}$	Vibrational contribution to the entropy	J/K/mol
27	`Sel`	$S_{\rm el}$	Electronic contribution to the entropy	J/K/mol
26	`Cv_vib`	$C_{v,{\rm vib}}$	Vibrational contribution to the constant-volume heat capacity	J/K/mol
27	`Cv_el`	$C_{v,{\rm el}}$	Electronic contribution to the constant-volume heat capacity	J/K/mol

The specific formulas and the method with which these quantities are calculated are given in the gibbs2 article. If several phases are given in the input, properties are calculated for each phase in turn. In the .eos file, each phase is indicated by a # Phase line. The pressure and temperature points at which the thermodynamic properties are calculated in the .eos file can be changed using the PRESSURE and TEMPERATURE keywords. By default, 100 temperature points are used between 0 and 1.5 times the minimum Debye temperature.

Plots of thermodynamic properties (all_p.gnu and all_t.gnu files)

Two gnuplot scripts (in the above example mgo_all_p.gnu and mgo_all_t.gnu) can be used to make simple plots of the information contained in the .eos file. To make the plots, use:

gnuplot mgo_all_p.gnu
gnuplot mgo_all_t.gnu

This creates a good number of pdf files, each containing the evolution of a thermodynamic property from the .eos file with either pressure at constant temperature (mgo_p_xx.pdf) or temperature at constant pressure (mgo_t_xx.pdf). The xx identifiers correspond to the column in the .eos file. For example, mgo_p_03.pdf contains the V(p) isotherms, because volume is column 3 in the .eos file. Similarly, mgo_t_09.pdf contains the $B_T(T)$ curves at constant pressure.

Phase Transitions (tpstab, dgtp, ptrans files)

If more than one phase is given in the input, gibbs2 can calculate stability diagrams and transition pressures from the data for the different phases by determining which has lowest free energy at a given temperature and pressure. When comparing phases with different number of atoms in the cell it is essential that the energy and volume of each phase correspond to a number of atoms equal to NAT times their value for the Z keyword.

When more than one phase is present, gibbs2 automatically calculates the static transition pressures, which are written to the output:

* Static transition pressures (linear interpolation)
#        Pressure range (GPa)      Stable phase
        0.0000 -->     503.1471          b1
      503.1471 -->     600.0000          b2

In addition, a plot is generated for the static enthalpy vs. pressure diagram, ending in _dH.gnu (the gnuplot script) and _dH.aux (the auxiliary data for the plot). This plot presents the difference in enthalpy between all phases and the first phase in the input as a function of applied pressure. Note this plot does not contain vibrational effects.

Several additional files are written that contain the temperature- and pressure-dependence of phase stability:

A file with extension .tpstab is created that gives the stable phase (i.e. the one that has minimum Gibbs free energy) as a function of pressure and temperature. The columns are, in order: temperature, pressure, identity of the stable phase, Gibbs free energy, volume, and isothermal bulk modulus.
The file with extension .dgtp gives the difference in Gibbs free energy between all the phases and the first phase, at all the temperatures and pressures.
The evolution of the transition pressures with temperature is represented in the .ptrans file and its associated gnuplot script (_ptrans.gnu). This is essentially the pressure-temperature phase diagram.

Optional keywords

SET ROOT root.s

Sets the base name for the files generated by gibbs2 to root.s. Default: the naem fo the input file without extension.

SET NOEFIT

This command instructs gibbs2 not to generate the static energy plot (efit).

SET NOPLOTDH

Do not generate the static enthalpy vs. pressure plots.

SET NOTRANS

Do not calculate phase transition pressures.

SET ERRORBAR|ERROR_BAR|ERRORBARS|ERROR_BARS

Calculate error bars for the dynamic properties. Only applies to average polynomial EOS.

SET WRITELEVEL {0|1|2}

Sets the verbosity level. It can be 0 (no output files except stdout), 1 (only .eos and .eos_static), or 2 (all files written). The default is 2.