<jupyter_start><jupyter_text>Main features of the program and more important functions implementedThe 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 expansion1. 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 properties3. Specific heat at costant volume and at constant pressure;4. Entropy;5. Helmholtz free energy;6. Gibbs free energy;7. 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 analysis8. 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](https://pubs.geoscienceworld.org/msa/ammin/article-abstract/99/7/1449/46740/ab-initio-thermodynamic-and-thermophysical?redirectedFrom=fulltext) 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](https://qm-thermodynamics.readthedocs.io/en/main/_static/Dispersion.html). LO-TO splittingThe impact of the LO-TO splitting on properties can also be evaluated. See the tutorial [here](https://qm-thermodynamics.readthedocs.io/en/main/_static/LO_TO_splitting.html). AnharmonicityThe 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](https://qm-thermodynamics.readthedocs.io/en/main/_static/anharm.html). The required data (scans of the anharmonic potential) can be obtained by using the [CRYSTAL17](https://www.crystal.unito.it/index.php) 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](https://qm-thermodynamics.readthedocs.io/en/main/_static/basic_eos_tutorial.html).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).<jupyter_code>%run bm3_thermal_2.py<jupyter_output>This is BMx-QM program, version 2.4.0 - 18/02/2021 
Run time:  2021-05-24 09:48:46.664774

File quick_start.txt found in the master folder
Input files in 'pyrope' folder

Input file

STATIC
pyrope_static.dat
VOLUME
volume.dat
FREQ
pyrope_freq_2.dat
EXP
experimental.txt
LO
LO.txt
FITVOL
POLY
725. 769. 16 3.
FU
4 20
SET
0 1 4 5 8 10 15 21
CP
0.  2. -2. 1. -1. -0.5
ALPHA
0. -2  -1 -0.5
KIEFFER
100. 100. 100.
END

-------- End of input file -------

Frequencies corrected for LO-TO splitting

Static BM3 EoS

Bulk Modulus: 173.95 (0.07) GPa
Kp:             4.30 (0.04)
V0:           756.1932 (0.00) A^3
E0:            -1.14182286e+04 (1.56e-06) hartree


Kieffer model for the acoustic branches activated<jupyter_text>Dealing with the frequenciesIn 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 temperature2. 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:- **No fit of frequencies as function of the cell volume**: this option can be used when a relatively large set of frequency/volume data is available, and in cases for which the *numerical noise* possibly affecting the *ab initio* computation of the frequencies, as they change with the cell volume, is small; otherwise fits can be used:    - **Polynomial fits** of frequencies as function of volume or- **Spline fits** of frequencies as function of volumeIn the current example, a polynomial fit is requested:~~~~~~~~~~FITVOLPOLY725. 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:<jupyter_code>eos_temp(300, prt=False)<jupyter_output>** BM3 fit **

EoS at the temperature of 300.00 K

Bulk Modulus: 163.17   (0.00) GPa
Kp:             4.11   (0.00)
V0:           767.4477 (0.000) A^3

Fix status: False
Fit of frequencies is active
Polynomial fit: degree  3
Volume range: 725.00 - 769.00, points=16<jupyter_text>The other possible choice for the frequency fitting is the *spline* fit. This can be required in *input.txt* with:~~~~~~~~~~~~~FITVOLSPLINE725. 769. 16 3 1.~~~~~~~~~~~~~The last float number is the *smoothness* parameter. Information about the spline fit can be found [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.UnivariateSpline.htmlscipy.interpolate.UnivariateSpline); 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:* POLY fit: deg=4, npoint=16* SPLINE fit: deg=3, smooth=5., npoint=16As an example:<jupyter_code>set_spline(3, 10.)<jupyter_output>Volume range 725.5483 - 772.8196 defined for 'SPLINE' fit<jupyter_text>**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***:<jupyter_code>info.show()<jupyter_output>Current settings and results

Static data            ** min, max volumes: 725.7483, 772.6196; points: 8
Frequency volume range ** min, max volumes: 725.7483, 772.6196; points: 22
Selected freq. sets    ** min, max volumes: 725.7483, 772.6196; points: 8
Frequency sets:           [0, 1, 4, 5, 8, 10, 15, 21]

Fit of frequencies     ** type: spline, degree: 3, smooth: 10.0
                          min, max volumes: 725.5483, 772.8196; points 16

** Static EoS (BM3) **
K0: 173.95 GPa, Kp: 4.30, V0: 756.1932 A^3

** BM3 EoS from the last computation, at the temperature of 300.00 K **
K0: 163.17 GPa, Kp: 4.11, V0: 767.4477 A^3
Kp not fixed

No zone center excluded modes
No off center excluded modes

Kieffer model on; frequencies 100.00 100.00 100.00 cm^-1

Frequencies corrected for LO-TO splitting.


**** Volume driver for volume_dir function ****
Delta: 2.0; degree:  2; left: 2.0; right: 2.0, Kp_fix: True; t_max: 500.00
EoS shift: 0.0; Quad_shrink:  4; T_dump: 0.0; Dump fact.: 1.0, T_last [...]<jupyter_text>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:<jupyter_code>set_volume_range(725., 769., prt=True)<jupyter_output>Volume range 725.0000 - 769.0000 defined for 'SPLINE' fit<jupyter_text>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:<jupyter_code>eos_temp(300,prt=False)<jupyter_output>** BM3 fit **

EoS at the temperature of 300.00 K

Bulk Modulus: 163.17   (0.00) GPa
Kp:             4.11   (0.00)
V0:           767.4477 (0.000) A^3

Fix status: False
Fit of frequencies is active
Spline fit: degree  3, smooth: 10.0
Volume range: 725.00 - 769.00, points=16<jupyter_text>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:<jupyter_code>fit_off()<jupyter_output><empty_output><jupyter_text>To see the effect on the equation of state:<jupyter_code>eos_temp(300, prt=False)<jupyter_output>** BM3 fit **

EoS at the temperature of 300.00 K

Bulk Modulus: 160.00   (1.56) GPa
Kp:             4.70   (0.40)
V0:           767.6477 (0.066) A^3

Fix status: False
Fitting is off<jupyter_text>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:<jupyter_code>bulk_dir(300)<jupyter_output>Warning: frequency fit is off; use of poly or spline fits is mandatory for bulk_dir<jupyter_text>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:<jupyter_code>set_poly(3)
bulk_dir(300)<jupyter_output>BM3 EoS from P(V) fit

K0:    163.16   (0.00) GPa
Kp:      4.11   (0.00)    
V0:  767.4479   (0.00) A^3<jupyter_text>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:<jupyter_code>check_poly_list([0, 2, 36])<jupyter_output><empty_output><jupyter_text>Take a look at [this](https://qm-thermodynamics.readthedocs.io/en/main/_static/Dealing_with_the_frequencies.html) tutorial to see other commands of the same type.   Equation of stateAs already seen above, the two commands for optimizing the equation of state at a given temperature are:* ***eos_temp(tt)***: optimization of a *volume-integrated F(V)* equation* ***bulk_dir(tt)***: optimization of a *P(V)* equationthe 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):<jupyter_code>static(plot=True)<jupyter_output>Static BM3 EoS

Bulk Modulus: 173.95 (0.07) GPa
Kp:             4.30 (0.04)
V0:           756.1932 (0.00) A^3
E0:            -1.14182286e+04 (1.56e-06) hartree<jupyter_text>A 4^th order Birch-Murnaghan EoS can also be optimized; see the relevant section in [basic_eos_tutorial](https://qm-thermodynamics.readthedocs.io/en/main/_static/basic_eos_tutorial.html) 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:<jupyter_code>set_fix(4.5)
bulk_dir(300)<jupyter_output>BM3 EoS from P(V) fit

K0:    161.32   (0.14) GPa
Kp:      4.50   (0.00)    
V0:  767.5211   (0.02) A^3<jupyter_text>**Note** that some commands can override the directive **set_fix** by providing the optional keyword argument ***fix***: precisely, if ***fix=0.***, *Kp* is optimized:<jupyter_code>bulk_dir(300, fix=0.)<jupyter_output>BM3 EoS from P(V) fit

K0:    163.16   (0.00) GPa
Kp:      4.11   (0.00)    
V0:  767.4479   (0.00) A^3<jupyter_text>The command ***reset_fix*** deactivates the ***set_fix*** directive:<jupyter_code>reset_fix()
bulk_dir(300)<jupyter_output>BM3 EoS from P(V) fit

K0:    163.16   (0.00) GPa
Kp:      4.11   (0.00)    
V0:  767.4479   (0.00) A^3<jupyter_text>**Note** however that in:<jupyter_code>bulk_dir(300, fix=4.5)<jupyter_output>BM3 EoS from P(V) fit

K0:    161.32   (0.14) GPa
Kp:      4.50   (0.00)    
V0:  767.5211   (0.02) A^3<jupyter_text>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:<jupyter_code>reset_fix()
bulk_serie(300,600, 12, degree=1)
set_fix(4.5)
bulk_serie(300,600, 12, degree=1)<jupyter_output>Results from the fit (from high to low order)
[-3.71e-02  1.73e+02]<jupyter_text>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 expansionThe 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:<jupyter_code>fix_status()<jupyter_output>Fix status: True
Kp fixed at 4.50<jupyter_text>To inquire about the current fit status of the frequencies, the command ***fit_status*** must be given:<jupyter_code>fit_status()<jupyter_output>Fit of frequencies is active
Polynomial fit: degree  3
Volume range: 725.00 - 769.00, points=16<jupyter_text>The thermal expansion at *T = 300K* and *P = 0 GPa* is:<jupyter_code>thermal_exp_p(300, 0)<jupyter_output>Thermal expansion: 3.14e-05 K^-1
Bulk modulus:        161.36 GPa
Pressure:              0.00 GPa
Volume:              767.5342 A^3<jupyter_text>In this calculation *Kp* (which is an intermediate result of the algorithm) can be optimized by using the optional keyword ***fix***:<jupyter_code>thermal_exp_p(300, 0, fix=0.)<jupyter_output>Thermal expansion: 3.10e-05 K^-1
Bulk modulus:        163.17 GPa
Pressure:              0.00 GPa
Volume:              767.4477 A^3<jupyter_text>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:  ~~~~~~~~~~~~~~ALPHA0. -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.<jupyter_code>alpha_serie(20, 600, 18, 0, prt=False, fit=True)<jupyter_output>Warning: volume out of range; reduce temperature<jupyter_text>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 heatThe 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*:<jupyter_code>entropy_p(300, 0)<jupyter_output>Entropy:  279.52 J/mol K
Specific heat (at constant volume):  329.60 J/mol K
Pressure:  0.00 GPa; Volume 767.5342 A^3<jupyter_text>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$:<jupyter_code>cp(300,0, prt=True)<jupyter_output>Cp: 335.10, Cv: 329.60, Cp/Cv: 1.01669, alpha: 3.136e-05, K: 161.36<jupyter_text>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):<jupyter_code>cp_serie(40,400,12, 0, prt=True)<jupyter_output>Temp.     Cp
 40.00  13.47
 72.73  55.65
105.45 108.29
138.18 159.84
170.91 205.99
203.64 246.68
236.36 281.45
269.09 310.64
301.82 335.47
334.55 357.36
367.27 376.24
400.00 394.73<jupyter_text>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):<jupyter_code>cp_serie(40,400,18, 0, prt=False, graph=False)<jupyter_output><empty_output><jupyter_text>Such values of the coefficient can be saved in a variable like:<jupyter_code>cp_coeff=cp_serie(40,400,18, 0, prt=False, graph=False)<jupyter_output><empty_output><jupyter_text>and used to compute $C_P$ at any temperature by using the polynomial fitting function ***cp_fun***:<jupyter_code>print(cp_coeff)
round(cp_fun(300, *cp_coeff),2)<jupyter_output>[ 1.17314930e+03 -1.60123258e-04 -6.97102652e+05 -4.07579977e-02
  9.32201845e+04 -1.93065830e+04]<jupyter_text>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 splittingIn 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):<jupyter_code>lo.off()<jupyter_output>LO-TO splitting not taken into account<jupyter_text>To see the effect of the correction on the EoS parameters in this case:<jupyter_code>reset_fix()
