Coverage for /builds/debichem-team/python-ase/ase/io/siesta_output.py: 97.76%

1import numpy as np 

2 

3from ase.units import Bohr, Ry 

4 

5 

6class OutputReader: 

7 def __init__(self, prefix, directory, bandpath=None): 

8 self.prefix = prefix 

9 self.directory = directory 

10 self.bandpath = bandpath 

11 

12 def read_results(self): 

13 results = {} 

14 results['n_grid_point'] = self.read_number_of_grid_points() 

15 results.update(self.read_energy()) 

16 results.update(self.read_forces_stress()) 

17 results.update(self.read_eigenvalues()) 

18 results.update(self.read_kpoints()) 

19 results['dipole'] = self.read_dipole() 

20 

21 if self.bandpath is not None and len(self.bandpath.kpts): 

22 results['bandstructure'] = self.read_bands(self.bandpath) 

23 

24 return results 

25 

26 def _prefixed(self, extension): 

27 return self.directory / f'{self.prefix}.{extension}' 

28 

29 def read_bands(self, bandpath): 

30 fname = self._prefixed('bands') 

31 with open(fname) as fd: 

32 kpts, energies, efermi = read_bands_file(fd) 

33 return resolve_band_structure(bandpath, kpts, energies, efermi) 

34 

35 def read_number_of_grid_points(self): 

36 """Read number of grid points from SIESTA's text-output file. """ 

37 

38 fname = self.directory / f'{self.prefix}.out' 

39 with open(fname) as fd: 

40 for line in fd: 

41 line = line.strip().lower() 

42 if line.startswith('initmesh: mesh ='): 

43 return [int(word) for word in line.split()[3:8:2]] 

44 

45 raise RuntimeError 

46 

47 def read_energy(self): 

48 """Read energy from SIESTA's text-output file. 

49 """ 

50 text = self._prefixed('out').read_text().lower() 

51 

52 assert 'final energy' in text 

53 lines = iter(text.split('\n')) 

54 

55 # Get the energy and free energy the last time it appears 

56 for line in lines: 

57 has_energy = line.startswith('siesta: etot =') 

58 if has_energy: 

59 energy = float(line.split()[-1]) 

60 line = next(lines) 

61 # XXX dangerous, this should test the string in question. 

62 free_energy = float(line.split()[-1]) 

63 

64 return {'energy': energy, 'free_energy': free_energy} 

65 

66 def read_forces_stress(self): 

67 """Read the forces and stress from the FORCE_STRESS file. 

68 """ 

69 fname = self.directory / 'FORCE_STRESS' 

70 with open(fname) as fd: 

71 lines = fd.readlines() 

72 

73 stress_lines = lines[1:4] 

74 stress = np.empty((3, 3)) 

75 for i in range(3): 

76 line = stress_lines[i].strip().split(' ') 

77 line = [s for s in line if len(s) > 0] 

78 stress[i] = [float(s) for s in line] 

79 

80 results = {} 

81 results['stress'] = np.array( 

82 [stress[0, 0], stress[1, 1], stress[2, 2], 

83 stress[1, 2], stress[0, 2], stress[0, 1]]) 

84 

85 results['stress'] *= Ry / Bohr**3 

86 

87 start = 5 

88 results['forces'] = np.zeros((len(lines) - start, 3), float) 

89 for i in range(start, len(lines)): 

90 line = [s for s in lines[i].strip().split(' ') if len(s) > 0] 

91 results['forces'][i - start] = [float(s) for s in line[2:5]] 

92 

93 results['forces'] *= Ry / Bohr 

94 return results 

95 

96 def read_eigenvalues(self): 

97 """ A robust procedure using the suggestion by Federico Marchesin """ 

98 

99 file_name = self._prefixed('EIG') 

100 try: 

101 with open(file_name) as fd: 

102 fermi_energy = float(fd.readline()) 

103 n, num_hamilton_dim, nkp = map(int, fd.readline().split()) 

