Anharmonicity¶

Computation of energy levels of the anharmonic oscillator and of its contribution to the vibrational Helmholtz function¶

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:

  • file anharm_path.txt: it must be present in the root folder (the one containing the program). The content of anharm_path.txt is the name of the folder containing the data (usually, the folder relative to the phase to be investigated). Such name is assigned to the abs_path variable;
  • file input_anharm.txt: this file must be present in the folder specified in anharm_path.txt. The content of this file is:
    • folder name where SCAN data from CRYSTAL are stored (path variable);
    • output file name (it will be stored in the folder specified in the file anharm_path.txt);
    • minimim, maximum temperatures and number of points where the anharmonic Helmholtz function will be computed;
    • order of the power series used to fit the Helmholtz free energy as a function of V and T. The unit of the computed free energy is the hartree.

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):

  • volumes.dat: it contains the volumes at which the SCANMODE's where done together with the harmonic frequencies computed by CRYSTAL. If not both 0., the last two columns, specifies the minimum and maximum q's to select. Volumes of the primitive cell in cubic $\mathring{A}$; frequencies in $cm^{-1}$;
  • vect.dat: eigenvectors of the normal mode: one column for each volume, in the same order as specified in the volumes.dat file;
  • input.txt: it contains the names of the files where the Q-energies from the SCANMODE's are stored, as they are copied and pasted from the CRYSTAL output;
  • files whose names are written in the input.txt file.

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 folder quarzo_2020/scan_a contains all the required input files;
  • the output file will be quarzo_2020/anharm_a.dat;
  • the Helmholtz function will be fitted on the 24 points in the $[100, 1100\,K]$ temperature range;
  • the power of the V, T power serie is 4.

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:

  1. q value of the coordinate along the normal mode (in unit of $q_{max}$);
  2. actual energy (hartree) computed by CRYSTAL at the specified q value;
  3. energy (in hartree) computed in the limit of the harmonic approximation, from the harmonic $\omega$ frequency computed by CRYSTAL

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:

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

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

  • scan_name: name of the scan file
  • vol: unit cell volume
  • h_freq: harmonic frequency
  • red: reduced mass
  • ..

For instance:

In [3]:
anh[0].vol
Out[3]:
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:

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

In [5]:
anh[0].par
Out[5]:
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:

In [6]:
energy_anh(0)

The first 5 eigenvalues (energy in J) of such H_matrix are:

In [7]:
anh[0].vals[0:5]
Out[7]:
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:

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

In [9]:
anh[0].qmax
Out[9]:
0.1785349265550021

this means that the corresponding q values, in the unit used in the scan data file, are:

In [10]:
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).

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

In [12]:
v_index=5
print("Volume: %4.3f (A3)" % anh[v_index].vol)
Volume: 124.910 (A3)
In [13]:
start_fit(v_index)
In [14]:
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

  • compute the partition function and the Helmholtz free energy ($F$) values on the V, T grid defined by the input list of volumes and the temperature list, also specified in input
  • fit the power series, up to the specified degree, to the numerically computed $F(V,T)$ values

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.

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

How to make use of the anharmonic correction in the main program for the computation of the EoS, and other thermodynamic properties¶

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,

  • the mode number;
  • the Brillouin flag (BF): if it is an optical mode (BF=0) or an off-center vibrational mode (BF=1; computed by the supercell approach);
  • the degeneracy of the 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

File structure¶

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.