eos_temp(300,prt=False)<jupyter_output>** BM3 fit **

EoS at the temperature of 300.00 K

Bulk Modulus: 163.02   (0.00) GPa
Kp:             4.11   (0.00)
V0:           767.5001 (0.000) A^3

Fix status: False
Fit of frequencies is active
Polynomial fit: degree  3
Volume range: 725.00 - 769.00, points=16<jupyter_text>which is to be compared with:<jupyter_code>lo.on()
eos_temp(300, prt=False)<jupyter_output>Frequencies corrected for LO-TO splitting

 ** BM3 fit **

EoS at the temperature of 300.00 K

Bulk Modulus: 163.17   (0.00) GPa
Kp:             4.11   (0.00)
V0:           767.4477 (0.000) A^3

Fix status: False
Fit of frequencies is active
Polynomial fit: degree  3
Volume range: 725.00 - 769.00, points=16<jupyter_text>The bulk modulus is not significantly affected by such LO-TO correction; however thermodynamic properties can be changed:<jupyter_code>lo.off()
cp(300,0,prt=True)<jupyter_output>LO-TO splitting not taken into account

Cp: 335.85, Cv: 330.30, Cp/Cv: 1.01681, alpha: 3.135e-05, K: 163.02<jupyter_text>to be compared with:<jupyter_code>lo.on()
cp(300,0,prt=True)<jupyter_output>Frequencies corrected for LO-TO splitting

