phono3py --fc3 --fc2 --dim="2 2 2" --loglevel=2 --mesh="64 64 64" -c POSCAR-unitcell --ga="0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8 9 9 9 10 10 10 11 11 11 12 12 12 13 13 13 14 14 14 15 15 15 16 16 16 17 17 17 18 18 18 19 19 19 20 20 20 21 21 21 22 22 22 23 23 23 24 24 24 25 25 25 26 26 26 27 27 27 28 28 28 29 29 29 30 30 30 31 31 31 32 32 32" --lw --bi="1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18"

  1 #!/usr/bin/env python
  2
  3 import os
  4 import sys
  5 import numpy as np
  6 import matplotlib.pyplot as plt
  7
  8 lw_dir = "linewidth"
  9 freq_file = "m646464.out"
 10 temp = 10 # temperature
 11 lw_factor = 10.0 # magnified factor
 12
 13 def readOutput(file_name):
 14     with open(file_name, 'rb') as f:
 15         all_lines = f.readlines()
 16
 17     all_found = 0
 18     for i in range(len(all_lines)):
 19         if all_found == 3:
 20             break
 21         elif "Mesh sampling" in all_lines[i]:
 22             mesh = map(int, all_lines[i].split()[3:6])
 23             all_found += 1
 24         elif "Band indices" in all_lines[i]:
 25             band_index = map(int, all_lines[i].replace('[','').replace(']','').split()[2:])
 26             all_found += 1
 27         elif "Grid point to" in all_lines[i]:
 28             gp_index = map(int, all_lines[i+1].split()[:])
 29             all_found += 1
 30
 31     return mesh, band_index, gp_index
 32
 33 def getLW(grid_index, temp, mesh, n_band):
 34     lw_grid = []
 35     mesh_file_name = ''.join(map(str, mesh))
 36     for i in range(n_band):
 37         file_name = os.path.join(lw_dir, 'linewidth-m'+mesh_file_name+'-g'+str(grid_index)+'-b'+str(i+1)+'.dat')
 38         with open(file_name, 'rb') as f:
 39             all_lines = f.readlines()
 40
 41         for line in all_lines:
 42             # add space and 0000 to avoid alias
 43             if (" "+str(float(temp))+"0000") in line:
 44                 lw_grid.append(float(line.split()[-1]))
 45
 46     if len(lw_grid) != n_band:
 47         print("Error!\n")
 48         sys.exit(1)
 49
 50     return lw_grid
 51
 52 def getFreq(grid_index, n_band):
 53     freq_grid = []
 54     with open(freq_file, 'rb') as f:
 55        all_lines = f.readlines()
 56
 57     for i in range(len(all_lines)):
 58         if " Grid point" in all_lines[i] and int(all_lines[i].split()[3]) == grid_index:
 59             # may have different lines for this data block
 60             for j in range(i+5, i+9):
 61                 freq_grid.extend(map(float, all_lines[j].strip().lstrip('[').rstrip(']').split()[:]))
 62             if ']' not in all_lines[i+8]:
 63                 freq_grid.extend(map(float, all_lines[i+9].strip().lstrip('[').rstrip(']').split()[:]))
 64
 65     if len(freq_grid) != n_band:
 66         print("Error!\n")
 67         sys.exit(1)
 68
 69     return freq_grid
 70
 71
 72 # Main routine
 73 mesh, band_index, gp_index = readOutput(freq_file)
 74
 75 freq = []
 76 lw = []
 77 n_band = len(band_index)
 78 x = range(len(gp_index))
 79
 80 for grid_index in gp_index:
 81     lw_tmp = getLW(grid_index, temp, mesh, n_band)
 82     lw.append(lw_tmp)
 83
 84     freq_tmp = getFreq(grid_index, n_band)
 85     freq.append(freq_tmp)
 86
 87 lw_upper = np.add(freq, 0.5*lw_factor*np.array(lw))
 88 lw_lower = np.add(freq, -0.5*lw_factor*np.array(lw))
 89
 90 freq_plot = zip(*freq)
 91 lw_plot = zip(*lw)
 92 lw_upper_plot = zip(*lw_upper)
 93 lw_lower_plot = zip(*lw_lower)
 94
 95 fig,ax = plt.subplots()
 96 for i in range(n_band):
 97     ax.plot(x, freq_plot[i], color='black')
 98     ax.plot(x, lw_upper_plot[i], x, lw_lower_plot[i], linestyle='--', linewidth=1, color='red')
 99     ax.fill_between(x, lw_upper_plot[i], lw_lower_plot[i], where=np.array(lw_upper_plot[i]) >= np.array(lw_lower_plot[i]), facecolor=None, interpolate=True)
100     #ax.fill_between(x, lw_upper_plot[i], lw_lower_plot[i], where=np.array(lw_upper_plot[i]) >= np.array(lw_lower_plot[i]), markerfacecolor="None", interpolate=True)
101 ax.set_title('ZrH2 phonon band structure with linewidth')
102 ax.set_xlabel('Kpoints')
103 ax.set_ylabel('Frequency(THz)')
104 plt.show()