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