New functions and new ways to compute thermal expansion as a function of temperature are implemented in the current release of the program; in addition, new classes are implemented to set several parameters controlling the numerical evaluation of derivatives with respect to the volume and temperature. In particular, have a check to the documentation of the classes
on the index section of the site.
This page illustrates the new features working on the data for pyrope.
Let's load and run the program in the usual way:
%matplotlib inline
%run bm3_thermal_2
This is BMx-QM program, version 2.4.3 - 03/09/2021 Run time: 2021-09-03 22:35:16.876999 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 SPLINE 725. 773. 16 2. 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 KIEFFER 73 73 130 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 Kieffer model for the acoustic branches activated
In order to compute thermal expansion, four methods can be used. Two of them compute the cell volume (V) in a specified range of temperature (T); the logarithm of the V's are then computed and the log(V)/T values are fitted by using a polynomial function or a spline. The last step is to compute the first derivative of the fit to get the thermal expansion $\alpha$ from its definition:
$$\alpha(T) = \frac{d \log(V)}{d T} = \frac{1}{V} \frac{d V}{d T}$$The two methods following this strategy differ by the functions used to get the cell volume at each temperature: one method uses the volume_dir function; the second one uses the volume_from_F function:
the volume_dir function in more general in that it can compute the equilibrium volume at any temperature and any pressure, but it requires the calculation of the pressure as the derivative of the Helmholtz function with respect to V; the computed volume is the one at which the total pressure of the crystal equals the imposed external pressure;
the second method uses the volume_from_F function and computes the volume as the one that minimizes the Helmholtz free energy at the given temperature (and pressure zero); this is equivalent to the minimization of the Gibbs free energy, as the pressure is zero.
Of the other two methods,
one uses the alpha_dir function that computes the volume in a small range of temperatures around the specified temperature, and subsequently makes a numerical derivative of the V(T) values (fitted by a given polynomials); whereas the other method
computes the thermal expansion from the $K\alpha$ product, where $K$ (the bulk modulus) is from an EoS fit at the specified temperature, and
$$K\alpha = \left(\frac{\partial P}{\partial T}\right)_V$$
is computed from the numerical derivative of a polynomial representation of the P(T) function (the pressures at different temperatures are computed by using the corresponding EoS)
The method (4) heavily relies on the quality of the equation of state, whereas method (3) is EoS-free as the pressures needed for the calculation are directly computed as first derivatives of the Helmholtz function, through the pressure_dir function (note however that the static pressure contribution is always from a static EoS).
Appropriate values of the parameters controlling the quality of fits and derivates involved can be set by the control driver function (check the documentation provided); the default parameters should be suitable in the majority of cases (so, no commands must be issued), but in some other cases (the presence of soft-modes producing large thermal expansions, for instance) it might be useful to get acquainted with the volume_ctrl, volume_F_ctrl and delta_ctrl drivers:
volume_ctrl.set_all(upgrade_shift=True)
delta_ctrl.adaptive_on()
delta_ctrl.adaptive_set(dmax=50)
Global information of (nearly) all of the parameters and settings involved in this particular case can be retrieved by the usual info.show command:
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: spline, degree: 2, smooth: 2.0 min, max volumes: 725.0000, 773.0000; points 16 *** Static EoS (BM3) *** K0: 173.89 GPa, Kp: 4.26, V0: 756.2048 A^3 No zone center excluded modes No off center excluded modes Kieffer model on; frequencies 73.00 73.00 130.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 Upgrade shift: True **** 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
Let's compute the thermal expansion of pyrope in the [100, 400 K] T-range (12 points) by using all of the four methods illustrated above. The command is alpha_dir_v; the optional parameter comp is set to True (the default is False) to get the computation done with all the methods available; note that if comp is set to False, only the method (1), described above, is used.
alpha_dir_v(100, 400, 12, comp=True)
Summary of the input parameters T range: 100.0, 400.0 K, Num. of points: 12 Type of Log(V) fit: spline, degree: 4, smooth: 0.0010 Compare with other methods to compute alpha: True Fit alpha values to a power serie: False It is advised to use polynomial fits for 'dir' calculations Fit of frequencies is active Spline fit: degree 2, smooth: 2.0 Volume range: 725.00 - 773.00, points=16
Temp V (1) (2) (3) (4) 100.0 763.7409 1.50e-05 1.50e-05 1.47e-05 1.48e-05 127.3 764.0887 1.85e-05 1.85e-05 1.85e-05 1.85e-05 154.5 764.5081 2.16e-05 2.16e-05 2.16e-05 2.16e-05 181.8 764.9875 2.42e-05 2.42e-05 2.43e-05 2.42e-05 209.1 765.5181 2.65e-05 2.65e-05 2.65e-05 2.65e-05 236.4 766.0935 2.86e-05 2.86e-05 2.85e-05 2.84e-05 263.6 766.7088 3.03e-05 3.03e-05 3.03e-05 3.01e-05 290.9 767.3605 3.20e-05 3.20e-05 3.20e-05 3.17e-05 318.2 768.0456 3.34e-05 3.34e-05 3.35e-05 3.30e-05 345.5 768.7622 3.49e-05 3.49e-05 3.49e-05 3.45e-05 372.7 769.5086 3.63e-05 3.63e-05 3.63e-05 3.57e-05 400.0 770.2841 3.77e-05 3.77e-05 3.76e-05 3.69e-05 (1) from V(T) fit (2) from V(T) from F fit (3) from the definition ('dir' computation) (4) From (dP/dT)_V and EoS
In the figure above, the continuous line corresponds to the method (1); the full circles corresponds to the method (2); star and plus markers (barely visible in the specific case, as they overlap with the circles) corresponds to methods (3) and (4) respectively.
It is possible to compute the thermal expansion in a range of temperature (by the method (1) in the example below), fit it to a power serie and associate the coefficients of the fit to a given mineral in the database. First, compute the thermal expansion and make the fit, saving the coefficients in a given array (alpha_fit, for instance)
alpha_fit=alpha_dir_v(100, 400, 12, comp=False, fit=True)
Summary of the input parameters T range: 100.0, 400.0 K, Num. of points: 12 Type of Log(V) fit: spline, degree: 4, smooth: 0.0010 Compare with other methods to compute alpha: False Fit alpha values to a power serie: True Trim applied to T and alpha values for the power serie fit: 0.0
Temp V Alpha 100.0 763.7409 1.50e-05 127.3 764.0887 1.85e-05 154.5 764.5081 2.16e-05 181.8 764.9875 2.42e-05 209.1 765.5181 2.65e-05 236.4 766.0935 2.86e-05 263.6 766.7088 3.03e-05 290.9 767.3605 3.20e-05 318.2 768.0456 3.34e-05 345.5 768.7622 3.49e-05 372.7 769.5086 3.63e-05 400.0 770.2841 3.77e-05
Second, load the alpha_fit coefficients in the database: for pyrope (py label) just issue the command:
py.load_alpha(alpha_fit, power_a)
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: +9.0687e-05 +0.0 -1.3196e-01 -2.0 +8.4219e-03 -1.0 -1.4677e-03 -0.5
Alternatively, in a single step, the phase can be specified by the optional argument phase:
alpha_dir_v(100, 400, 12, comp=False, fit=True, phase='py')
Summary of the input parameters T range: 100.0, 400.0 K, Num. of points: 12 Type of Log(V) fit: spline, degree: 4, smooth: 0.0010 Compare with other methods to compute alpha: False Fit alpha values to a power serie: True Trim applied to T and alpha values for the power serie fit: 0.0 Power serie coefficient uploaded for the phase py
Temp V Alpha 100.0 763.7409 1.50e-05 127.3 764.0887 1.85e-05 154.5 764.5081 2.16e-05 181.8 764.9875 2.42e-05 209.1 765.5181 2.65e-05 236.4 766.0935 2.86e-05 263.6 766.7088 3.03e-05 290.9 767.3605 3.20e-05 318.2 768.0456 3.34e-05 345.5 768.7622 3.49e-05 372.7 769.5086 3.63e-05 400.0 770.2841 3.77e-05 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: +9.0687e-05 +0.0 -1.3196e-01 -2.0 +8.4219e-03 -1.0 -1.4677e-03 -0.5
Other functions for getting the coefficients of the power serie of the thermal expansion are available; for instance, for method (4):
alpha_fit=alpha_serie(100, 400, 12, 0, prt=False)
py.load_alpha(alpha_fit, power_a)
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: +8.3868e-05 +0.0 -1.0987e-01 -2.0 +6.9188e-03 -1.0 -1.2732e-03 -0.5
Or the function alpha_dir_serie implementing the method (3):
alpha_fit=alpha_dir_serie(100, 400, 12,0)
py.load_alpha(alpha_fit, power_a)
py.info()
Temp. Alpha 100.0 1.47e-05 127.3 1.85e-05 154.5 2.16e-05 181.8 2.43e-05 209.1 2.65e-05 236.4 2.85e-05 263.6 3.03e-05 290.9 3.20e-05 318.2 3.35e-05 345.5 3.49e-05 372.7 3.63e-05 400.0 3.76e-05 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: +9.1652e-05 +0.0 -1.5809e-01 -2.0 +9.0435e-03 -1.0 -1.5157e-03 -0.5
Or yet the volume_from_F_serie implementing the method (2):
alpha_fit=volume_from_F_serie(100, 400, 12, expansion=True, fit_alpha=True, export_alpha_fit=True)
py.load_alpha(alpha_fit, power_a)
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: +9.0668e-05 +0.0 -1.3181e-01 -2.0 +8.4162e-03 -1.0 -1.4671e-03 -0.5
In summary, the four methods to compute thermal expansion are:
method (1): thermal expansion from the log(V(T)) derivative by using the volume_dir function; this method is implemented in the alpha_dir_v (with argument comp=False)
method (2): thermal expansion from the log(V(T)) derivative by using the volume_from_F function; this method is implemented in the volume_from_F_serie function
method (3): thermal expansion from the local derivative of V(T) with respect to T; this method is implemented in the alpha_dir_serie function
method (4): thermal expansion from $(dP/dT)_V$ and EoS; this method is implemented in the alpha_serie function
All of the four methods can also be applied with the single alpha_dir_v function by setting to True the comp parameter.