Cp: 335.09, Cv: 329.64, Cp/Cv: 1.01654, alpha: 3.104e-05, K: 163.17<jupyter_text>The Kieffer's model for the acoustic branchesThe 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:<jupyter_code>kieffer.on()
kieffer.freq(300., 350., 400.)
kieffer.plot()<jupyter_output>Kieffer correction on<jupyter_text>At any temperature, the contribution to $F$ from the acoustic branches can be recovered with the method ***kieffer.get_value***; for instance:<jupyter_code>kieffer.get_value(300)<jupyter_output><empty_output><jupyter_text>**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.:<jupyter_code>eos_temp(300, prt=False)<jupyter_output>** BM3 fit **

EoS at the temperature of 300.00 K

Bulk Modulus: 163.17   (0.00) GPa
Kp:             4.11   (0.00)
V0:           767.4477 (0.000) A^3

Fix status: False
Fit of frequencies is active
Polynomial fit: degree  3
Volume range: 725.00 - 769.00, points=16<jupyter_text>On the other hand, other properties can be significantly affected:<jupyter_code>cp(300,0,prt=True)<jupyter_output>Cp: 334.13, Cv: 328.68, Cp/Cv: 1.01658, alpha: 3.104e-05, K: 163.17<jupyter_text>to be compared with:<jupyter_code>kieffer.off()
cp(300,0,prt=True)<jupyter_output>Kieffer correction off

