Coverage for ase / io / gpaw_out.py: 83.80%
216 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-21 15:52 +0000
1# fmt: off
3import re
5import numpy as np
7from ase.atoms import Atoms
8from ase.calculators.singlepoint import (
9 SinglePointDFTCalculator,
10 SinglePointKPoint,
11)
14def index_startswith(lines: list[str], string: str) -> int:
15 for i, line in enumerate(lines):
16 if line.strip().startswith(string):
17 return i
18 raise ValueError
21def index_pattern(lines: list[str], pattern: str) -> int:
22 repat = re.compile(pattern)
23 for i, line in enumerate(lines):
24 if repat.match(line):
25 return i
26 raise ValueError
29def read_forces(lines: list[str],
30 ii: int,
31 atoms: Atoms) -> tuple[list[tuple[float, float, float]], int]:
32 f = []
33 for i in range(ii + 1, ii + 1 + len(atoms)):
34 try:
35 x, y, z = lines[i].split()[-3:]
36 f.append((float(x), float(y), float(z)))
37 except (ValueError, IndexError) as m:
38 raise OSError(f'Malformed GPAW log file: {m}')
39 return f, i
42def read_stresses(lines: list[str],
43 ii: int,) -> tuple[list[tuple[float, float, float]], int]:
44 s = []
45 for i in range(ii + 1, ii + 4):
46 try:
47 x, y, z = lines[i].split()[-3:]
48 s.append((float(x), float(y), float(z)))
49 except (ValueError, IndexError) as m:
50 raise OSError(f'Malformed GPAW log file: {m}')
51 return s, i
54def find_blocks(lines: list[str]) -> list[list[str]]:
55 blocks = []
56 i1 = 0
57 for i2, line in enumerate(lines):
58 if line == 'positions:\n':
59 if i1 > 0:
60 blocks.append(lines[i1:i2])
61 i1 = i2
62 blocks.append(lines[i1:])
63 return blocks
66def read_gpaw_out(fileobj, index): # -> Union[Atoms, List[Atoms]]:
67 """Read text output from GPAW calculation."""
68 lines = [line.lower() for line in fileobj.readlines()]
70 try:
71 ii = index_startswith(lines, 'total charge:')
72 except ValueError:
73 q = None
74 else:
75 q = float(lines[ii].split()[2])
77 images = []
78 for i, block in enumerate(find_blocks(lines)):
79 # XXX wrong to pass the charge here, it could differ in each block!
80 atoms = process_block(block, q, is_first=not images)
81 if atoms is None:
82 break
84 images.append(atoms)
86 if len(images) == 0:
87 raise OSError('Corrupted GPAW-text file!')
89 return images[index]
92def process_block(lines, q, is_first):
93 try:
94 i = lines.index('unit cell:\n')
95 except ValueError:
96 pass
97 else:
98 if lines[i + 2].startswith(' -'):
99 del lines[i + 2] # old format
100 cell: list[float | list[float]] = []
101 pbc = []
102 for line in lines[i + 2:i + 5]:
103 words = line.split()
104 if len(words) == 5: # old format
105 cell.append(float(words[2]))
106 pbc.append(words[1] == 'yes')
107 else: # new format with GUC
108 cell.append([float(word) for word in words[3:6]])
109 pbc.append(words[2] == 'yes')
111 symbols = []
112 positions = []
113 magmoms = []
114 for line in lines[1:]:
115 words = line.split()
116 if len(words) < 5:
117 break
118 n, symbol, x, y, z = words[:5]
119 symbols.append(symbol.split('.')[0].title())
120 positions.append([float(x), float(y), float(z)])
121 if len(words) > 5:
122 mom = float(words[-1].rstrip(')'))
123 magmoms.append(mom)
124 if len(symbols):
125 atoms = Atoms(symbols=symbols, positions=positions,
126 cell=cell, pbc=pbc)
127 else:
128 atoms = Atoms(cell=cell, pbc=pbc)
130 try:
131 ii = index_pattern(lines, '\\d+ k-point')
132 word = lines[ii].split()
133 kx = int(word[2])
134 ky = int(word[4])
135 kz = int(word[6])
136 bz_kpts = (kx, ky, kz)
137 ibz_kpts = int(lines[ii + 1].split()[0])
138 except (ValueError, TypeError, IndexError):
139 bz_kpts = None
140 ibz_kpts = None
142 try:
143 i = index_startswith(lines, 'energy contributions relative to') + 2
144 except ValueError:
145 e = energy_contributions = None
146 else:
147 energy_contributions = {}
148 while i < len(lines):
149 fields = lines[i].split(':')
150 if len(fields) == 2:
151 name = fields[0]
152 energy = float(fields[1])
153 energy_contributions[name] = energy
154 if name in ['zero kelvin', 'extrapolated']:
155 e = energy
156 break
157 i += 1
158 else: # no break
159 raise ValueError
161 try:
162 ii = index_pattern(lines, '(fixed )?fermi level(s)?:')
163 except ValueError:
164 eFermi = None
165 else:
166 fields = lines[ii].split()
167 try:
168 def strip(string):
169 for rubbish in '[],':
170 string = string.replace(rubbish, '')
171 return string
172 eFermi = [float(strip(fields[-2])),
173 float(strip(fields[-1]))]
174 except ValueError:
175 eFermi = float(fields[-1])
177 # read Eigenvalues and occupations
178 ii1 = ii2 = 1e32
179 try:
180 ii1 = index_startswith(lines, 'band eigenvalues occupancy')
181 except ValueError:
182 pass
183 try:
184 ii2 = index_startswith(lines, 'band eigenvalues occupancy')
185 except ValueError:
186 pass
187 ii = min(ii1, ii2)
188 if ii == 1e32:
189 kpts = None
190 else:
191 ii += 1
192 words = lines[ii].split()
193 vals = []
194 while len(words) > 2:
195 vals.append([float(w) for w in words])
196 ii += 1
197 words = lines[ii].split()
198 vals = np.array(vals).transpose()
199 kpts = [SinglePointKPoint(1, 0, 0)]
200 kpts[0].eps_n = vals[1]
201 kpts[0].f_n = vals[2]
202 if vals.shape[0] > 3:
203 kpts.append(SinglePointKPoint(1, 1, 0))
204 kpts[1].eps_n = vals[3]
205 kpts[1].f_n = vals[4]
206 # read dipole moment
207 try:
208 ii = index_startswith(lines, 'dipole moment:')
209 except ValueError:
210 dipole = None
211 else:
212 line = lines[ii]
213 for x in '()[],':
214 line = line.replace(x, '')
215 dipole = np.array([float(c) for c in line.split()[2:5]])
217 try:
218 ii = index_startswith(lines, 'local magnetic moments')
219 except ValueError:
220 pass
221 else:
222 magmoms = []
223 for j in range(ii + 1, ii + 1 + len(atoms)):
224 line = lines[j]
225 if '#' in line: # new GPAW format
226 magmom = line.split()[-4].split(']')[0]
227 else:
228 magmom = line.split()[-1].rstrip(')')
229 magmoms.append(float(magmom))
231 try:
232 ii = lines.index('forces in ev/ang:\n')
233 except ValueError:
234 f = None
235 else:
236 f, i = read_forces(lines, ii, atoms)
238 try:
239 ii = lines.index('stress tensor:\n')
240 except ValueError:
241 stress_tensor = None
242 else:
243 stress_tensor, i = read_stresses(lines, ii)
245 try:
246 parameters = {}
247 ii = index_startswith(lines, 'vdw correction:')
248 except ValueError:
249 name = 'gpaw'
250 else:
251 name = lines[ii - 1].strip()
252 # save uncorrected values
253 parameters.update({
254 'calculator': 'gpaw',
255 'uncorrected_energy': e,
256 })
257 line = lines[ii + 1]
258 assert line.startswith('energy:')
259 e = float(line.split()[-1])
260 f, i = read_forces(lines, ii + 3, atoms)
262 if not is_first and e is None:
263 return None
265 if q is not None and len(atoms) > 0:
266 n = len(atoms)
267 atoms.set_initial_charges([q / n] * n)
269 if magmoms:
270 atoms.set_initial_magnetic_moments(magmoms)
271 else:
272 magmoms = None
273 if e is not None or f is not None:
274 calc = SinglePointDFTCalculator(atoms, energy=e, forces=f,
275 dipole=dipole, magmoms=magmoms,
276 efermi=eFermi,
277 bzkpts=bz_kpts, ibzkpts=ibz_kpts,
278 stress=stress_tensor)
279 calc.name = name
280 calc.parameters = parameters
281 if energy_contributions is not None:
282 calc.energy_contributions = energy_contributions
283 if kpts is not None:
284 calc.kpts = kpts
285 atoms.calc = calc
287 return atoms