Update file

Current version: 1.0.0 - 03/03/2021

Significant Changes to some functions of the older versions

The main change concerns the possibility to extrapolate the specific heat at constant volume ($C_P$) in the high temperature range, by imposing the Dulong-Petit limit for $C_V$ at some chosen high temperature value. The extrapolation is used to get high temperature values of $C_P$ (specific heat at constant pressure) and $\alpha$ (thermal expansion).

The affected functions are

  • cp_serie
  • alpha_serie
  • upload_mineral

Documentation is provided below, with reference to the pyrope mineral.

In [1]:
%matplotlib inline
%run bm3_thermal_2.py


This is BM3-thermal program, version 1.0.0 - 03/03/2021 
Run time:  2021-03-04 16:34:46.277770

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. 773. 12 3
FU
4 20
SET
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
TEMP
300
CP
0. 2. -2. 1. -1. -0.5 -3.
ALPHA
0. -2  -1 -0.5
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

The thermodynamic information concerning pyrope as stored in the Holland & Powell database, are:

In [2]:
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

These data are used to compute the equilibrium curve for the reaction

$${\rm 3/2 \ Mg_2Si_2O_6 + Al_2O_3 \ \leftrightarrow \ Mg_3Al_2Si_3O_{12}}$$

(entirely from Holland & Powell data)

In [3]:
equilib(300,1000,12)

  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    
  363.6   3.46     6483.589     17.830       -0.534       -33.38    
  427.3   3.25     7659.187     17.926       -0.542       -33.10    
  490.9   3.04     8770.490     17.866       -0.549       -32.53    
  554.5   2.83     9787.076     17.649       -0.557       -31.69    
  618.2   2.63    10694.223     17.299       -0.565       -30.62    
  681.8   2.44    11485.327     16.845       -0.573       -29.40    
  745.5   2.26    12158.208     16.310       -0.581       -28.07    
  809.1   2.09    12713.185     15.713       -0.589       -26.67    
  872.7   1.92    13152.013     15.070       -0.597       -25.23    
  936.4   1.77    13477.272     14.393       -0.605       -23.78    
 1000.0   1.62    13692.003     13.692       -0.613       -22.33    



Average Clapeyron Slope (from Delta S/Delta V): -29.20 bar/K
Clapeyron slope (from a linear fit of the P/T curve): -29.53 bar/K

Now, thermodynamic data for pyrope are computed from quantum-mechanical properties as specified in the input file of the program. For the purpose, the upload_mineral function is used:

In [7]:
help(upload_mineral)

Help on function upload_mineral in module __main__:

upload_mineral(tmin, tmax, points=12, HT_lim=0.0, deg=1, g_deg=1, model=1, mqm='py', b_dir=False, blk_dir=False, extra_alpha=True, volc=False)
    Prepares data to be uploaded in the mineral database.
    
    Args:
        tmin:   minimum temperature for fit K0, Cp and alpha
        tmax:   maximum temperature for fit
        points: number of points in the T range for fit
        mqm:    name of the mineral, as specified in the internal 
                database,
        b_dir:  if True, the alpha_dir_serie function is used for the
                computation of thermal expansion
        blk_dir:  if True, the bulk_dir_serie function is used
                to compute the bulk modulus as a function of T;
                if False, the function bulk_serie is used.
        HT_lim: Temperature at which the Dulong-Petit limit for Cv
                is supposed to be reached (default 0.: no Dulong-Petit
                                           model)
        model:  Used in the HT_limit estimation of Cv; Einstein model
                for Cv(T) with one frequency (default model=1), or with
                2 frequencies (model=2)
        deg:    Used in the HT limit estimation of alpha (relevant if
                HT_lim > 0; default: 1)
        g_deg:  Used in the HT limit estimation of Cp (relevant if
                HT_lim > 0.; default 1)
        volc:   if True, V0 is set at the value found in the database
                (default: False)
        
                
    Example:
        upload_mineral(300,800,16,mqm='coe', b_dir=True)

Here we issue the command by specifying

  • tmin=200;
  • tmax=450 (tmax corresponds to a cell volume close to the maximum one at which static energies and frequencies are effectively computed);
  • $C_P$ and $\alpha$ will be computed at 12 points in the $[tmin, tmax]$ range;
  • HT_lim=4000 (the temperature at which the DP limit is supposed to be reached);
  • model=2: a two-temperature Einstein model is used to fit the computed $C_V$ data; the Einstein model fitted in the $[tmin, tmax]$ range (augmented with the DP datum at HT_lim) is then used to extrapolate $C_V$ at higher temperature;
  • g_deg=1: $C_P/C_V$ is fitted as a linear function of $T$; this fit, done in the $[tmin, tmax]$ range, is then used to get $C_P$ from $C_V$ at higher temperatures, up to HT_lim;
  • deg=1: the $\gamma/VK_T$ function is fitted as a linear function of $T$ in the $[tmin, tmax]$ range; such fit is then used to get an extrapolation of $\alpha$ at high temperatures up to HT_lim;
  • extra_alpha=True: the $\alpha$ extrapolation is required.

The other parameters are identical to those of the older upload_mineral function.

In [4]:
upload_mineral(200,450,12,HT_lim=4000, deg=1, g_deg=1, model=2, extra_alpha=True, volc=True)

 ** BM3 fit **

EoS at the temperature of 298.15 K

Bulk Modulus: 159.45   (0.14) GPa
Kp:             4.89   (0.03)
V0:           767.6153 (0.006) A^3

Fix status: False
Fit of frequencies is active
Polynomial fit: degree  3
Volume range: 725.00 - 773.00, points=12

