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:
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.
The impact of the LO-TO splitting on properties can also be evaluated. See the tutorial here.
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).
%run bm3_thermal_2.py
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
In the calculation of the equation of state it is possible to choose among several options; the first two choices are:
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:
eos_temp(300, prt=False)
** 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
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:
set_spline(3, 10.)
Volume range 725.5483 - 772.8196 defined for 'SPLINE' fit
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:
info.show()
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 0.0
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:
set_volume_range(725., 769., prt=True)
Volume range 725.0000 - 769.0000 defined for '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:
eos_temp(300,prt=False)
** 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
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:
fit_off()
To see the effect on the equation of state:
eos_temp(300, prt=False)
** 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
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:
bulk_dir(300)
Warning: frequency fit is off; use of poly or spline fits is mandatory for bulk_dir
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:
set_poly(3)
bulk_dir(300)
BM3 EoS from P(V) fit K0: 163.16 (0.00) GPa Kp: 4.11 (0.00) V0: 767.4479 (0.00) A^3
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:
check_poly_list([0, 2, 36])
Take a look at this tutorial to see other commands of the same type.
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):
static(plot=True)
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
A 4^th order Birch-Murnaghan EoS can also be optimized; see the relevant section in basic_eos_tutorial for details.
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:
set_fix(4.5)
bulk_dir(300)
BM3 EoS from P(V) fit K0: 161.32 (0.14) GPa Kp: 4.50 (0.00) V0: 767.5211 (0.02) A^3
Note that some commands can override the directive set_fix by providing the optional keyword argument fix: precisely, if fix=0., Kp is optimized:
bulk_dir(300, fix=0.)
BM3 EoS from P(V) fit K0: 163.16 (0.00) GPa Kp: 4.11 (0.00) V0: 767.4479 (0.00) A^3
The command reset_fix deactivates the set_fix directive:
reset_fix()
bulk_dir(300)
BM3 EoS from P(V) fit K0: 163.16 (0.00) GPa Kp: 4.11 (0.00) V0: 767.4479 (0.00) A^3
Note however that in:
bulk_dir(300, fix=4.5)
BM3 EoS from P(V) fit K0: 161.32 (0.14) GPa Kp: 4.50 (0.00) V0: 767.5211 (0.02) A^3
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:
reset_fix()
bulk_serie(300,600, 12, degree=1)
Results from the fit (from high to low order) [-3.48e-02 1.74e+02]
set_fix(4.5)
bulk_serie(300,600, 12, degree=1)
Results from the fit (from high to low order) [-3.71e-02 1.73e+02]
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.
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:
fix_status()
Fix status: True Kp fixed at 4.50
To inquire about the current fit status of the frequencies, the command fit_status must be given:
fit_status()
Fit of frequencies is active Polynomial fit: degree 3 Volume range: 725.00 - 769.00, points=16
The thermal expansion at T = 300K and P = 0 GPa is:
thermal_exp_p(300, 0)
Thermal expansion: 3.14e-05 K^-1 Bulk modulus: 161.36 GPa Pressure: 0.00 GPa Volume: 767.5342 A^3
In this calculation Kp (which is an intermediate result of the algorithm) can be optimized by using the optional keyword fix:
thermal_exp_p(300, 0, fix=0.)
Thermal expansion: 3.10e-05 K^-1 Bulk modulus: 163.17 GPa Pressure: 0.00 GPa Volume: 767.4477 A^3
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.
alpha_serie(20, 600, 18, 0, prt=False, fit=True)
Warning: volume out of range; reduce temperature
array([ 6.80649721e-05, -9.01011534e-03, 2.59751556e-03, -7.83666189e-04])
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.
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:
entropy_p(300, 0)
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
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$:
cp(300,0, prt=True)
Cp: 335.10, Cv: 329.60, Cp/Cv: 1.01669, alpha: 3.136e-05, K: 161.36
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):
cp_serie(40,400,12, 0, prt=True)
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
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):
cp_serie(40,400,18, 0, prt=False, graph=False)
array([ 1.17314930e+03, -1.60123258e-04, -6.97102652e+05, -4.07579977e-02, 9.32201845e+04, -1.93065830e+04])
Such values of the coefficient can be saved in a variable like:
cp_coeff=cp_serie(40,400,18, 0, prt=False, graph=False)
and used to compute $C_P$ at any temperature by using the polynomial fitting function cp_fun:
print(cp_coeff)
round(cp_fun(300, *cp_coeff),2)
[ 1.17314930e+03 -1.60123258e-04 -6.97102652e+05 -4.07579977e-02 9.32201845e+04 -1.93065830e+04]
334.83
In the last example, the Python round command has been used to approximate the printed value returned by cp_fun to 2 decimal digits.
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):
lo.off()
LO-TO splitting not taken into account
To see the effect of the correction on the EoS parameters in this case:
reset_fix()
eos_temp(300,prt=False)
** 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
which is to be compared with:
lo.on()
eos_temp(300, prt=False)
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
The bulk modulus is not significantly affected by such LO-TO correction; however thermodynamic properties can be changed:
lo.off()
cp(300,0,prt=True)
LO-TO splitting not taken into account Cp: 335.85, Cv: 330.30, Cp/Cv: 1.01681, alpha: 3.135e-05, K: 163.02
to be compared with:
lo.on()
cp(300,0,prt=True)
Frequencies corrected for LO-TO splitting Cp: 335.09, Cv: 329.64, Cp/Cv: 1.01654, alpha: 3.104e-05, K: 163.17
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:
kieffer.on()
kieffer.freq(300., 350., 400.)
kieffer.plot()
Kieffer correction on
At any temperature, the contribution to $F$ from the acoustic branches can be recovered with the method kieffer.get_value; for instance:
kieffer.get_value(300)
array(-26.37303978)
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.:
eos_temp(300, prt=False)
** 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
On the other hand, other properties can be significantly affected:
cp(300,0,prt=True)
Cp: 334.13, Cv: 328.68, Cp/Cv: 1.01658, alpha: 3.104e-05, K: 163.17
to be compared with:
kieffer.off()
cp(300,0,prt=True)
Kieffer correction off Cp: 328.94, Cv: 323.49, Cp/Cv: 1.01685, alpha: 3.104e-05, K: 163.17
Functions exist to compute and to plot the Gibbs free energy $G=F+PV$ as function of $T$ or $P$; for instance:
set_fix(4.39)
gibbs_serie_t(200,400,10,0., v0=py.v0, g0=py.g0)
Temp (K) Gibbs energy (J/mol) 200.00 -5913502.3 222.22 -5917192.6 244.44 -5921461.8 266.67 -5926302.1 288.89 -5931703.0 311.11 -5937651.9 333.33 -5944135.5 355.56 -5951139.7 377.78 -5958650.3 400.00 -5966653.0
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):
py.info()
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
To compute the $G$ energy in the range $[0, 10GPa]$ of pressure (10 points), at T=300K, use the function gibbs_serie_p:
gibbs_serie_p(0, 10, 10, 300., v0=py.v0, g0=py.g0)
Press (GPa) Gibbs energy (J/mol) 0.00 -5934609.7 1.11 -5809264.5 2.22 -5684735.1 3.33 -5560993.9 4.44 -5438015.2 5.56 -5315774.7 6.67 -5194249.6 7.78 -5073418.5 8.89 -4953261.2 10.00 -4833758.5
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:
equilib(300,600,12, prod=['py',1], rea=['ens', 1.5, 'cor',1])
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
Average Clapeyron Slope (from Delta S/Delta V): -32.64 bar/K Clapeyron slope (from a linear fit of the P/T curve): -32.75 bar/K
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:
reset_fix()
set_volume_range(725,772)
info.show()
upload_mineral(200,600, mqm='py', volc=True)
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 K Bulk Modulus: 163.22 (0.00) GPa Kp: 4.11 (0.00) V0: 767.4041 (0.000) A^3 Fix status: False Fit of frequencies is active Polynomial fit: degree 3 Volume range: 725.00 - 772.00, points=16 Results from the fit (from high to low order) [-3.48e-02 1.74e+02]
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:
py.info()
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
The new equilibrium curve for the reaction, in the P/T space, is then:
equilib(300,600,12, prod=['py',1], rea=['ens', 1.5, 'cor',1])
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
Average Clapeyron Slope (from Delta S/Delta V): -25.14 bar/K Clapeyron slope (from a linear fit of the P/T curve): -25.10 bar/K
which is in very good agreement with that derived by using the Holland & Powell data for pyrope.