104 _ee = np.split( 

105 np.array(fd.read().split()).astype(float), nkp) 

106 except OSError: 

107 return {} 

108 

109 n_spin = 1 if num_hamilton_dim > 2 else num_hamilton_dim 

110 ksn2e = np.delete(_ee, 0, 1).reshape([nkp, n_spin, n]) 

111 

112 eig_array = np.empty((n_spin, nkp, n)) 

113 eig_array[:] = np.inf 

114 

115 for k, sn2e in enumerate(ksn2e): 

116 for s, n2e in enumerate(sn2e): 

117 eig_array[s, k, :] = n2e 

118 

119 assert np.isfinite(eig_array).all() 

120 return {'eigenvalues': eig_array, 'fermi_energy': fermi_energy} 

121 

122 def read_kpoints(self): 

123 """ Reader of the .KP files """ 

124 

125 fname = self._prefixed('KP') 

126 try: 

127 with open(fname) as fd: 

128 nkp = int(next(fd)) 

129 kpoints = np.empty((nkp, 3)) 

130 kweights = np.empty(nkp) 

131 

132 for i in range(nkp): 

133 line = next(fd) 

134 tokens = line.split() 

135 numbers = np.array(tokens[1:]).astype(float) 

136 kpoints[i] = numbers[:3] 

137 kweights[i] = numbers[3] 

138 except OSError: 

139 return {} 

140 

141 return {'kpoints': kpoints, 'kpoint_weights': kweights} 

142 

143 def read_dipole(self): 

144 """Read dipole moment. """ 

145 dipole = np.zeros([1, 3]) 

146 with open(self._prefixed('out')) as fd: 

147 for line in fd: 

148 if line.rfind('Electric dipole (Debye)') > -1: 

149 dipole = np.array([float(f) for f in line.split()[5:8]]) 

150 # debye to e*Ang 

151 return dipole * 0.2081943482534 

152 

153 

154def read_bands_file(fd): 

155 efermi = float(next(fd)) 

156 next(fd) # Appears to be max/min energy. Not important for us 

157 header = next(fd) # Array shape: nbands, nspins, nkpoints 

158 nbands, nspins, nkpts = np.array(header.split()).astype(int) 

159 

160 # three fields for kpt coords, then all the energies 

161 ntokens = nbands * nspins + 3 

162 

163 # Read energies for each kpoint: 

164 data = [] 

165 for _ in range(nkpts): 

166 line = next(fd) 

167 tokens = line.split() 

168 while len(tokens) < ntokens: 

169 # Multirow table. Keep adding lines until the table ends, 

170 # which should happen exactly when we have all the energies 

171 # for this kpoint. 

172 line = next(fd) 

173 tokens += line.split() 

174 assert len(tokens) == ntokens 

175 values = np.array(tokens).astype(float) 

176 data.append(values) 

177 

178 data = np.array(data) 

179 assert len(data) == nkpts 

180 kpts = data[:, :3] 

181 energies = data[:, 3:] 

182 energies = energies.reshape(nkpts, nspins, nbands) 

183 assert energies.shape == (nkpts, nspins, nbands) 

184 return kpts, energies, efermi 

185 

186 

187def resolve_band_structure(path, kpts, energies, efermi): 

188 """Convert input BandPath along with Siesta outputs into BS object.""" 

189 # Right now this function doesn't do much. 

190 # 

191 # Not sure how the output kpoints in the siesta.bands file are derived. 

192 # They appear to be related to the lattice parameter. 

193 # 

194 # We should verify that they are consistent with our input path, 

195 # but since their meaning is unclear, we can't quite do so. 

196 # 

197 # Also we should perhaps verify the cell. If we had the cell, we 

198 # could construct the bandpath from scratch (i.e., pure outputs). 

199 from ase.spectrum.band_structure import BandStructure 

200 ksn2e = energies 

201 skn2e = np.swapaxes(ksn2e, 0, 1) 

202 bs = BandStructure(path, skn2e, reference=efermi) 

203 return bs