Thermodynamics of kyanite is here evaluated from frequencies and static energies evaluated ab initio, and processed within the framework offered by the Quasi-Harmonic approximation (QHA)
The modified Kieffer model is used to estimate the contribution of the acoustic modes to the Helmholtz free energy
The following command is a directive issued in order to show plots directly in the current notebook, instead of being done in a separate windows
%matplotlib inline
In the following cell, the program is loaded; the input file is shown and the Kieffer model is applied to estimate the free energy contribution from the acoustic modes, in a given temperature range. The $F_{acoustic}(T)$ function, in numerical form, is stored in a stack from which values of $F$, at any temperature in the range, are recovered by numerical interpolation of the stored values by the method get_value of the kieffer class
%run bm3_thermal_2.py
This is BMx-QM program, version 2.4.0 - 18/02/2021 Run time: 2021-05-15 19:25:58.725099 File quick_start.txt found in the master folder Input files in 'kyanite' folder Input file TITLE Kyanite STATIC static.txt VOLUME volumes.txt FREQ frequencies.txt EXP experimental.txt FITVOL POLY 283. 302. 16 2 FU 4 8 SET 0 1 2 3 4 5 6 7 CP 0. 1. -2. -0.5 -1. -3. ALPHA 0. -1. -2. -0.5 KIEFFER 200. 300. 300. END -------- End of input file ------- Warning ** Volume and frequencies lists have been sorted indexing: [7 6 5 4 3 2 1 0] Static BM3 EoS Bulk Modulus: 184.19 (0.03) GPa Kp: 3.95 (0.01) V0: 298.0788 (0.00) A^3 E0: -4.59916109e+03 (3.01e-07) hartree Kieffer model for the acoustic branches activated
The info.show command displays various parameters set in the input.txt file
info.show()
Kyanite Current settings and results Static data ** min, max volumes: 280.1477, 303.9919; points: 11 Frequency volume range ** min, max volumes: 283.6984, 301.4423; points: 8 Selected freq. sets ** min, max volumes: 283.6984, 301.4423; points: 8 Frequency sets: [0, 1, 2, 3, 4, 5, 6, 7] Fit of frequencies ** type: poly, degree: 2 min, max volumes: 283.0000, 302.0000; points 16 ** Static EoS (BM3) ** K0: 184.19 GPa, Kp: 3.95, V0: 298.0788 A^3 No zone center excluded modes No off center excluded modes Kieffer model on; frequencies 200.00 300.00 300.00 cm^-1 **** 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
A 3^rd order Birch-Murnaghan EoS (BM3) is computed at 300K; this is done by fitting a volume integrated form of the EoS to the $F(V)$ numerical values.
eos_temp(300)
** BM3 fit ** EoS at the temperature of 300.00 K Bulk Modulus: 176.51 (0.01) GPa Kp: 4.00 (0.00) V0: 301.4421 (0.000) A^3 Fix status: False Fit of frequencies is active Polynomial fit: degree 2 Volume range: 283.00 - 302.00, points=16
Volume-Pressure list at 300.00 K 283.000 12.64 284.267 11.64 285.533 10.67 286.800 9.71 288.067 8.77 289.333 7.86 290.600 6.96 291.867 6.08 293.133 5.22 294.400 4.37 295.667 3.55 296.933 2.74 298.200 1.95 299.467 1.18 300.733 0.42 302.000 -0.33
The same result is obtained by fitting the BM3 EoS to the $P(V)$ set of values at 300K; pressures $P$ at each volume are computed as derivatives of $F$ with respect to $V$, at constant temperature:
bulk_dir(300)
BM3 EoS from P(V) fit K0: 176.50 (0.01) GPa Kp: 4.00 (0.00) V0: 301.4423 (0.00) A^3
A comparison of computed $C_P$ and $S$ with available experimental data is possible; $K_p$ is fixed at the value estimated at 300K
set_fix(4.00)
compare_exp()
Warning: volume out of range; reduce temperature Temp Cp exp Cp calc Del Cp S exp S calc Del S 298.15 121.60 123.26 1.66 82.80 83.25 0.45 300.00 122.28 123.85 1.57 83.55 84.02 0.46 350.00 137.79 138.21 0.42 103.63 104.23 0.60 400.00 149.27 149.42 0.15 122.82 123.45 0.63 450.00 158.16 158.28 0.12 140.93 141.58 0.65 500.00 165.27 165.08 -0.18 157.98 158.63 0.66 550.00 171.10 171.26 0.16 174.01 174.67 0.65 600.00 175.99 175.89 -0.10 189.12 189.75 0.64 Average error on Cp: 0.55; maximum error: 1.66 Average error on S: 0.59; maximum error: 0.66
Warning: issue on volume repeated 363 times
In the following, the equilibrium kyanite <--> andalusite will be analysed. To begin with, the equilibrium curve in the P/T space is calculated by using the Holland & Powell (H&P) thermodynamics data. Such data for kyanite are shown by the method info of the ky instance of the class mineral:
ky.info()
Mineral: kyanite K0: 159.00 GPa, Kp: 4.00, dK0/dT: -0.0238 GPa/K, V0: 4.4140 J/bar G0: -2442541.00 J/mol, S0: 83.50 J/mol K EoS type: m Cp coefficients and powers: +2.7940e+02 +0.0 -7.1240e-03 +1.0 -2.0556e+06 -2.0 -2.2894e+03 -0.5 Alpha coefficients and powers: +4.0400e-05 +0.0 -4.0400e-04 -0.5
The equilibrium curve is computed by the command equilib in the 500 - 700K temperature range (12 points)
equilib(500, 700, 12, prod=['ky',1], rea=['andal',1])
T (K) P (GPa) DH(J/mol) DS (J/mol K) DV (J/bar) Slope (bar/K) 500.0 0.04 -4678.303 -9.357 -0.742 12.61 518.2 0.06 -4841.963 -9.344 -0.742 12.59 536.4 0.08 -5004.423 -9.330 -0.742 12.57 554.5 0.10 -5165.743 -9.315 -0.742 12.55 572.7 0.13 -5325.980 -9.299 -0.742 12.53 590.9 0.15 -5485.183 -9.283 -0.742 12.50 609.1 0.17 -5643.400 -9.265 -0.742 12.48 627.3 0.20 -5800.673 -9.247 -0.743 12.45 645.5 0.22 -5957.041 -9.229 -0.743 12.43 663.6 0.24 -6112.542 -9.211 -0.743 12.40 681.8 0.26 -6267.208 -9.192 -0.743 12.37 700.0 0.29 -6421.072 -9.173 -0.743 12.35
Average Clapeyron Slope (from Delta S/Delta V): 12.49 bar/K Clapeyron slope (from a linear fit of the P/T curve): 12.49 bar/K
Now, thermodynamics data for kyanite are computed from quantum-mechanical (QM) calculations. These regard:
The volume at the reference state is the experimental one: as it is well known, quantum-mechanical calculations generally overestimate cell volumes; this fact leads to an overestimation of the $PV$ contribution to the Gibbs free energy ($G=F+PV$) which generally becomes significant at high pressures.
The $G$ free energy in the reference ($G_0$) is also from the H&P database: the $G$ energy computed at the QM level is some absolute energy of the mineral phase, whereas the corresponding quantity from the H&P database is the energy of formation from the elements. In order to properly fix the zero of the energy scale (that is, to put both kyanite and andalusite on the same energy scale), the $G_0$ must be properly set.
All of such calculations, together with the loading of such quantities in the instance ky of the mineral class, are performed by the upload_mineral command:
upload_mineral(250,700,mqm='ky', volc=True)
** BM3 fit ** EoS at the temperature of 298.15 K Bulk Modulus: 176.54 (0.01) GPa Kp: 4.00 (0.00) V0: 301.4329 (0.000) A^3 Fix status: False Fit of frequencies is active Polynomial fit: degree 2 Volume range: 283.00 - 302.00, points=16 Results from the fit (from high to low order) [-1.87e-02 1.82e+02]
Warning: volume out of range; reduce temperature
Now, ky.info will show the new data stored in the database, for kyanite:
ky.info()
Mineral: kyanite K0: 176.54 GPa, Kp: 4.00, dK0/dT: -0.0187 GPa/K, V0: 4.4140 J/bar G0: -2442541.00 J/mol, S0: 83.25 J/mol K EoS type: bm Cp coefficients and powers: -2.9986e+02 +0.0 +1.2169e-01 +1.0 +1.2437e+07 -2.0 +2.0340e+04 -0.5 -2.7762e+05 -1.0 +1.0000e+00 -3.0 Alpha coefficients and powers: +3.0117e-05 +0.0 -6.5918e-03 -1.0 +3.1400e-01 -2.0 +8.5304e-05 -0.5
equilib(500, 700, 12, prod=['ky',1], rea=['andal',1])
T (K) P (GPa) DH(J/mol) DS (J/mol K) DV (J/bar) Slope (bar/K) 500.0 0.04 -4685.225 -9.370 -0.742 12.63 518.2 0.06 -4851.442 -9.362 -0.742 12.61 536.4 0.08 -5018.134 -9.356 -0.742 12.61 554.5 0.11 -5185.431 -9.351 -0.742 12.60 572.7 0.13 -5353.409 -9.347 -0.742 12.59 590.9 0.15 -5522.089 -9.345 -0.742 12.59 609.1 0.18 -5691.441 -9.344 -0.742 12.59 627.3 0.20 -5861.383 -9.344 -0.742 12.59 645.5 0.22 -6031.784 -9.345 -0.742 12.59 663.6 0.24 -6202.473 -9.346 -0.743 12.59 681.8 0.27 -6373.234 -9.347 -0.743 12.59 700.0 0.29 -6543.817 -9.348 -0.743 12.59
Average Clapeyron Slope (from Delta S/Delta V): 12.60 bar/K Clapeyron slope (from a linear fit of the P/T curve): 12.59 bar/K
The impact of the contribution form acoustic modes for the present case is well shown by switching off the Kieffer model and by repeating the computation:
kieffer.off()
upload_mineral(250, 700, mqm='ky', volc=True)
ky.info()
equilib(500, 700, 12, prod=['ky',1], rea=['andal',1])
Kieffer correction off ** BM3 fit ** EoS at the temperature of 298.15 K Bulk Modulus: 176.54 (0.01) GPa Kp: 4.00 (0.00) V0: 301.4329 (0.000) A^3 Fix status: False Fit of frequencies is active Polynomial fit: degree 2 Volume range: 283.00 - 302.00, points=16 Results from the fit (from high to low order) [-1.87e-02 1.82e+02]
Mineral: kyanite K0: 176.54 GPa, Kp: 4.00, dK0/dT: -0.0187 GPa/K, V0: 4.4140 J/bar G0: -2442541.00 J/mol, S0: 77.23 J/mol K EoS type: bm Cp coefficients and powers: +1.7667e+02 +0.0 -1.1546e-02 +1.0 +4.8148e+06 -2.0 +3.4013e+03 -0.5 -9.1446e+04 -1.0 +1.0000e+00 -3.0 Alpha coefficients and powers: +3.0620e-05 +0.0 -6.2078e-03 -1.0 +2.9092e-01 -2.0 +5.8886e-05 -0.5 T (K) P (GPa) DH(J/mol) DS (J/mol K) DV (J/bar) Slope (bar/K) 500.0 0.25 -9175.306 -18.351 -0.739 24.83 518.2 0.29 -9608.680 -18.543 -0.739 25.10 536.4 0.34 -10045.599 -18.729 -0.739 25.36 554.5 0.38 -10486.302 -18.910 -0.738 25.61 572.7 0.43 -10931.050 -19.086 -0.738 25.86 590.9 0.48 -11380.109 -19.259 -0.738 26.11 609.1 0.53 -11833.751 -19.429 -0.737 26.35 627.3 0.57 -12292.242 -19.596 -0.737 26.59 645.5 0.62 -12755.847 -19.763 -0.737 26.83 663.6 0.67 -13224.821 -19.928 -0.736 27.06 681.8 0.72 -13699.412 -20.092 -0.736 27.30 700.0 0.77 -14179.858 -20.257 -0.736 27.54
Average Clapeyron Slope (from Delta S/Delta V): 26.21 bar/K Clapeyron slope (from a linear fit of the P/T curve): 26.22 bar/K
Indeed, the equilibrium curve is shifted at significantly higher pressures and the Clapeyron slope significantly increases.
The frequencies used in the Kieffer model can be set by the freq method of the kieffer class (or the input file, under the KIEFFER keyword):
kieffer.on()
kieffer.freq(200, 250, 300)
upload_mineral(250, 700, mqm='ky',volc=True)
ky.info()
equilib(500, 700, 12, prod=['ky',1], rea=['andal',1])
Kieffer correction on ** BM3 fit ** EoS at the temperature of 298.15 K Bulk Modulus: 176.54 (0.01) GPa Kp: 4.00 (0.00) V0: 301.4329 (0.000) A^3 Fix status: False Fit of frequencies is active Polynomial fit: degree 2 Volume range: 283.00 - 302.00, points=16 Results from the fit (from high to low order) [-1.87e-02 1.82e+02]
Mineral: kyanite K0: 176.54 GPa, Kp: 4.00, dK0/dT: -0.0187 GPa/K, V0: 4.4140 J/bar G0: -2442541.00 J/mol, S0: 83.59 J/mol K EoS type: bm Cp coefficients and powers: +2.7928e+02 +0.0 -3.7895e-02 +1.0 +3.2153e+06 -2.0 +6.2067e+00 -0.5 -5.4019e+04 -1.0 +1.0000e+00 -3.0 Alpha coefficients and powers: +3.0494e-05 +0.0 -6.2966e-03 -1.0 +2.9597e-01 -2.0 +6.5213e-05 -0.5 T (K) P (GPa) DH(J/mol) DS (J/mol K) DV (J/bar) Slope (bar/K) 500.0 0.03 -4502.591 -9.005 -0.742 12.13 518.2 0.05 -4656.308 -8.986 -0.742 12.11 536.4 0.07 -4809.587 -8.967 -0.742 12.08 554.5 0.09 -4962.872 -8.949 -0.742 12.05 572.7 0.12 -5116.624 -8.934 -0.742 12.03 590.9 0.14 -5271.317 -8.921 -0.743 12.01 609.1 0.16 -5427.427 -8.911 -0.743 12.00 627.3 0.18 -5585.429 -8.904 -0.743 11.99 645.5 0.20 -5745.791 -8.902 -0.743 11.98 663.6 0.23 -5908.977 -8.904 -0.743 11.99 681.8 0.25 -6075.437 -8.911 -0.743 11.99 700.0 0.27 -6245.614 -8.922 -0.743 12.01
Average Clapeyron Slope (from Delta S/Delta V): 12.03 bar/K Clapeyron slope (from a linear fit of the P/T curve): 12.02 bar/K