/scripts/Inelastic/dos/load_phonon.py

https://github.com/wdzhou/mantid · Python · 185 lines · 110 code · 38 blank · 37 comment · 31 complexity · 19c2c900a8b6a2d1d377465d57c60ccc MD5 · raw file

  1. #pylint: disable=redefined-builtin
  2. from __future__ import (absolute_import, division, print_function)
  3. from six.moves import range
  4. import re
  5. import numpy as np
  6. import dos.load_helper as load_helper
  7. element_isotope = dict()
  8. def parse_phonon_file(file_name, record_eigenvectors):
  9. """
  10. Read frequencies from a <>.phonon file
  11. @param file_name - file path of the file to read
  12. @return the frequencies, infra red and raman intensities and weights of frequency blocks
  13. """
  14. file_data = {}
  15. # Get regex strings from load_helper
  16. header_regex = re.compile(load_helper.PHONON_HEADER_REGEX)
  17. eigenvectors_regex = re.compile(load_helper.PHONON_EIGENVEC_REGEX)
  18. block_count = 0
  19. frequencies, ir_intensities, raman_intensities, weights, q_vectors, eigenvectors = [], [], [], [], [], []
  20. data_lists = (frequencies, ir_intensities, raman_intensities)
  21. with open(file_name, 'rU') as f_handle:
  22. file_data.update(_parse_phonon_file_header(f_handle))
  23. while True:
  24. line = f_handle.readline()
  25. # Check we've reached the end of file
  26. if not line:
  27. break
  28. # Check if we've found a block of frequencies
  29. header_match = header_regex.match(line)
  30. if header_match:
  31. block_count += 1
  32. weight, q_vector = load_helper._parse_block_header(header_match, block_count)
  33. weights.append(weight)
  34. q_vectors.append(q_vector)
  35. # Parse block of frequencies
  36. for line_data in _parse_phonon_freq_block(f_handle,
  37. file_data['num_branches']):
  38. for data_list, item in zip(data_lists, line_data):
  39. data_list.append(item)
  40. vector_match = eigenvectors_regex.match(line)
  41. if vector_match:
  42. if record_eigenvectors:
  43. # Parse eigenvectors for partial dos
  44. vectors = _parse_phonon_eigenvectors(f_handle,
  45. file_data['num_ions'],
  46. file_data['num_branches'])
  47. eigenvectors.append(vectors)
  48. else:
  49. # Skip over eigenvectors
  50. for _ in range(file_data['num_ions'] * file_data['num_branches']):
  51. line = f_handle.readline()
  52. if not line:
  53. raise IOError("Bad file format. Uexpectedly reached end of file.")
  54. frequencies = np.asarray(frequencies)
  55. ir_intensities = np.asarray(ir_intensities)
  56. raman_intensities = np.asarray(raman_intensities)
  57. warray = np.repeat(weights, file_data['num_branches'])
  58. eigenvectors = np.asarray(eigenvectors)
  59. file_data.update({
  60. 'frequencies': frequencies,
  61. 'ir_intensities': ir_intensities,
  62. 'raman_intensities': raman_intensities,
  63. 'weights': warray,
  64. 'q_vectors':q_vectors,
  65. 'eigenvectors': eigenvectors
  66. })
  67. return file_data, element_isotope
  68. #----------------------------------------------------------------------------------------
  69. def _parse_phonon_file_header(f_handle):
  70. """
  71. Read information from the header of a <>.phonon file
  72. @param f_handle - handle to the file.
  73. @return List of ions in file as list of tuple of (ion, mode number)
  74. """
  75. file_data = {'ions': []}
  76. while True:
  77. line = f_handle.readline()
  78. if not line:
  79. raise IOError("Could not find any header information.")
  80. if 'Number of ions' in line:
  81. file_data['num_ions'] = int(line.strip().split()[-1])
  82. elif 'Number of branches' in line:
  83. file_data['num_branches'] = int(line.strip().split()[-1])
  84. elif 'Unit cell vectors' in line:
  85. file_data['unit_cell'] = _parse_phonon_unit_cell_vectors(f_handle)
  86. elif 'Fractional Co-ordinates' in line:
  87. if file_data['num_ions'] is None:
  88. raise IOError("Failed to parse file. Invalid file header.")
  89. # Extract the mode number for each of the ion in the data file
  90. for _ in range(file_data['num_ions']):
  91. line = f_handle.readline()
  92. line_data = line.strip().split()
  93. species = line_data[4]
  94. ion = {'species': species}
  95. ion['fract_coord'] = np.array([float(line_data[1]), float(line_data[2]), float(line_data[3])])
  96. ion['isotope_number'] = float(line_data[5])
  97. element_isotope[species] = float(line_data[5])
  98. # -1 to convert to zero based indexing
  99. ion['index'] = int(line_data[0]) - 1
  100. ion['bond_number'] = len([i for i in file_data['ions'] if i['species'] == species]) + 1
  101. file_data['ions'].append(ion)
  102. if 'END header' in line:
  103. if file_data['num_ions'] is None or file_data['num_branches'] is None:
  104. raise IOError("Failed to parse file. Invalid file header.")
  105. return file_data
  106. #----------------------------------------------------------------------------------------
  107. def _parse_phonon_freq_block(f_handle, num_branches):
  108. """
  109. Iterator to parse a block of frequencies from a .phonon file.
  110. @param f_handle - handle to the file.
  111. """
  112. for _ in range(num_branches):
  113. line = f_handle.readline()
  114. line_data = line.strip().split()[1:]
  115. line_data = [float(x) for x in line_data]
  116. yield line_data
  117. #----------------------------------------------------------------------------------------
  118. def _parse_phonon_unit_cell_vectors(f_handle):
  119. """
  120. Parses the unit cell vectors in a .phonon file.
  121. @param f_handle Handle to the file
  122. @return Numpy array of unit vectors
  123. """
  124. data = []
  125. for _ in range(3):
  126. line = f_handle.readline()
  127. line_data = line.strip().split()
  128. line_data = [float(x) for x in line_data]
  129. data.append(line_data)
  130. return np.array(data)
  131. #----------------------------------------------------------------------------------------
  132. def _parse_phonon_eigenvectors(f_handle, num_ions, num_branches):
  133. vectors = []
  134. for _ in range(num_ions * num_branches):
  135. line = f_handle.readline()
  136. if not line:
  137. raise IOError("Could not parse file. Invalid file format.")
  138. line_data = line.strip().split()
  139. vector_componets = line_data[2::2]
  140. vector_componets = [float(x) for x in vector_componets]
  141. vectors.append(vector_componets)
  142. return np.asarray(vectors)
  143. #----------------------------------------------------------------------------------------