Coverage for /builds/ase/ase/ase/io/siesta_output.py: 97.76%

1# fmt: off 

2 

3import numpy as np 

4 

5from ase.units import Bohr, Ry 

6 

7 

8class OutputReader: 

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

10 self.prefix = prefix 

11 self.directory = directory 

12 self.bandpath = bandpath 

13 

14 def read_results(self): 

15 results = {} 

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

17 results.update(self.read_energy()) 

18 results.update(self.read_forces_stress()) 

19 results.update(self.read_eigenvalues()) 

20 results.update(self.read_kpoints()) 

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

22 

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

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

25 

26 return results 

27 

28 def _prefixed(self, extension): 

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

30 

31 def read_bands(self, bandpath): 

32 fname = self._prefixed('bands') 

33 with open(fname) as fd: 

34 kpts, energies, efermi = read_bands_file(fd) 

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

36 

37 def read_number_of_grid_points(self): 

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

39 

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

41 with open(fname) as fd: 

42 for line in fd: 

43 line = line.strip().lower() 

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

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

46 

47 raise RuntimeError 

48 

49 def read_energy(self): 

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

51 """ 

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

53 

54 assert 'final energy' in text 

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

56 

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

58 for line in lines: 

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

60 if has_energy: 

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

62 line = next(lines) 

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

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

65 

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

67 

68 def read_forces_stress(self): 

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

70 """ 

71 fname = self.directory / 'FORCE_STRESS' 

72 with open(fname) as fd: 

73 lines = fd.readlines() 

74 

75 stress_lines = lines[1:4] 

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

77 for i in range(3): 

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

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

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

81 

82 results = {} 

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

84 [stress[0, 0], stress[1, 1], stress[2, 2], 

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

86 

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

88 

89 start = 5 

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

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

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

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

94 

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

96 return results 

97 

98 def read_eigenvalues(self): 

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

100 

101 file_name = self._prefixed('EIG') 

102 try: 

103 with open(file_name) as fd: 

104 fermi_energy = float(fd.readline()) 

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

106 _ee = np.split( 

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

108 except OSError: 

109 return {} 

110 

111 n_spin = 1 if num_hamilton_dim > 2 else num_hamilton_dim 

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

113 

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

115 eig_array[:] = np.inf 

116 

117 for k, sn2e in enumerate(ksn2e): 

118 for s, n2e in enumerate(sn2e): 

119 eig_array[s, k, :] = n2e 

120 

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

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

123 

124 def read_kpoints(self): 

125 """ Reader of the .KP files """ 

126 

127 fname = self._prefixed('KP') 

128 try: 

129 with open(fname) as fd: 

130 nkp = int(next(fd)) 

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

132 kweights = np.empty(nkp) 

133 

134 for i in range(nkp): 

135 line = next(fd) 

136 tokens = line.split() 

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

138 kpoints[i] = numbers[:3] 

139 kweights[i] = numbers[3] 

140 except OSError: 

141 return {} 

142 

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

144 

145 def read_dipole(self): 

146 """Read dipole moment. """ 

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

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

149 for line in fd: 

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

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

152 # debye to e*Ang 

153 return dipole * 0.2081943482534 

154 

155 

156def read_bands_file(fd): 

157 efermi = float(next(fd)) 

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

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

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

161 

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

163 ntokens = nbands * nspins + 3 

164 

165 # Read energies for each kpoint: 

166 data = [] 

167 for _ in range(nkpts): 

168 line = next(fd) 

169 tokens = line.split() 

170 while len(tokens) < ntokens: 

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

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

173 # for this kpoint. 

174 line = next(fd) 

175 tokens += line.split() 

176 assert len(tokens) == ntokens 

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

178 data.append(values) 

179 

180 data = np.array(data) 

181 assert len(data) == nkpts 

182 kpts = data[:, :3] 

183 energies = data[:, 3:] 

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

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

186 return kpts, energies, efermi 

187 

188 

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

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

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

192 # 

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

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

195 # 

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

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

198 # 

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

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

201 from ase.spectrum.band_structure import BandStructure 

202 ksn2e = energies 

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

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

205 return bs