Read ab initio data from ALAMODE output files

Posted on In , Word count in article: 1.8k

In order to perform a fully ab initio calculation, for example, the VRMC simulation that we did in the last post, ab initio data, such as specific heat (spectral), phonon dispersion, and phonon lifetime, is needed. ALAMODE is a powerful tool to calculate these data, but the official documentation is not very clear about the formats...

Basically, all we need is in the *.result file, which contains basic information about the system, phonon frequencies, branches, and also linewidths. But you will also need the *.self_isotope file, which contains the self-energy due to the isotope scattering (to obtain that file, please add self_isotope = 2 rather than 1 in the anphon input file).

Example of these files (right click and save as):

Contents of output files

The *.result file

General information

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
## General information
#SYSTEM
2 1
265.904
#END SYSTEM
#KPOINT
36 36 36
1240
     1:   0.000000e+00   0.000000e+00   0.000000e+00    0.000021
     2:   0.000000e+00   0.000000e+00   2.777778e-02    0.000171
     ...
  1239:   2.222222e-01  -5.000000e-01  -2.777778e-01    0.000257
  1240:   2.500000e-01  -5.000000e-01  -2.500000e-01    0.000129
#END KPOINT
#CLASSICAL
0
#END CLASSICAL
#FCSXML
si333.xml
#END  FCSXML
#SMEARING
-1
10
#END SMEARING
#TEMPERATURE
200 400 20
#END TEMPERATURE
##END General information

As shown above, the first section is the general information about the system. Each subsection is separated by #END and #. The followings are the information we need:

  • #SYSTEM: the second line (265.904) in this subsection is the volume of the primitive cell (in Bohr\(^3\)).
  • #KPOINT: the first line (36 36 36) is the number of \(k\)-points in each direction, and their product is thus the total number of \(k\)-points. The second line (1240) is the irreducible number of \(k\)-points. The following lines are thus the coordinates of those irreducible \(k\)-points and their weights.
  • #TEMPERATURE: The line (200 400 20) is the temperature range and step size set by the anphon input. For this case, their will be 11 temperatures, and the result file will provide linewidths of all irreducible modes at each temperature.

Phonon frequency

1
2
3
4
5
6
7
8
##Phonon Frequency
#K-point (irreducible), Branch, Omega (cm^-1)
     1     1    1.09737e-10
     1     2    1.09737e-10
     ...
  1240     5        461.494
  1240     6        461.494
##END Phonon Frequency

This section contains frequencies of all irreducible modes. For each line, the first number is the irreducible \(k\)-point index (from 1 to 1240), the second number is the branch index (from 1 to 6), and the third number is the frequency. Note that the unit is cm\(^{-1}\), and you need to convert it to Rad/s by multiplying \(\frac{2\pi}{c}\times 1\times 10^{2}\), where \(c\) is the speed of light in m/s.

Phonon relaxation time

1
2
3
4
5
6
7
8
9
10
11
##Phonon Relaxation Time
#GAMMA_EACH
1 1
1
              0              0              0
              0
              0
              ...
       0.674943
       0.707127
#END GAMMA_EACH

This section contains group velocities and linewidths of all modes at each temperature, and it is organized by all irreducible modes.

In this section, there are \(n_{k, \text{irr}}\times n_s\) subsections (here is \(1240\times 6\)). Each subsection starts with #GAMMA_EACH and ends with #END GAMMA_EACH. Let's take a look at one of them:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#GAMMA_EACH
13 5
8
       -106.711       -106.711        105.928
       -106.711        105.928       -106.711
        105.928        105.928        105.928
        106.711        106.711       -105.928
       -105.928       -105.928       -105.928
        106.711       -105.928        106.711
       -105.928        106.711        106.711
        105.928       -106.711       -106.711
        1.03999
        1.11157
        1.18524
         1.2606
        1.33733
         1.4152
        1.49402
        1.57363
        1.65393
         1.7348
        1.81616
#END GAMMA_EACH

Here, the first line (13 5) represents the irreducible \(k\)-point index 13 and the branch index 5. The second line (8) is the number equivalent modes of this irreducible mode. In the following 8 lines, there are 3 numbers in each line, which are the three components of the group velocity of that equivalent mode. After that, there are 11 lines, which are the linewidths of that equivalent mode at each temperature.

The *.self_isotope file

The isotope scattering is not included in the *.result file, but it is very important for the calculation of the thermal conductivity. In the *.self_isotope file, there are the self-energy due to the isotope scattering of all irreducible modes, in the following format:

1
2
3
4
5
6
7
8
# Irreducible k point  :       48 (  24)
## xk =               0      0.0277778      -0.194444
     48    1        84.2973     0.00050785
     48    2        87.3591    0.000612916
     48    3        190.352      0.0128404
     48    4        484.132       0.411506
     48    5         492.31      0.0912846
     48    6        493.993      0.0783658

