The program supports the calculation of the effects of phonon dispersion on thermodynamic properties. The required input data can be obtained by any program like CRYSTAL17 implementing a supercell approach to get vibrational frequencies of the off-center modes.
The example below refers to the case of the $\alpha$-quartz, for which phonon frequencies were computed for a 2x2x2 supercell: in addition to the frequencies at the center of the Brillouin zone (BZ), 7 other sets of frequencies were obtained corresponding to the (0 0 1/2), (0 1/2 0), (1/2 0 0), (0 1/2 1/2), (1/2 0 1/2), (1/2 1/2 0) and (1/2 1/2 1/2) K-points of the BZ (the DISPERSI keyword was specified in the CRYSTAL input file). Supercell frequencies were computed for several values of the unit cell volume.
In order to require phonon dispersion correction, the keyword DISP must be present in the input file input.txt; such keyword must be followed by the names of two files containing the relevant data. In the present case:
DISP
quartz_disp.dat
quartz_disp_info.dat
The file quartz_disp.dat contains the set of frequencies, each column referring to a specific unit cell volume; the first column of the file does contain the degeneracy of each mode. Along each column, the frequencies are ordered according to the following hierarchy:
the frequencies of the modes belonging to each one of the K-points of the BZ, that is, the (0 0 1/2) modes first, then the (0 1/2 0) modes... and, lastly, the (1/2 1/2 1/2) modes
the symmetry of the modes
In this way, all the frequencies along a given row of the input file refer to the same mode, as it changes as a function of volume. The example below is a section of the quartz_disp.dat, where the frequencies of the (0 0 1/2) K-point are shown, calculated for three different values of the unit cell volume:
1 209.578 211.24 210.95
1 401.8457 387.08 395.83
1 790.0472 776.39 783.82
1 1120.1813 1128.20 1124.21
1 56.7594 47.51 53.18
1 171.7892 181.46 175.85
1 393.5102 372.92 385.14
1 779.9419 762.05 772.48
1 1061.4419 1055.16 1059.07
2 130.2962 105.66 121.46
2 188.5704 193.30 190.35
2 263.4401 247.93 257.05
2 428.0019 422.10 426.00
2 458.9878 452.47 455.91
2 580.3685 563.36 574.12
2 794.6851 776.13 786.05
2 1080.7245 1074.68 1078.27
2 1222.3513 1222.03 1222.12
The file quartz_disp_info.dat contains information about the structure of the quartz_disp.dat file, which are required by the program to correctly interpret the data. The file in this example is:
7 3 5. 1200.
18 27 27 27 27 27 27
1 1 1 1 1 1 1
114.860943 121.000435 117.551141 ...
We have four rows whose content is:
Notes:
Having prepared all of the required input, the program is here started:
%matplotlib inline
%run bm3_thermal_2
This is BMx-QM program, version 2.4.0 - 18/02/2021 Run time: 2021-05-14 17:18:57.783259 File quick_start.txt found in the master folder Input files in 'quarzo_2020' folder Input file TITLE Quartz 2020 STATIC static.dat VOLUME volumes.dat FREQ frequencies.dat FITVOL POLY 92. 126. 16 3 SET 0 1 2 3 4 5 6 7 8 9 10 11 FU 3 3 LO LO.txt DISP quartz_disp.dat quartz_disp_info.dat END -------- End of input file ------- Phonon dispersion correction activated by the DISP keyword in the input file; the contribution to the entropy and to the specific heat is taken into account. Warning ** Volume and frequencies lists have been sorted indexing: [ 0 1 2 3 11 8 10 4 9 5 6 7] Frequencies corrected for LO-TO splitting Static BM3 EoS Bulk Modulus: 43.32 (0.25) GPa Kp: 4.73 (0.09) V0: 114.7420 (0.02) A^3 E0: -1.31891502e+03 (1.52e-05) hartree Instructions from quick_start file will be executed: disp.free_fit_ctrl(degree=6) disp.eos_on() volume_ctrl.set_all(degree=2, upgrade_shift=True, quad_shrink=8) -------------------------------------------------- disp.free_fit_ctrl(degree=6) V,T-fit of the Helmholtz free energy contribution from the off-centered modes Power of the fit: 6 Mean error: 1.06e-04 Maximum error: 3.59e-04
disp.eos_on() Phonon dispersion correction for bulk_dir computation volume_ctrl.set_all(degree=2, upgrade_shift=True, quad_shrink=8)
<Figure size 480x320 with 0 Axes>
Note the following commands:
disp.free_fit_ctrl(degree=6)
This command invokes a method of the disp class: free_fit_ctrl that, in the present case, sets to 6 the degree of the F(V,T) power series.
disp.eos_on()
It requests the computation of F(V,T) function (power series up to the 6-th degree as specified above).
volume_ctrl.set_all(degree=2, upgrade_shift=True, quad_shrink=8)
It invokes a volume driver that sets appropriate numerical parameters for the calculation of the volume at given temperature and pressure, (un)compressibilities and thermal expansion, in a direct way, without the recourse to any equation of state.
For more, see the help on the volume_control_class in the index section of this site.
As an example, we calculate here the bulk modulus $K$ as a function of temperature in the range $[50, 700K]$ (24 T values in the range), at the zero pressure. The calculation is performed according to the definition
$$K = -V\left(\frac{\partial P}{\partial V}\right)_T$$where P(V), at any given T, is directly computed as
$$P=-\left(\frac{\partial F}{\partial V}\right)_T$$To the purpose, the function bulk_modulus_p_serie is used
bulk_modulus_p_serie(50, 700, 24, 0, noeos=True, fit=True, type='spline', deg=3, smooth=2)
To compute the bulk modulus at a given temperature, from an EoS fit (3-th order Birch-Murnaghan EoS), the function bulk_dir can be used:
bulk_dir(300)
BM3 EoS from P(V) fit K0: 40.37 (0.14) GPa Kp: 5.00 (0.02) V0: 116.2296 (0.02) A^3
Above, the computation was done with the phonon dispersion effects included. We can now switch off the phonon dispersion correction and recompute the fit:
disp.eos_off()
bulk_dir(300)
No phonon dispersion correction for bulk_dir computation BM3 EoS from P(V) fit K0: 38.79 (0.23) GPa Kp: 5.18 (0.04) V0: 116.4281 (0.03) A^3
In this case, the phonon dispersion had a relatively small impact on the computed bulk modulus $K_0$ at T=300 K.
Entropy and specific heat are particularly sensitive to phonon dispersion effects. With dispersion included:
entropy_p(298,0)
cp(298,0,prt=True)
Entropy: 39.63 J/mol K Specific heat (at constant volume): 41.80 J/mol K Pressure: 0.00 GPa; Volume 116.6736 A^3 Cp: 42.22, Cv: 41.80, Cp/Cv: 1.01016, alpha: 4.080e-05, K: 36.57
With no dispersion:
disp.off()
entropy_p(298,0)
cp(298,0, prt=True)
Dispersion correction off Entropy: 27.11 J/mol K Specific heat (at constant volume): 35.50 J/mol K Pressure: 0.00 GPa; Volume 116.6736 A^3 Cp: 35.92, Cv: 35.50, Cp/Cv: 1.01197, alpha: 4.080e-05, K: 36.57
The computation with dispersion included agrees with the experimental data by Richet et al. (1982)