A new class (thermal_expansion_class) has been defined which is an interface to the calculation of thermal expansion by using three different methods.
Let's see how it works. As usual, run an instance of the bm3_thermal_2 module (the example concerns the mineral pyrope):
%matplotlib inline
%run bm3_thermal_2.py
This is BMx-QM program, version 2.5.0 - 27/10/2021 Run time: 2021-10-28 17:35:22.856883 File quick_start.txt found in the master folder Input files in 'pyrope' folder Input file STATIC static_new.dat VOLUME volume_agosto_2021.dat FREQ freq_agosto_2021.dat EXP experimental.txt LO LO.txt FITVOL POLY 725. 773. 16 2. FU 4 20 SET 0 1 2 3 4 5 6 CP 0. 2. -2. 1. -1. -0.5 ALPHA 0. -2 -1 -0.5 END -------- End of input file ------- Frequencies corrected for LO-TO splitting Static BM3 EoS Bulk Modulus: 173.89 (0.07) GPa Kp: 4.26 (0.02) V0: 756.2048 (0.00) A^3 E0: -1.14182291e+04 (1.89e-06) hartree
Make some preliminary calculation of the equation of state; for instance
eos_temp(298.15, prt=False)
** BM3 fit ** EoS at the temperature of 298.15 K Bulk Modulus: 159.30 (0.13) GPa Kp: 4.96 (0.03) V0: 767.5264 (0.006) A^3 Fix status: False Fit of frequencies is active Polynomial fit: degree 2 Volume range: 725.00 - 773.00, points=16
Set suitable parameters for the calculation of numerical derivatives of the pressure with respect to the temperature. This can be done by using the adaptive scheme defined in the delta_ctrl driver:
delta_ctrl.adaptive_on()
Print a summary of the current parameters and conditions adopted for the calculations:
info.show()
Current settings and results
Static data ** min, max volumes: 718.3435, 771.2740; points: 10
Frequency volume range ** min, max volumes: 725.7483, 772.6196; points: 7
Selected freq. sets ** min, max volumes: 725.7483, 772.6196; points: 7
Frequency sets: [0, 1, 2, 3, 4, 5, 6]
Fit of frequencies ** type: poly, degree: 2
min, max volumes: 725.0000, 773.0000; points 16
*** Static EoS (BM3) ***
K0: 173.89 GPa, Kp: 4.26, V0: 756.2048 A^3
** BM3 EoS from the last computation, at the temperature of 298.15 K **
K0: 159.30 GPa, Kp: 4.96, V0: 767.5264 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
Upgrade shift: False
**** Volume driver for volume_from_F function ****
In addition to the attributes set in the parent volume_control_class:
shift: 0.0, flag: False, upgrade_shift: False
**** Numerical T-derivatives driver class (delta_ctrl) ****
Adaptive scheme active:
T_min, T_max: 50.0, 1000.0 K
Delta_min, Delta_max: 10.0, 50.0 K
Degree: 4
N. of points 8
Now let's calculate the thermal expansion. The method thermal_expansion.info provides some details concerning the different algorithms:
thermal_expansion.info()
Methods currently implemented
k_alpha_dir: computes alpha from the product K*alpha, through the
derivative of P with respect to T, at constant V
At any T and P, K and P are directly computed from
the Helmholtz free energy function derivatives. No EoS
is involved at any step;
k_alpha_eos: same as k_alpha_dir, but pressures and bulk moduli
are computed from an EoS;
alpha_dir: the computation is perfomed through the derivative
of the unit cell volume with respect to V; volumes are
calculated without reference to any EoS, by the function
volume_dir.
The function to be used is thermal_expansion.compute:
help(thermal_expansion.compute)
Help on method compute in module __main__:
compute(tt, pp, method='default', fix=0, prt=False) method of __main__.thermal_expansion_class instance
Thermal expansion at a specific temperature and pressure
Args:
tt: temperature (K)
pp: pressure (GPa)
method: 3 methods are currently implemented ('k_alpha_dir',
'k_alpha_eos' and 'alpha_dir'); default 'k_alpha_dir'
fix: relevant for method 'k_alpha_eos' (default 0., Kp not fixed)
prt: relevant for method 'k_alpha_eos'; it controls printout
(default False)
At T=298.15 and P=0 GPa, the thermal expansion of pyrope is estimated (default method k_alpha_dir) to be:
thermal_expansion.compute(298.15, 0, prt=True)
Thermal expansion: 3.24e-05 (K^-1) Bulk modulus: 158.06 (GPa) Volume: 767.5381 (A^3)
By the method k_alpha_eos, the computation provides:
thermal_expansion.compute(298.15, 0, method='k_alpha_eos', prt=True)
Thermal expansion: 3.22e-05 K^-1 Bulk modulus: 159.30 GPa Pressure: 0.00 GPa Volume: 767.5264 A^3
Note that with method k_alpha_dir, the computation of bulk modulus in done by the bulk_modulus_p function, with the argument noeos set to True, and the equilibrium volume is from the volume_dir function.
On the contrary, with the method k_alpha_eos, bulk modulus and equilibrium volume come from an EoS fit at the select temperature. Pressures at constant volume are computed from the refined EoS at the different temperature in an interval centered on the selected temperature.
The third method alpha_dir is the least efficient in term of computing time, but it does not require the estimation of the bulk modulus by whatever method: it computes $\alpha$ as
$$\alpha = \frac{1}{V} \left(\frac{\partial V}{\partial T}\right)_P$$The volume at each temperature in an interval centered on the wanted temperature, is compute by the volume_dir function. In this case, the algorithm provides:
thermal_expansion.compute(298.15, 0, method='alpha_dir', prt=True)
Thermal expansion: 3.24e-05 K^-1
Thermal expansion can also be computed in a range of temperatures, and optionally fitted by a power series whose parameters can be uploaded in the internal thermodynamic database if a phase label is provided. This can be done by using the compute_serie method:
help(thermal_expansion.compute_serie)
Help on method compute_serie in module __main__:
compute_serie(tmin, tmax, pressure=0, nt=0, fit='default', tex='default', title='default', save='default', phase='default', method='default', prt=True, fix=0) method of __main__.thermal_expansion_class instance
Thermal expansion in a T-range
Args:
tmin, tmax: minimum and maximum temperature in the range
pressure: pressure (GPa); default 0
nt: number of points in the T-range; if nt=0, the default is chosen (12)
method: one of the three methods currently implemented
fit: if True, a power series fit is performed
phase: if fit is True and a phase name is specified (label), the data
from the power series fit are loaded in the internal database
fix: relevant for the method 'k_alpha_eos'; if fix is not 0.,
Kp is fixed at the specified value
title: if True, a title of the plot is provided
tex: if tex is True, laTeX formatting is provided
prt: relevant for the method 'k_alpha_eos'
save: if True, the plot is saved in a file
Note:
if save is True and method is 'k_alpha_eos', the name of the file
where the plot is saved is controlled by the plot.name and plot.ext variables.
The file resolution is controlled by the plot.dpi variable.
The appropriate parameters can be set by the set_param method
of the plot instance of the plot_class class (in the plot.py module)
Example:
>>> plot.set_param(dpi=200, name='alpha_k_eos_serie')
>>> thermal_expansion.compute_serie(100, 500, method='k_alpha_eos', save=True)
Here the example for pyrope (default method 'k_alpha_dir')
thermal_expansion.compute_serie(150, 400, fit=True, phase='py')
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: +1.0438e-04 +0.0 -3.9410e-01 -2.0 +1.5243e-02 -1.0 -2.0496e-03 -0.5
Note that, in this case, all of the thermodynamic parameters for pyrope are read from an external database and only the thermal expansion data are upgraded to the new computed values.
By writing a few lines of python code, it is for instance possible to compare the thermal expansion computed through two different algorithms like k_alpha_dir and k_alpha_eos. The following lines computes the coefficients of the power series from the two methods:
tmin, tmax = 100, 450
alpha_fit1=thermal_expansion.compute_serie(tmin, tmax, fit=True)
alpha_fit2=thermal_expansion.compute_serie(tmin, tmax, method='k_alpha_eos',fit=True)
The coefficients of the power series optimized for the two cases are then used to compute $\alpha$ by using the alpha_fun function that accepts a list of temperature (t_list) and the list of coefficients as arguments:
t_list=np.linspace(tmin, tmax, 50)
alpha_list_1=np.array([alpha_fun(it, *alpha_fit1) for it in t_list])
alpha_list_2=np.array([alpha_fun(it, *alpha_fit2) for it in t_list])
The two $\alpha$ curves are then plotted together by using the plot.multi method:
x=[t_list, t_list]
y=[alpha_list_1, alpha_list_2]
style=['k-', 'r--']
label=['k\_alpha\_dir', 'k\_alpha\_eos']
plot.set_param(fsize=14, tsize=12)
plot.multi(x,y,style,label, xlab='Temperature (K)', ylab=r'$\alpha$ (K$^{-1}$)', title='Thermal expansion', tex=True)
Use the help function to get the: documentation concerning the plotting facilities:
help(plot.multi)
Help on method multi in module plot:
multi(x, y, sl, lgl, title='', xlab='default', ylab='default', fsize=0, tsize=0, save=False, dpi=0, name='default', ext='default', tex='default') method of plot.plot_class instance
Plot multiple curves.
Args:
x: array of independent variable arrays
y: array of dependent variable arrays
sl: array of styles (one for each curve)
lgl: array of curve labels (for the legend)
Note:
See the documentation for the 'simple' method for information
about the other keyword arguments.
Example:
>>> x=np.linspace(-2, 2, 24)
>>> y1=x**2
>>> y2=x**3
>>> x=[x,x]
>>> y=[y1,y2]
>>> style=['k-', 'r+']
>>> label=[r'Y=X$^2$', r'Y=X$^3$']
>>> plot.multi(x,y, style, label, tex=True)
help(plot.simple)
Help on method simple in module plot:
simple(x, y, title='', style='default', xlab='default', ylab='default', fsize=0, tsize=0, save=False, dpi=0, name='default', ext='default', tex='default') method of plot.plot_class instance
Plots a single (x, y) curve.
Args:
x: independent variable (array of data)
y: dependent variable (array of data)
title: plot title (default no title '')
style: style of the plot (default self.style)
xlab: X axis label (default self.xlab)
ylab: Y axis label (default self.ylab)
fsize: label and title size (default self.fsize)
tsize: tick size (default self.tsize)
save: if True, the plot will be saved on a file (default False)
dpi: resolution of the saved file (default self.dpi)
name: name of the saved file (default self.name)
ext: extension (and format) of the saved file (default self.ext)
tex: if True, tex fonts and format are used
Note:
The arguments dpi, name and ext are relevant only if save is True
Examples:
>>> plot.simple(x,y, style='r-', xlab='T (K)', ylab='Alpha (K^-1)', fsize=12)
>>> plot.set_param(fsize=12, style='r-')
>>> plot.simple(x,y, xlab='T (K)', ylab='alpha (K^-1)')