With reference to the soft-mode driving the high temperature $\alpha-\beta$ transition in quartz, which is intrinsecally strongly anharmonic, the use of the program anharm will be here discussed.
The input required to run the program is summarized in the following:
The output file contains the powers of the fitting polynomial together with the optimized coefficients for the reconstruction of the Helmholtz free energy as a function of V and T in the specified ranges. Volume ranges are from the input files found in the specified folder (path variable).
Files required to be found in the specified folder (path variable):
In the present case, the file anharm_path.txt file contains the folder name quarzo_2020.
The content of the file input_anharm.txt (in the quarzo_2020 folder) is:
scan_a
anharm_a.dat
100 1100 24
4
This means that:
The volumes.dat file content (corresponding to 7 different data-sets, at different unit cell volumes, from the CRYSTAL/SCAN computations) is:
114.860943 215.6839 -8. 8.
109.884957 248.0118 -8. 8.
104.662515 274.6727 -8. 8.
97.900625 296.8036 -8. 8.
121.000240 161.2544 -6. 15.
124.910088 104.4002 -5.5 11.5
122.901896 138.5263 -5.5 13.5
The input.txt file content is:
scan_1.dat
scan_2.dat
scan_3.dat
scan_4.dat
scan_5.dat
scan_6.dat
scan_7.dat
The structure of each scan_i.dat file (in the example is the scan_1.dat file):
-6.0000 -0.1318890832359E+04 -0.1318897402955E+04
-5.6000 -0.1318894314601E+04 -0.1318899682885E+04
.
.
11.6000 -0.1318885734175E+04 -0.1318848974101E+04
12.0000 -0.1318884517817E+04 -0.1318844335623E+04
The columns of such file (copied and pasted from the CRYSTAL mail output file) are:
The content of the vect.dat file is:
0.0319 0.0348 0.0393 0.0461 0.0305 0.0338 0.0313
-0.0184 -0.0201 -0.0227 -0.0266 -0.0176 -0.0195 -0.0181
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0.0368 0.0401 0.0454 0.0532 0.0353 0.0390 0.0362
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
-0.0319 -0.0348 -0.0393 -0.0461 -0.0305 -0.0338 -0.0313
-0.0184 -0.0201 -0.0227 -0.0266 -0.0176 -0.0195 -0.0181
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
-0.0182 -0.0191 -0.0177 -0.0107 -0.0154 -0.0125 -0.0142
-0.0527 -0.0483 -0.0446 -0.0398 -0.0618 -0.0772 -0.0669
0.0391 0.0309 0.0215 0.0069 0.0514 0.0677 0.0569
0.0547 0.0513 0.0475 0.0399 0.0612 0.0731 0.0650
0.0106 0.0076 0.0070 0.0106 0.0175 0.0278 0.0212
0.0391 0.0309 0.0215 0.0069 0.0514 0.0677 0.0569
-0.0365 -0.0323 -0.0298 -0.0291 -0.0458 -0.0606 -0.0508
0.0421 0.0407 0.0376 0.0292 0.0442 0.0494 0.0457
0.0391 0.0309 0.0215 0.0069 0.0514 0.0677 0.0569
0.0365 0.0323 0.0298 0.0291 0.0458 0.0606 0.0508
0.0421 0.0407 0.0376 0.0292 0.0442 0.0494 0.0457
-0.0391 -0.0309 -0.0215 -0.0069 -0.0514 -0.0677 -0.0569
0.0182 0.0191 0.0177 0.0107 0.0154 0.0125 0.0142
-0.0527 -0.0483 -0.0446 -0.0398 -0.0618 -0.0772 -0.0669
-0.0391 -0.0309 -0.0215 -0.0069 -0.0514 -0.0677 -0.0569
-0.0547 -0.0513 -0.0475 -0.0399 -0.0612 -0.0731 -0.0650
0.0106 0.0076 0.0070 0.0106 0.0175 0.0278 0.0212
-0.0391 -0.0309 -0.0215 -0.0069 -0.0514 -0.0677 -0.0569
Each column of the file (from the CRYSTAL main output file) represents the components of eigenvector (normal mode) considered, at a given unit cell volume. It is used to compute the value of the reduced mass of the mode, at each volume.
Let us load and run the program anharm.py:
%matplotlib inline
%run anharm.py
Run time: 2021-05-07 11:22:49.523691 This is Anharm, version 1.1 - 16/07/2020 Basic usage: helm_fit() Number of data SCAN's: 7: Volumes: [114.860943 109.884957 104.662515 97.900625 121.00024 124.910088 122.901896]
The program will immediately load all the input files except the scan data file. A list of the unit cell volumes considered is printed. The function helm_fit() will run the program and will produce the desired output (file anharm_a.dat); however one can also proceed step by step, as it will be shown in the following. At the step 1, the set_up() function is used:
set_up()
set_up read all of the scan data files and saves all the relevant quantities in the array anh. Each element of such an array (the anh array is an instance of a class: the defined anh_class), one for each volume, has several attributes, for instance:
For instance:
anh[0].vol
114.860943
is the value of the unit cell volume (in $\mathring{A}^{3}$) corresponding to the first data set.
The start_fit function makes the quartic and the quadratic fit of the $E(q)$ curve numerically obtained by the CRYSTAL SCANMODE computation. The (mandatory) argument to the start_fit function is the index referring to the volume that has to be considered:
start_fit(0)
Indeed, this is the fit corresponding to the volume stored in the variable anh[0].vol. The parameters resulting from such fit are stored in the anh[0].par variable:
anh[0].par
array([ 2.17177972e-05, 1.52464168e-02, -5.10495380e-03, 4.45067540e-04])
Such parameters are then used to compute the matrix representation of the anharmonic Hamiltonian in the space of the eigenfunction of the harmonic oscillator (having anh[0].h_freq frequency). The H matrix is then diagonalized and its eigenvalues (sorted by value) are stored in the variable anh[0].vals. All this is performed by the energy_anh function which, once again, requires the volume index as mandatory argument:
energy_anh(0)
The first 5 eigenvalues (energy in J) of such H_matrix are:
anh[0].vals[0:5]
array([2.12667613e-21, 6.36466990e-21, 1.05794276e-20, 1.47705944e-20, 1.89378012e-20])
The energy values ($\epsilon_i$) corresponding to the vibrational levels of the anharmonic oscillator are then used to numerically compute the partition function at a given temperature T:
$$ Z = \sum_{i=0}^N e^{-\epsilon_i/kT} $$The number ($N$) of levels to be included depends upon the temperature and a check for the convergence of the $Z$ function against $N$ must be performed. The function check_partition is used for the purpose. The arguments of the function are the volume index and the temperature; for instance, for the volume we are considering, at T=1000 K:
check_partition(0,1000)
Threshold for convergence (on the variation of Z): 0.0050 Convergence reached at the 19 level Minimum and maximum q values: -0.93, 1.36
That is: convercence is reached with 19 energy levels included in the calculation; the minimum and maximum q values (in unit of $q_{max}$) are -0.93, 1.36. As qmax is:
anh[0].qmax
0.1785349265550021
this means that the corresponding q values, in the unit used in the scan data file, are:
qmn=-0.93/anh[0].qmax
qmx=1.36/anh[0].qmax
print(qmn.round(2), qmx.round(2))
-5.21 7.62
To visualize the vibrational level structure (up the 19-th level), over the $E(q)$ curve plotted in the [qmn, qmx] range, the function plot_levels can be used. The (mandatory) arguments to the function are the volume index and the number of levels to be plotted. Optional arguments are qmin and qmax (if not specified, the whole q interval will be used).
plot_levels(0, 19, qmin=qmn, qmax=qmx)
Energy levels occupation (probabilities) at several temperatures in the 300.0 - 1000.0 interval 300.0 475.0 650.0 825.0 1000.0 63.9 47.4 37.4 30.7 26.1 23.0 24.8 23.3 21.2 19.2 8.3 13.1 14.6 14.6 14.1 3.0 6.9 9.1 10.1 10.4 1.1 3.7 5.7 7.0 7.7 0.4 1.9 3.6 4.9 5.7 0.2 1.0 2.3 3.4 4.2 0.1 0.6 1.4 2.4 3.2 0.0 0.3 0.9 1.7 2.3 0.0 0.2 0.6 1.2 1.8 0.0 0.1 0.4 0.8 1.3 0.0 0.0 0.2 0.6 1.0 0.0 0.0 0.2 0.4 0.7 0.0 0.0 0.1 0.3 0.6 0.0 0.0 0.1 0.2 0.4 0.0 0.0 0.0 0.1 0.3 0.0 0.0 0.0 0.1 0.2 0.0 0.0 0.0 0.1 0.2 0.0 0.0 0.0 0.1 0.1
A table with occupation values (probability of occupation) of each vibrational level against the temperature is also printed.
At an higher cell volume, more levels are required to reach convergence:
v_index=5
print("Volume: %4.3f (A3)" % anh[v_index].vol)
Volume: 124.910 (A3)
start_fit(v_index)
energy_anh(v_index)
check_partition(v_index, 900)
Threshold for convergence (on the variation of Z): 0.0050 Convergence reached at the 34 level Minimum and maximum q values: -1.26, 2.73
In order to
the function helm_fit can be used. The (optional) argument given is the temperature at which the convergence of the $Z$ function is checked. helm_fit will perform all the required $E(q)$ fits; computes the vibrational energies as eigenvalues of the H matrix; computes the $F(V,T)$ values and eventually makes the $F$ fit.
helm_fit(temp=900)
The file anharm_a.dat does contains the parameters of such a fit:
0 0 -1.3161037711227497
0 1 0.0002901312706983102
1 0 0.04826796015579944
0 2 8.331363275061264e-09
1 1 -8.607400365677583e-06
2 0 -0.0006613495408658537
0 3 2.073483891555499e-12
1 2 -2.465509273634077e-10
2 1 8.522339541870452e-08
3 0 4.014481368419415e-06
0 4 -1.6198649836458328e-15
1 3 3.242965795628521e-14
2 2 7.877944085208924e-13
3 1 -2.802389398877478e-10
4 0 -9.110525332100122e-09
END
This file can be used as input to the program computing the thermodynamics of the considered phase.
In order to use the anharm_a.dat file to perform the anharmonic correction in the main bm3_thermal program, you should insert the ANH keyword in the input file (input.txt) of the main program; the keyword must be followed by an integer number that specifies the number of modes corrected for anharmonicity and, for each mode,
Here is the example for quartz where three anharmonic modes were taken into account (modes n. 0, 8 and 4). The mode 4 is an off-center one (BF=1). Mode 8 is degenerate (its multiciplicity is equal to 2).
ANH
3
0 0 1
8 0 2
4 1 1
The file name anharm_a.dat must be written in the anh.txt file, which is read by the bm3_thermal_2 program; such file must be present in the folder containing the data of the considered phase.
In the main folder (the one containing the bm3_thermal_2 program), the file anharm_path.txt must be present and it should contain the name of the folder referring to the considered phase.
In the example of quartz, three anharm files are produced (e.g. anharm_a.dat, anharm_b.dat, anharm_c.dat), one for each anharmonic mode. Such three names must be written in the anh.txt in the quartz folder.