Coverage for /builds/kinetik161/ase/ase/io/octopus/output.py: 48.50%
167 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1import re
3import numpy as np
5from ase.units import Angstrom, Bohr, Debye, Hartree, eV
8class OctopusIOError(IOError):
9 pass # Cannot find output files
12def read_eigenvalues_file(fd):
13 unit = None
15 for line in fd:
16 m = re.match(r'Eigenvalues\s*\[(.+?)\]', line)
17 if m is not None:
18 unit = m.group(1)
19 break
20 line = next(fd)
21 assert line.strip().startswith('#st'), line
23 # fermilevel = None
24 kpts = []
25 eigs = []
26 occs = []
28 for line in fd:
29 m = re.match(r'#k.*?\(\s*(.+?),\s*(.+?),\s*(.+?)\)', line)
30 if m:
31 k = m.group(1, 2, 3)
32 kpts.append(np.array(k, float))
33 eigs.append({})
34 occs.append({})
35 else:
36 m = re.match(r'\s*\d+\s*(\S+)\s*(\S+)\s*(\S+)', line)
37 if m is None:
38 m = re.match(r'Fermi energy\s*=\s*(\S+)\s*', line)
39 assert m is not None
40 # We can also return the fermilevel but so far we just read
41 # it from the static/info instead.
42 # fermilevel = float(m.group(1))
43 else:
44 spin, eig, occ = m.group(1, 2, 3)
46 if not eigs:
47 # Only initialized if kpoint header was written
48 eigs.append({})
49 occs.append({})
51 eigs[-1].setdefault(spin, []).append(float(eig))
52 occs[-1].setdefault(spin, []).append(float(occ))
54 nkpts = len(kpts)
55 nspins = len(eigs[0])
56 nbands = len(eigs[0][spin])
58 kptsarr = np.array(kpts, float)
59 eigsarr = np.empty((nkpts, nspins, nbands))
60 occsarr = np.empty((nkpts, nspins, nbands))
62 arrs = [eigsarr, occsarr]
64 for arr in arrs:
65 arr.fill(np.nan)
67 for k in range(nkpts):
68 for arr, lst in [(eigsarr, eigs), (occsarr, occs)]:
69 arr[k, :, :] = [lst[k][sp] for sp
70 in (['--'] if nspins == 1 else ['up', 'dn'])]
72 for arr in arrs:
73 assert not np.isnan(arr).any()
75 eigsarr *= {'H': Hartree, 'eV': eV}[unit]
76 return kptsarr, eigsarr, occsarr
79def read_static_info_stress(fd):
80 stress_cv = np.empty((3, 3))
82 headers = next(fd)
83 assert headers.strip().startswith('T_{ij}')
84 for i in range(3):
85 line = next(fd)
86 tokens = line.split()
87 vec = np.array(tokens[1:4]).astype(float)
88 stress_cv[i] = vec
89 return stress_cv
92def read_static_info_kpoints(fd):
93 for line in fd:
94 if line.startswith('List of k-points'):
95 break
97 tokens = next(fd).split()
98 assert tokens == ['ik', 'k_x', 'k_y', 'k_z', 'Weight']
99 bar = next(fd)
100 assert bar.startswith('---')
102 kpts = []
103 weights = []
105 for line in fd:
106 # Format: index kx ky kz weight
107 m = re.match(r'\s*\d+\s*(\S+)\s*(\S+)\s*(\S+)\s*(\S+)', line)
108 if m is None:
109 break
110 kxyz = m.group(1, 2, 3)
111 weight = m.group(4)
112 kpts.append(kxyz)
113 weights.append(weight)
115 ibz_k_points = np.array(kpts, float)
116 k_point_weights = np.array(weights, float)
117 return dict(ibz_k_points=ibz_k_points, k_point_weights=k_point_weights)
120def read_static_info_eigenvalues(fd, energy_unit):
122 values_sknx = {}
124 nbands = 0
125 fermilevel = None
126 for line in fd:
127 line = line.strip()
128 if line.startswith('#'):
129 continue
130 if not line[:1].isdigit():
131 m = re.match(r'Fermi energy\s*=\s*(\S+)', line)
132 if m is not None:
133 fermilevel = float(m.group(1)) * energy_unit
134 break
136 tokens = line.split()
137 nbands = max(nbands, int(tokens[0]))
138 energy = float(tokens[2]) * energy_unit
139 occupation = float(tokens[3])
140 values_sknx.setdefault(tokens[1], []).append((energy, occupation))
142 nspins = len(values_sknx)
143 if nspins == 1:
144 val = [values_sknx['--']]
145 else:
146 val = [values_sknx['up'], values_sknx['dn']]
147 val = np.array(val, float)
148 nkpts, remainder = divmod(len(val[0]), nbands)
149 assert remainder == 0
151 eps_skn = val[:, :, 0].reshape(nspins, nkpts, nbands)
152 occ_skn = val[:, :, 1].reshape(nspins, nkpts, nbands)
153 eps_skn = eps_skn.transpose(1, 0, 2).copy()
154 occ_skn = occ_skn.transpose(1, 0, 2).copy()
155 assert eps_skn.flags.contiguous
156 d = dict(nspins=nspins,
157 nkpts=nkpts,
158 nbands=nbands,
159 eigenvalues=eps_skn,
160 occupations=occ_skn)
161 if fermilevel is not None:
162 d.update(fermi_level=fermilevel)
163 return d
166def read_static_info_energy(fd, energy_unit):
167 def get(name):
168 for line in fd:
169 if line.strip().startswith(name):
170 return float(line.split('=')[-1].strip()) * energy_unit
171 return dict(energy=get('Total'), free_energy=get('Free'))
174def read_static_info(fd):
175 results = {}
177 def get_energy_unit(line): # Convert "title [unit]": ---> unit
178 return {'[eV]': eV, '[H]': Hartree}[line.split()[1].rstrip(':')]
180 for line in fd:
181 if line.strip('*').strip().startswith('Brillouin zone'):
182 results.update(read_static_info_kpoints(fd))
183 elif line.startswith('Eigenvalues ['):
184 unit = get_energy_unit(line)
185 results.update(read_static_info_eigenvalues(fd, unit))
186 elif line.startswith('Energy ['):
187 unit = get_energy_unit(line)
188 results.update(read_static_info_energy(fd, unit))
189 elif line.startswith('Stress tensor'):
190 assert line.split()[-1] == '[H/b^3]'
191 stress = read_static_info_stress(fd)
192 stress *= Hartree / Bohr**3
193 results.update(stress=stress)
194 elif line.startswith('Total Magnetic Moment'):
195 if 0:
196 line = next(fd)
197 values = line.split()
198 results['magmom'] = float(values[-1])
200 line = next(fd)
201 assert line.startswith('Local Magnetic Moments')
202 line = next(fd)
203 assert line.split() == ['Ion', 'mz']
204 # Reading Local Magnetic Moments
205 mag_moment = []
206 for line in fd:
207 if line == '\n':
208 break # there is no more thing to search for
209 line = line.replace('\n', ' ')
210 values = line.split()
211 mag_moment.append(float(values[-1]))
213 results['magmoms'] = np.array(mag_moment)
214 elif line.startswith('Dipole'):
215 assert line.split()[-1] == '[Debye]'
216 dipole = [float(next(fd).split()[-1]) for i in range(3)]
217 results['dipole'] = np.array(dipole) * Debye
218 elif line.startswith('Forces'):
219 forceunitspec = line.split()[-1]
220 forceunit = {'[eV/A]': eV / Angstrom,
221 '[H/b]': Hartree / Bohr}[forceunitspec]
222 forces = []
223 line = next(fd)
224 assert line.strip().startswith('Ion')
225 for line in fd:
226 if line.strip().startswith('---'):
227 break
228 tokens = line.split()[-3:]
229 forces.append([float(f) for f in tokens])
230 results['forces'] = np.array(forces) * forceunit
231 elif line.startswith('Fermi'):
232 tokens = line.split()
233 unit = {'eV': eV, 'H': Hartree}[tokens[-1]]
234 eFermi = float(tokens[-2]) * unit
235 results['efermi'] = eFermi
237 if 'ibz_k_points' not in results:
238 results['ibz_k_points'] = np.zeros((1, 3))
239 results['k_point_weights'] = np.ones(1)
240 if 0: # 'efermi' not in results:
241 # Find HOMO level. Note: This could be a very bad
242 # implementation with fractional occupations if the Fermi
243 # level was not found otherwise.
244 all_energies = results['eigenvalues'].ravel()
245 all_occupations = results['occupations'].ravel()
246 args = np.argsort(all_energies)
247 for arg in args[::-1]:
248 if all_occupations[arg] > 0.1:
249 break
250 eFermi = all_energies[arg]
251 results['efermi'] = eFermi
253 return results