Cp: 328.94, Cv: 323.49, Cp/Cv: 1.01685, alpha: 3.104e-05, K: 163.17<jupyter_text>Thermodynamics and mineral equilibriaFunctions exist to compute and to plot the Gibbs free energy $G=F+PV$ as function of $T$ or $P$; for instance:<jupyter_code>set_fix(4.39)
gibbs_serie_t(200,400,10,0., v0=py.v0, g0=py.g0)<jupyter_output><empty_output><jupyter_text>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):<jupyter_code>py.info()<jupyter_output>Mineral: pyrope

K0: 173.70 GPa, Kp: 4.00, dK0/dT: -0.0261 GPa/K, V0: 11.3180 J/bar
G0: -5934105.00 J/mol, S0: 266.30 J/mol K

EoS type: bm 
Cp coefficients and powers:
+6.3350e+02  +0.0
-5.1961e+06  -2.0
-4.3152e+03  -0.5

Alpha coefficients and powers:
+4.3600e-05  +0.0
-4.3600e-04  -0.5<jupyter_text>To compute the $G$ energy in the range $[0, 10GPa]$ of pressure (10 points), at T=300K, use the function ***gibbs_serie_p***:<jupyter_code>gibbs_serie_p(0, 10, 10, 300., v0=py.v0, g0=py.g0)<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>equilib(300,600,12, prod=['py',1], rea=['ens', 1.5, 'cor',1])<jupyter_output>T (K)  P (GPa)  DH(J/mol)  DS (J/mol K)  DV (J/bar)  Slope (bar/K)
300.0    3.67    5306.981      17.690       -0.527       -33.55    
327.3    3.58    5805.002      17.738       -0.530       -33.46    
354.5    3.49    6313.443      17.807       -0.533       -33.40    
381.8    3.40    6823.242      17.870       -0.536       -33.32    
409.1    3.31    7327.933      17.913       -0.539       -33.21    
436.4    3.22    7822.834      17.927       -0.543       -33.04    
463.6    3.13    8304.525      17.912       -0.546       -32.81    
490.9    3.04    8770.490      17.866       -0.549       -32.53    
518.2    2.95    9218.875      17.791       -0.552       -32.20    
545.5    2.86    9648.317      17.689       -0.556       -31.82    
572.7    2.78   10057.827      17.561       -0.559       -31.40    
600.0    2.69   10446.696      17.411       -0.563       -30.95<jupyter_text>**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:<jupyter_code>reset_fix()
set_volume_range(725,772)
info.show()
upload_mineral(200,600, mqm='py', volc=True)<jupyter_output>Current settings and results