in which the first line is the irreducible \(k\)-point index and the number of equivalent modes, the second line is the coordinates of that irreducible \(k\)-point. After that, there are \(n_s\) lines, each of which contains the irreducible \(k\)-point index, the branch index, the frequency, and the isotope scattering linewidth.

Read data with Python

Use the class alamode_material to read the data:

alamode_material.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
"""
The class alamode_material is used to read the output files of ALAMODE
to facilitate ab initio calculations.

Author: :wq

License: CC BY-NC-SA 4.0 DEED (https://creativecommons.org/licenses/by-nc-sa/4.0/deed.en)
"""

import numpy as np
from scipy import constants as const

class alamode_material():
    """
    Example of use:
        >>> material = alamode_material('si333.result')
    Or with isotope scattering:
        >>> material = alamode_material('si333.result', 'si333.self_isotope')
    Attributes (all in SI units):
        self.name:              str, name of the material
        self.V_cell:            float, volume of the unit cell
        self.k_mesh:            np.array[3], mesh of k points
        self.nk:                int, total number of k points
        self.nk_irr:            int, number of irreducible k points
        self.ns:                int, number of phonon branches
        self.k_irr:             np.array[nk_irr], irreducible k points
        self.weight_irr:        np.array[nk_irr], weight of irreducible k points
        self.T_array:           np.array[:], temperature array
        self.omega:             np.array[nk_irr, ns], phonon frequency
        self.vg:                np.array[nk, ns, 3], group velocity
        self.full2irr:          np.array[nk], index of irreducible k points
        self.equiv_list:        list[nk_irr] of list, equivalent k points
        self.gamma:             np.array[nk_irr, ns, len(T_array)], self-energy (including isotope scattering)
        self.gamma_isotope:     np.array[nk_irr, ns], self-energy from isotope scattering
        self.scatt_time:        np.array[nk_irr, ns, len(T_array)], scattering time
        self.cv:                np.array[nk_irr, ns, len(T_array)], spectral specific heat
        self.cv_total:          np.array[len(T_array)], total volumetric specific heat
        self.conductivity_xx:   np.array[len(T_array)], thermal conductivity in xx direction
        self.conductivity_yy:   np.array[len(T_array)], thermal conductivity in yy direction
        self.conductivity_zz:   np.array[len(T_array)], thermal conductivity in zz direction
    """
    def __init__(self, name='My Material', result_path=None, isotope_path=None):
        self.name = name
        self.read_result(result_path)

        if result_path is None:
            raise ValueError('Please provide the path of the result file.')
        
        if isotope_path is not None:
            self.read_isotope(isotope_path)
            for i in range(len(self.T_array)):
                self.gamma[:, :, i] += self.gamma_isotope
        self.scatt_time = np.divide(1, 2*self.gamma, out=np.zeros_like(self.gamma), where=self.gamma!=0)

        self.calculate_cv()
        self.calculate_conductivity()

        return
    
    def read_result(self, result_path):
        with open(result_path, 'r') as f:
            lines = f.readlines()
        kayer_to_radpersecond = 2 * const.pi * const.c * 1e2

        line_number = 0
        while line_number < len(lines):
            if '#SYSTEM' in lines[line_number]:
                line_number+=2
                self.V_cell = float(lines[line_number]) * const.value('Bohr radius')**3
            if '#KPOINT' in lines[line_number]:
                line_number += 1
                self.k_mesh = np.array(lines[line_number].split(), dtype=int)
                self.nk = np.prod(self.k_mesh)

                line_number += 1
                self.nk_irr = int(lines[line_number])

                self.k_irr = np.zeros((self.nk_irr, 3))
                self.weight_irr = np.zeros(self.nk_irr)

                for i in range(self.nk_irr):
                    line_number += 1
                    self.k_irr[i, 0], self.k_irr[i, 1], self.k_irr[i, 2], self.weight_irr[i] = lines[line_number].split()[1:]

            if '#TEMPERATURE' in lines[line_number]:
                line_number += 1

                T_start = float(lines[line_number].split()[0])
                T_end = float(lines[line_number].split()[1])
                T_step = float(lines[line_number].split()[2])
                self.T_array = np.arange(T_start, T_end + 0.5* T_step, T_step)
                T_number = len(self.T_array)

            if '##Phonon Frequency' in lines[line_number]:
                line_number += 2
                self.ns = 0
                while int(lines[line_number + self.ns].split()[0]) == 1:
                    self.ns += 1

                # Allocate all parameters we need
                self.omega = np.zeros((self.nk_irr, self.ns))
                self.vg = np.zeros((self.nk, self.ns, 3))
                vg_read_index = np.zeros(self.ns, dtype=int)
                self.full2irr = np.zeros(self.nk, dtype=int)
                self.equiv_list = [[] for _ in range(self.nk_irr)]
                self.gamma = np.zeros((self.nk_irr, self.ns, T_number))

                for _ in range(self.nk_irr*self.ns):
                    k_i = int(lines[line_number].split()[0])
                    s_i = int(lines[line_number].split()[1])
                    self.omega[k_i-1, s_i-1] = float(lines[line_number].split()[2]) * kayer_to_radpersecond
                    line_number += 1

            if '#GAMMA_EACH' in lines[line_number]:
                line_number += 1
                k_i = int(lines[line_number].split()[0])
                s_i = int(lines[line_number].split()[1])
                line_number += 1
                equiv_number = int(lines[line_number])
                for _ in range(equiv_number):
                    vg_index = vg_read_index[s_i-1]
                    if s_i == 1:
                        self.full2irr[vg_index] = k_i-1
                        self.equiv_list[k_i-1].append(vg_index)
                    
                    line_number += 1
                    self.vg[vg_index, s_i-1, 0], self.vg[vg_index, s_i-1, 1], self.vg[vg_index, s_i-1, 2] = lines[line_number].split()
                    vg_read_index[s_i-1] += 1
                for i in range(T_number):
                    line_number += 1
                    self.gamma[k_i-1, s_i-1, i] = float(lines[line_number]) * kayer_to_radpersecond

            line_number += 1
        
        return
    
    def read_isotope(self, isotope_path):
        self.gamma_isotope = np.zeros((self.nk_irr, self.ns))

        with open(isotope_path, 'r') as f:
            lines = f.readlines()
        kayer_to_radpersecond = 2 * const.pi * const.c * 1e2

        line_number = 0
        while line_number < len(lines):
            if '# Irreducible k point  :' in lines[line_number]:
                k_i = int(lines[line_number].replace('(',' ').replace(')',' ').split()[-2])
                dos = int(lines[line_number].replace('(',' ').replace(')',' ').split()[-1])
                if dos != len(self.equiv_list[k_i-1]):
                    raise ValueError('DoS in self_isotope file is not consistent with the result file.')
                else:
                    line_number += 2
                    for i in range(self.ns):
                        if int(lines[line_number].split()[0]) != k_i or int(lines[line_number].split()[1]) != i+1:
                            raise ValueError('Error in reading self_isotope file at line %d.'% line_number)
                        self.gamma_isotope[k_i-1, i] = float(lines[line_number].split()[-1]) * kayer_to_radpersecond
                        line_number += 1
                
            line_number += 1
        
        return
    
    def calculate_cv(self):
        self.cv = np.zeros((self.nk_irr, self.ns, len(self.T_array)))
        for i in range(len(self.T_array)):
            if self.T_array[i] < 1e-12:
                self.cv[:, :, i] = 0
                continue                
            x = const.hbar * self.omega / (const.k * self.T_array[i])
            self.cv[:, :, i] = const.k * (x / (2 * np.sinh(x/2)))**2 / self.V_cell
            
            self.cv_total = 0
            for i in range(self.nk_irr):
                self.cv_total += self.cv[i, :, :].sum(axis=0) * len(self.equiv_list[i]) / self.nk
        
        return

    def calculate_conductivity(self):
        self.conductivity_xx = np.zeros(len(self.T_array))
        self.conductivity_yy = np.zeros(len(self.T_array))
        self.conductivity_zz = np.zeros(len(self.T_array))
        for k_i in range(self.nk):
            for s_i in range(self.ns):
                self.conductivity_xx += 1/self.nk * self.cv[self.full2irr[k_i] ,s_i, :] * self.vg[k_i, s_i, 0]**2 * self.scatt_time[self.full2irr[k_i], s_i, :]
                self.conductivity_yy += 1/self.nk * self.cv[self.full2irr[k_i] ,s_i, :] * self.vg[k_i, s_i, 1]**2 * self.scatt_time[self.full2irr[k_i], s_i, :]
                self.conductivity_zz += 1/self.nk * self.cv[self.full2irr[k_i] ,s_i, :] * self.vg[k_i, s_i, 2]**2 * self.scatt_time[self.full2irr[k_i], s_i, :]
        
        return

Example of use:

1
2
3
4
5
6
7
8
from alamode_material import alamode_material
import numpy as np

Si = alamode_material('Si', 'si333.result', 'si333.self_isotope')

# Verify the specific heat and thermal conductivity values at 300K
print("Volumetric specific heat at 300 K: {:.2f} J/(m^3 K)".format(Si.cv_total[np.where(Si.T_array == 300)[0][0]]))
print("Thermal conductivity in xx direction at 300 K: {:.2f} W/(m K)".format(Si.conductivity_xx[np.where(Si.T_array == 300)[0][0]]))

Getting the following output:

1
2
Volumetric specific heat at 300 K: 1681023.27 J/(m^3 K)
Thermal conductivity in xx direction at 300 K: 143.18 W/(m K)

Note that the calculated thermal conductivity is slightly different from the value given by ALAMODE in *.kl file. I suppose that this is due to numerical errors. Luckily, the difference is very small, for Si, it is about 0.3 W/(m K) at 300 K.