Results from the fit (from high to low order)
[-3.98e-02  1.71e+02]

*** High temperature Cp estimation from the Dulong-Petit limit

 T limit: 4000.00

 ** BM3 fit **

EoS at the temperature of 298.15 K

Bulk Modulus: 159.45   (0.14) GPa
Kp:             4.89   (0.03)
V0:           767.6153 (0.006) A^3

Fix status: False
Fit of frequencies is active
Polynomial fit: degree  3
Volume range: 725.00 - 773.00, points=12
Kp fixed to 4.89

Einstein temperature
empirical estimation (from molar entropy): 538.88 K
result from fit:                           323.72, 1064.93 K
Dulong-Petit limit for Cv (T = 4000.00 K): 498.87 J/mol K

  T (K)  Cv "real"  Cv "fit"
 200.00   234.53     236.21 
 227.78   263.93     263.16 
 255.56   289.35     287.85 
 283.33   311.66     310.12 
 311.11   331.01     329.95 
 338.89   348.07     347.47 
 366.67   362.72     362.86 
 394.44   374.81     376.35 
 422.22   386.34     388.17 
 450.00   396.07     398.53 

Fit at T = 4000.00 K: Cv = 497.26 J/mol K

Gamma estimation (extrapolation from lower T values)

T:   470.00  Cv:   405.19  gamma: 1.030  Cp:   417.38
T:   705.33  Cv:   452.10  gamma: 1.047  Cp:   473.36
T:   940.67  Cv:   471.40  gamma: 1.064  Cp:   501.56
T:  1176.00  Cv:   480.93  gamma: 1.081  Cp:   519.84
T:  1411.33  Cv:   486.27  gamma: 1.098  Cp:   533.86
T:  1646.67  Cv:   489.55  gamma: 1.115  Cp:   545.76
T:  1882.00  Cv:   491.70  gamma: 1.132  Cp:   556.49
T:  2117.33  Cv:   493.19  gamma: 1.149  Cp:   566.53
T:  2352.67  Cv:   494.26  gamma: 1.166  Cp:   576.13
T:  2588.00  Cv:   495.05  gamma: 1.183  Cp:   585.45
T:  2823.33  Cv:   495.66  gamma: 1.200  Cp:   594.57
T:  3058.67  Cv:   496.13  gamma: 1.216  Cp:   603.54
T:  3294.00  Cv:   496.51  gamma: 1.233  Cp:   612.41
T:  3529.33  Cv:   496.81  gamma: 1.250  Cp:   621.21
T:  3764.67  Cv:   497.06  gamma: 1.267  Cp:   629.94
T:  4000.00  Cv:   497.26  gamma: 1.284  Cp:   638.63

---- High temperature estimation of the thermal expansion coefficient ----


 ** BM3 fit **

EoS at the temperature of 298.15 K

Bulk Modulus: 159.45   (0.03) GPa
Kp:             4.89   (0.00)
V0:           767.6153 (0.004) A^3

Fix status: True
Kp fixed at 4.89
Fit of frequencies is active
Polynomial fit: degree  3
Volume range: 725.00 - 773.00, points=12

Dulong-Petit limit of Cv 498.87 (J/mol K)
Einstein's temperature: 538.88 (K)
Gamma(V) fit from already stored data

Gamma/VK fit of degree 1
Alpha constrained at the high temperature limit: 4.80e-05 K^-1

 -------------------------------

The computed thermodynamic data, at the QM level, for pyrope are:

In [8]:
py.info()

Mineral: pyrope

K0: 159.45 GPa, Kp: 4.89, dK0/dT: -0.0398 GPa/K, V0: 11.3180 J/bar
G0: -5934105.00 J/mol, S0: 265.94 J/mol K

EoS type: bm 
Cp coefficients and powers:
+3.4537e+02  +0.0
-3.9308e-07  +2.0
+1.1340e+07  -2.0
+4.7702e-02  +1.0
-2.2814e+05  -1.0
+1.0457e+04  -0.5
+1.0000e+00  -3.0

Alpha coefficients and powers:
+4.5995e-05  +0.0
+6.9899e-01  -2.0
-1.1873e-02  -1.0
+3.1157e-04  -0.5

The computed equilibrium curve for the afore mentioned reaction is now:

In [6]:
equilib(300,1000,12)

  T (K)  P (GPa)  DH(J/mol)  DS (J/mol K)  DV (J/bar)  Slope (bar/K)
  300.0   3.61     3676.842     12.256       -0.546       -22.44    
  363.6   3.46     4617.715     12.699       -0.542       -23.41    
  427.3   3.31     5667.622     13.265       -0.539       -24.62    
  490.9   3.15     6872.608     14.000       -0.535       -26.15    
  554.5   2.98     8220.964     14.825       -0.532       -27.85    
  618.2   2.79     9687.742     15.671       -0.529       -29.61    
  681.8   2.60    11248.308     16.498       -0.527       -31.33    
  745.5   2.40    12882.143     17.281       -0.524       -32.97    
  809.1   2.18    14573.454     18.012       -0.522       -34.50    
  872.7   1.96    16310.781     18.689       -0.520       -35.93    
  936.4   1.72    18086.371     19.316       -0.518       -37.25    
 1000.0   1.48    19895.585     19.896       -0.517       -38.49    



Average Clapeyron Slope (from Delta S/Delta V): -30.38 bar/K
Clapeyron slope (from a linear fit of the P/T curve): -30.38 bar/K

where te relevant thermodynamic data for enstatite ($\rm Mg_2Si_2O_6$) and corundum ($\rm Al_2O_3$) are still from the Holland & Powell database