Static data            ** min, max volumes: 725.7483, 772.6196; points: 8
Frequency volume range ** min, max volumes: 725.7483, 772.6196; points: 22
Selected freq. sets    ** min, max volumes: 725.7483, 772.6196; points: 8
Frequency sets:           [0, 1, 4, 5, 8, 10, 15, 21]

Fit of frequencies     ** type: poly, degree: 3
                          min, max volumes: 725.0000, 772.0000; points 16

** Static EoS (BM3) **
K0: 173.95 GPa, Kp: 4.30, V0: 756.1932 A^3

** BM3 EoS from the last computation, at the temperature of 298.15 K **
K0: 163.22 GPa, Kp: 4.11, V0: 767.4041 A^3
Kp not fixed

No zone center excluded modes
No off center excluded modes

Kieffer model off

Frequencies corrected for LO-TO splitting.


**** Volume driver for volume_dir function ****
Delta: 2.0; degree:  2; left: 2.0; right: 2.0, Kp_fix: True; t_max: 500.00
EoS shift: 0.0; Quad_shrink:  4; T_dump: 0.0; Dump fact.: 1.0, T_last  0.0

 ** BM3 fit **

EoS at the temperature of 298.15 [...]<jupyter_text>**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:<jupyter_code>py.info()<jupyter_output>Mineral: pyrope

K0: 163.22 GPa, Kp: 4.11, dK0/dT: -0.0348 GPa/K, V0: 11.3180 J/bar
G0: -5934105.00 J/mol, S0: 265.65 J/mol K

EoS type: bm 
Cp coefficients and powers:
+1.8238e+03  +0.0
+2.9213e-04  +2.0
-1.0706e+06  -2.0
-7.6691e-01  +1.0
+1.6075e+05  -1.0
-3.1437e+04  -0.5

Alpha coefficients and powers:
+7.1343e-05  +0.0
-1.5895e-01  -2.0
+5.5067e-03  -1.0
-9.8576e-04  -0.5<jupyter_text>The new equilibrium curve for the reaction, in the *P/T* space, is then:<jupyter_code>equilib(300,600,12, prod=['py',1], rea=['ens', 1.5, 'cor',1])<jupyter_output>T (K)  P (GPa)  DH(J/mol)  DS (J/mol K)  DV (J/bar)  Slope (bar/K)
300.0    3.62    3673.161      12.244       -0.542       -22.60    
327.3    3.56    4081.270      12.471       -0.540       -23.08    
354.5    3.50    4486.566      12.654       -0.539       -23.48    
381.8    3.43    4902.077      12.839       -0.538       -23.88    
409.1    3.37    5334.273      13.039       -0.536       -24.31    
436.4    3.30    5785.629      13.259       -0.535       -24.78    
463.6    3.23    6256.265      13.494       -0.534       -25.27    
490.9    3.16    6745.010      13.740       -0.533       -25.79    
518.2    3.09    7250.095      13.991       -0.532       -26.33    
545.5    3.02    7769.615      14.244       -0.530       -26.86    
572.7    2.95    8301.832      14.495       -0.529       -27.39    
600.0    2.87    8845.366      14.742       -0.528       -27.92
{'fraction_non_alphanumeric': 0.10682718200538571, 'fraction_numerical': 0.0817044194519246, 'mean_word_length': 3.995410666244659, 'pattern_counts': {'":': 0, '<': 143, '<?xml version=': 0, '>': 143, 'https://': 9, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '45', 'dir': 'jupyter-structured-clean-dedup', 'id': '146246', 'license': 'MIT', 'n_tokens_mistral': 12029, 'n_tokens_neox': 10356, 'n_words': 4039, 'path': 'Main_features.ipynb', 'repo_name': 'Mauro-Prencipe/QM-thermodynamics', 'size': '31565'}