Coverage for /builds/hweiske/ase/ase/io/nwchem/nwreader.py: 65.57%

1import re

2from collections import OrderedDict

4import numpy as np

6from ase import Atoms

7from ase.calculators.singlepoint import (SinglePointDFTCalculator,

8 SinglePointKPoint)

9from ase.units import Bohr, Hartree

11from .parser import _define_pattern

13# Note to the reader of this code: Here and below we use the function

14# _define_pattern from parser.py in this same directory to compile

15# regular expressions. These compiled expressions are stored along with

16# an example string that the expression should match in a list that

17# is used during tests (test/nwchem/nwchem_parser.py) to ensure that

18# the regular expressions are still working correctly.

20# Matches the beginning of a GTO calculation

21_gauss_block = _define_pattern(

22 r'^[\s]+NWChem (?:SCF|DFT) Module\n$',

23 " NWChem SCF Module\n",

24)

27# Matches the beginning of a plane wave calculation

28_pw_block = _define_pattern(

29 r'^[\s]+\*[\s]+NWPW (?:PSPW|BAND|PAW|Band Structure) Calculation'

30 r'[\s]+\*[\s]*\n$',

31 " * NWPW PSPW Calculation *\n",

32)

35# Top-level parser

36def read_nwchem_out(fobj, index=-1):

37 """Splits an NWChem output file into chunks corresponding to

38 individual single point calculations."""

39 lines = fobj.readlines()

41 if index == slice(-1, None, None):

42 for line in lines:

43 if _gauss_block.match(line):

44 return [parse_gto_chunk(''.join(lines))]

45 if _pw_block.match(line):

46 return [parse_pw_chunk(''.join(lines))]

47 raise ValueError('This does not appear to be a valid NWChem '

48 'output file.')

50 # First, find each SCF block

51 group = []

52 atomslist = []

53 header = True

54 lastgroup = []

55 lastparser = None

56 parser = None

57 for line in lines:

58 group.append(line)

59 if _gauss_block.match(line):

60 next_parser = parse_gto_chunk

61 elif _pw_block.match(line):

62 next_parser = parse_pw_chunk

63 else:

64 continue

66 if header:

67 header = False

68 else:

69 atoms = parser(''.join(group))

70 if atoms is None and parser is lastparser:

71 atoms = parser(''.join(lastgroup + group))

72 if atoms is not None:

73 atomslist[-1] = atoms

74 lastgroup += group

75 else:

76 atomslist.append(atoms)

77 lastgroup = group

78 lastparser = parser

79 group = []

80 parser = next_parser

81 if not header:

82 atoms = parser(''.join(group))

83 if atoms is not None:

84 atomslist.append(atoms)

86 return atomslist[index]

89# Matches a geometry block and returns the geometry specification lines

90_geom = _define_pattern(

91 r'\n[ \t]+Geometry \"[ \t\S]+\" -> \"[ \t\S]*\"[ \t]*\n'

92 r'^[ \t-]+\n'

93 r'(?:^[ \t\S]*\n){3}'

94 r'^[ \t]+No\.[ \t]+Tag[ \t]+Charge[ \t]+X[ \t]+Y[ \t]+Z\n'

95 r'^[ \t-]+\n'

96 r'((?:^(?:[ \t]+[\S]+){6}[ \t]*\n)+)',

97 """\

99 Geometry "geometry" -> ""

100 -------------------------

101

102 Output coordinates in angstroms (scale by 1.889725989 to convert to a.u.)

103

104 No. Tag Charge X Y Z

105 ---- ---------------- ---------- -------------- -------------- --------------

106 1 C 6.0000 0.00000000 0.00000000 0.00000000

107 2 H 1.0000 0.62911800 0.62911800 0.62911800

108 3 H 1.0000 -0.62911800 -0.62911800 0.62911800

109 4 H 1.0000 0.62911800 -0.62911800 -0.62911800

110""", re.M)

111

112# Unit cell parser

113_cell_block = _define_pattern(r'^[ \t]+Lattice Parameters[ \t]*\n'

114 r'^(?:[ \t\S]*\n){4}'

115 r'((?:^(?:[ \t]+[\S]+){5}\n){3})',

116 """\

117 Lattice Parameters

118 ------------------

119

120 lattice vectors in angstroms (scale by 1.889725989 to convert to a.u.)

121

122 a1=< 4.000 0.000 0.000 >

123 a2=< 0.000 5.526 0.000 >

124 a3=< 0.000 0.000 4.596 >

125 a= 4.000 b= 5.526 c= 4.596

126 alpha= 90.000 beta= 90.000 gamma= 90.000

127 omega= 101.6

128""", re.M)

129

130

131# Parses the geometry and returns the corresponding Atoms object

132def _parse_geomblock(chunk):

133 geomblocks = _geom.findall(chunk)

134 if not geomblocks:

135 return None

136 geomblock = geomblocks[-1].strip().split('\n')

137 natoms = len(geomblock)

138 symbols = []

139 pos = np.zeros((natoms, 3))

140 for i, line in enumerate(geomblock):

141 line = line.strip().split()

142 symbols.append(line[1])

143 pos[i] = [float(x) for x in line[3:6]]

144

145 cellblocks = _cell_block.findall(chunk)

146 if cellblocks:

147 cellblock = cellblocks[-1].strip().split('\n')

148 cell = np.zeros((3, 3))

149 for i, line in enumerate(cellblock):

150 line = line.strip().split()

151 cell[i] = [float(x) for x in line[1:4]]

152 else:

153 cell = None

154 return Atoms(symbols, positions=pos, cell=cell)

155

156

157# GTO-specific parser stuff

158

159# Matches gradient block from a GTO calculation

160_gto_grad = _define_pattern(

161 r'^[ \t]+[\S]+[ \t]+ENERGY GRADIENTS[ \t]*[\n]+'

162 r'^[ \t]+atom[ \t]+coordinates[ \t]+gradient[ \t]*\n'

163 r'^(?:[ \t]+x[ \t]+y[ \t]+z){2}[ \t]*\n'

164 r'((?:^(?:[ \t]+[\S]+){8}\n)+)[ \t]*\n',

165 """\

166 UHF ENERGY GRADIENTS

167

168 atom coordinates gradient

169 x y z x y z

170 1 C 0.293457 -0.293457 0.293457 -0.000083 0.000083 -0.000083

171 2 H 1.125380 1.355351 1.125380 0.000086 0.000089 0.000086

172 3 H -1.355351 -1.125380 1.125380 -0.000089 -0.000086 0.000086

173 4 H 1.125380 -1.125380 -1.355351 0.000086 -0.000086 -0.000089

174

175""", re.M)

176

177# Energy parsers for a variety of different GTO calculations

178_e_gto = OrderedDict()

179_e_gto['tce'] = _define_pattern(

180 r'^[\s]+[\S]+[\s]+total energy \/ hartree[\s]+'

181 r'=[\s]+([\S]+)[\s]*\n',

182 " CCD total energy / hartree "

183 "= -75.715332545665888\n", re.M,

184)

185_e_gto['ccsd'] = _define_pattern(

186 r'^[\s]+Total CCSD energy:[\s]+([\S]+)[\s]*\n',

187 " Total CCSD energy: -75.716168566598569\n",

188 re.M,

189)

190_e_gto['tddft'] = _define_pattern(

191 r'^[\s]+Excited state energy =[\s]+([\S]+)[\s]*\n',

192 " Excited state energy = -75.130134499965\n",

193 re.M,

194)

195_e_gto['mp2'] = _define_pattern(

196 r'^[\s]+Total MP2 energy[\s]+([\S]+)[\s]*\n',

197 " Total MP2 energy -75.708800087578\n",

198 re.M,

199)

200_e_gto['mf'] = _define_pattern(

201 r'^[\s]+Total (?:DFT|SCF) energy =[\s]+([\S]+)[\s]*\n',

202 " Total SCF energy = -75.585555997789\n",

203 re.M,

204)

205

206

207# GTO parser

208def parse_gto_chunk(chunk):

209 atoms = None

210 forces = None

211 energy = None

212 dipole = None

213 quadrupole = None

214 for theory, pattern in _e_gto.items():

215 matches = pattern.findall(chunk)

216 if matches:

217 energy = float(matches[-1].replace('D', 'E')) * Hartree

218 break

219

220 gradblocks = _gto_grad.findall(chunk)

221 if gradblocks:

222 gradblock = gradblocks[-1].strip().split('\n')

223 natoms = len(gradblock)

224 symbols = []

225 pos = np.zeros((natoms, 3))

226 forces = np.zeros((natoms, 3))

227 for i, line in enumerate(gradblock):

228 line = line.strip().split()

229 symbols.append(line[1])

230 pos[i] = [float(x) for x in line[2:5]]

231 forces[i] = [-float(x) for x in line[5:8]]

232 pos *= Bohr

233 forces *= Hartree / Bohr

234 atoms = Atoms(symbols, positions=pos)

235

236 dipole, quadrupole = _get_multipole(chunk)

237

238 kpts = _get_gto_kpts(chunk)

239

240 if atoms is None:

241 atoms = _parse_geomblock(chunk)

242

243 if atoms is None:

244 return None

245

246 # SinglePointDFTCalculator doesn't support quadrupole moment currently

247 calc = SinglePointDFTCalculator(atoms=atoms,

248 energy=energy,

249 free_energy=energy, # XXX Is this right?

250 forces=forces,

251 dipole=dipole,

252 # quadrupole=quadrupole,

253 )

254 calc.kpts = kpts

255 atoms.calc = calc

256 return atoms

257

258

259# Extracts dipole and quadrupole moment for a GTO calculation

260# Note on the regex: Some, but not all, versions of NWChem

261# insert extra spaces in the blank lines. Do not remove the \s*

262# in between \n and \n

263_multipole = _define_pattern(

264 r'^[ \t]+Multipole analysis of the density[ \t\S]*\n'

265 r'^[ \t-]+\n\s*\n^[ \t\S]+\n^[ \t-]+\n'

266 r'((?:(?:(?:[ \t]+[\S]+){7,8}\n)|[ \t]*\n){12})',

267 """\

268 Multipole analysis of the density

269 ---------------------------------

270

271 L x y z total alpha beta nuclear

272 - - - - ----- ----- ---- -------

273 0 0 0 0 -0.000000 -5.000000 -5.000000 10.000000

274

275 1 1 0 0 0.000000 0.000000 0.000000 0.000000

276 1 0 1 0 -0.000001 -0.000017 -0.000017 0.000034

277 1 0 0 1 -0.902084 -0.559881 -0.559881 0.217679

278

279 2 2 0 0 -5.142958 -2.571479 -2.571479 0.000000

280 2 1 1 0 -0.000000 -0.000000 -0.000000 0.000000

281 2 1 0 1 0.000000 0.000000 0.000000 0.000000

282 2 0 2 0 -3.153324 -3.807308 -3.807308 4.461291

283 2 0 1 1 0.000001 -0.000009 -0.000009 0.000020

284 2 0 0 2 -4.384288 -3.296205 -3.296205 2.208122

285""", re.M)

286

287

288# Parses the dipole and quadrupole moment from a GTO calculation

289def _get_multipole(chunk):

290 matches = _multipole.findall(chunk)

291 if not matches:

292 return None, None

293 # This pulls the 5th column out of the multipole moments block;

294 # this column contains the actual moments.

295 moments = [float(x.split()[4]) for x in matches[-1].split('\n')

296 if x and not x.isspace()]

297 dipole = np.array(moments[1:4]) * Bohr

298 quadrupole = np.zeros(9)

299 quadrupole[[0, 1, 2, 4, 5, 8]] = [moments[4:]]

300 quadrupole[[3, 6, 7]] = quadrupole[[1, 2, 5]]

301 return dipole, quadrupole.reshape((3, 3)) * Bohr**2

302

303

304# MO eigenvalue and occupancy parser for GTO calculations

305_eval_block = _define_pattern(

306 r'^[ \t]+[\S]+ Final (?:Alpha |Beta )?Molecular Orbital Analysis[ \t]*'

307 r'\n^[ \t-]+\n\n'

308 r'(?:^[ \t]+Vector [ \t\S]+\n(?:^[ \t\S]+\n){3}'

309 r'(?:^(?:(?:[ \t]+[\S]+){5}){1,2}[ \t]*\n)+\n)+',

310 """\

311 ROHF Final Molecular Orbital Analysis

312 -------------------------------------

313

314 Vector 1 Occ=2.000000D+00 E=-2.043101D+01

315 MO Center= 1.1D-20, 1.5D-18, 1.2D-01, r^2= 1.5D-02

316 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function

317 ----- ------------ --------------- ----- ------------ ---------------

318 1 0.983233 1 O s

319

320 Vector 2 Occ=2.000000D+00 E=-1.324439D+00

321 MO Center= -2.1D-18, -8.6D-17, -7.1D-02, r^2= 5.1D-01

322 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function

323 ----- ------------ --------------- ----- ------------ ---------------

324 6 0.708998 1 O s 1 -0.229426 1 O s

325 2 0.217752 1 O s

326 """, re.M) # noqa: W291

327

328

329# Parses the eigenvalues and occupations from a GTO calculation

330def _get_gto_kpts(chunk):

331 eval_blocks = _eval_block.findall(chunk)

332 if not eval_blocks:

333 return []

334 kpts = []

335 kpt = _get_gto_evals(eval_blocks[-1])

336 if kpt.s == 1:

337 kpts.append(_get_gto_evals(eval_blocks[-2]))

338 kpts.append(kpt)

339 return kpts

340

341

342# Extracts MO eigenvalue and occupancy for a GTO calculation

343_extract_vector = _define_pattern(

344 r'^[ \t]+Vector[ \t]+([\S])+[ \t]+Occ=([\S]+)[ \t]+E=[ \t]*([\S]+)[ \t]*\n',

345 " Vector 1 Occ=2.000000D+00 E=-2.043101D+01\n", re.M,

346)

347

348

349# Extracts the eigenvalues and occupations from a GTO calculation

350def _get_gto_evals(chunk):

351 spin = 1 if re.match(r'[ \t\S]+Beta', chunk) else 0

352 data = []

353 for vector in _extract_vector.finditer(chunk):

354 data.append([float(x.replace('D', 'E')) for x in vector.groups()[1:]])

355 data = np.array(data)

356 occ = data[:, 0]

357 energies = data[:, 1] * Hartree

358

359 return SinglePointKPoint(1., spin, 0, energies, occ)

360

361

362# Plane wave specific parsing stuff

363

364# Matches the gradient block from a plane wave calculation

365_nwpw_grad = _define_pattern(

366 r'^[ \t]+[=]+[ \t]+Ion Gradients[ \t]+[=]+[ \t]*\n'

367 r'^[ \t]+Ion Forces:[ \t]*\n'

368 r'((?:^(?:[ \t]+[\S]+){7}\n)+)',

369 """\

370 ============= Ion Gradients =================

371 Ion Forces:

372 1 O ( -0.000012 0.000027 -0.005199 )

373 2 H ( 0.000047 -0.013082 0.020790 )

374 3 H ( 0.000047 0.012863 0.020786 )

375 C.O.M. ( -0.000000 -0.000000 -0.000000 )

376 ===============================================

377""", re.M)

378

379# Matches the gradient block from a PAW calculation

380_paw_grad = _define_pattern(

381 r'^[ \t]+[=]+[ \t]+Ion Gradients[ \t]+[=]+[ \t]*\n'

382 r'^[ \t]+Ion Positions:[ \t]*\n'

383 r'((?:^(?:[ \t]+[\S]+){7}\n)+)'

384 r'^[ \t]+Ion Forces:[ \t]*\n'

385 r'((?:^(?:[ \t]+[\S]+){7}\n)+)',

386 """\

387 ============= Ion Gradients =================

388 Ion Positions:

389 1 O ( -3.77945 -5.22176 -3.77945 )

390 2 H ( -3.77945 -3.77945 3.77945 )

391 3 H ( -3.77945 3.77945 3.77945 )

392 Ion Forces:

393 1 O ( -0.00001 -0.00000 0.00081 )

394 2 H ( 0.00005 -0.00026 -0.00322 )

395 3 H ( 0.00005 0.00030 -0.00322 )

396 C.O.M. ( -0.00000 -0.00000 -0.00000 )

397 ===============================================

398""", re.M)

399

400# Energy parser for plane wave calculations

401_nwpw_energy = _define_pattern(

402 r'^[\s]+Total (?:PSPW|BAND|PAW) energy'

403 r'[\s]+:[\s]+([\S]+)[\s]*\n',

404 " Total PSPW energy : -0.1709317826E+02\n",

405 re.M,

406)

407

408# Parser for the fermi energy in a plane wave calculation

409_fermi_energy = _define_pattern(

410 r'^[ \t]+Fermi energy =[ \t]+([\S]+) \([ \t]+[\S]+[ \t]*\n',

411 " Fermi energy = -0.5585062E-01 ( -1.520eV)\n", re.M,

412)

413

414

415# Plane wave parser

416def parse_pw_chunk(chunk):

417 atoms = _parse_geomblock(chunk)

418 if atoms is None:

419 return None

420

421 energy = None

422 efermi = None

423 forces = None

424 stress = None

425

426 matches = _nwpw_energy.findall(chunk)

427 if matches:

428 energy = float(matches[-1].replace('D', 'E')) * Hartree

429

430 matches = _fermi_energy.findall(chunk)

431 if matches:

432 efermi = float(matches[-1].replace('D', 'E')) * Hartree

433

434 gradblocks = _nwpw_grad.findall(chunk)

435 if not gradblocks:

436 gradblocks = _paw_grad.findall(chunk)

437 if gradblocks:

438 gradblock = gradblocks[-1].strip().split('\n')

439 natoms = len(gradblock)

440 symbols = []

441 forces = np.zeros((natoms, 3))

442 for i, line in enumerate(gradblock):

443 line = line.strip().split()

444 symbols.append(line[1])

445 forces[i] = [float(x) for x in line[3:6]]

446 forces *= Hartree / Bohr

447

448 if atoms.cell:

449 stress = _get_stress(chunk, atoms.cell)

450

451 ibz_kpts, kpts = _get_pw_kpts(chunk)

452

453 # NWChem does not calculate an energy extrapolated to the 0K limit,

454 # so right now, energy and free_energy will be the same.

455 calc = SinglePointDFTCalculator(atoms=atoms,

456 energy=energy,

457 efermi=efermi,

458 free_energy=energy,

459 forces=forces,

460 stress=stress,

461 ibzkpts=ibz_kpts)

462 calc.kpts = kpts

463 atoms.calc = calc

464 return atoms

465

466

467# Extracts stress tensor from a plane wave calculation

468_stress = _define_pattern(

469 r'[ \t]+[=]+[ \t]+(?:total gradient|E all FD)[ \t]+[=]+[ \t]*\n'

470 r'^[ \t]+S =((?:(?:[ \t]+[\S]+){5}\n){3})[ \t=]+\n',

471 """\

472 ============= total gradient ==============

473 S = ( -0.22668 0.27174 0.19134 )

474 ( 0.23150 -0.26760 0.23226 )

475 ( 0.19090 0.27206 -0.22700 )

476 ===================================================

477""", re.M)

478

479

480# Extract stress tensor from a plane wave calculation

481def _get_stress(chunk, cell):

482 stress_blocks = _stress.findall(chunk)

483 if not stress_blocks:

484 return None

485 stress_block = stress_blocks[-1]

486 stress = np.zeros((3, 3))

487 for i, row in enumerate(stress_block.strip().split('\n')):

488 stress[i] = [float(x) for x in row.split()[1:4]]

489 stress = (stress @ cell) * Hartree / Bohr / cell.volume

490 stress = 0.5 * (stress + stress.T)

491 # convert from 3x3 array to Voigt form

492 return stress.ravel()[[0, 4, 8, 5, 2, 1]]

493

494

495# MO/band eigenvalue and occupancy parser for plane wave calculations

496_nwpw_eval_block = _define_pattern(

497 r'(?:(?:^[ \t]+Brillouin zone point:[ \t]+[\S]+[ \t]*\n'

498 r'(?:[ \t\S]*\n){3,4})?'

499 r'^[ \t]+(?:virtual )?orbital energies:\n'

500 r'(?:^(?:(?:[ \t]+[\S]+){3,4}){1,2}[ \t]*\n)+\n{,3})+',

501 """\

502 Brillouin zone point: 1

503 weight= 0.074074

504 k =< 0.333 0.333 0.333> . <b1,b2,b3>

505 =< 0.307 0.307 0.307>

506

507 orbital energies:

508 0.3919370E+00 ( 10.665eV) occ=1.000

509 0.3908827E+00 ( 10.637eV) occ=1.000 0.4155535E+00 ( 11.308eV) occ=1.000

510 0.3607689E+00 ( 9.817eV) occ=1.000 0.3827820E+00 ( 10.416eV) occ=1.000

511 0.3544000E+00 ( 9.644eV) occ=1.000 0.3782641E+00 ( 10.293eV) occ=1.000

512 0.3531137E+00 ( 9.609eV) occ=1.000 0.3778819E+00 ( 10.283eV) occ=1.000

513 0.2596367E+00 ( 7.065eV) occ=1.000 0.2820723E+00 ( 7.676eV) occ=1.000

514

515 Brillouin zone point: 2

516 weight= 0.074074

517 k =< -0.000 0.333 0.333> . <b1,b2,b3>

518 =< 0.614 0.000 0.000>

519

520 orbital energies:

521 0.3967132E+00 ( 10.795eV) occ=1.000

522 0.3920006E+00 ( 10.667eV) occ=1.000 0.4197952E+00 ( 11.423eV) occ=1.000

523 0.3912442E+00 ( 10.646eV) occ=1.000 0.4125086E+00 ( 11.225eV) occ=1.000

524 0.3910472E+00 ( 10.641eV) occ=1.000 0.4124238E+00 ( 11.223eV) occ=1.000

525 0.3153977E+00 ( 8.582eV) occ=1.000 0.3379797E+00 ( 9.197eV) occ=1.000

526 0.2801606E+00 ( 7.624eV) occ=1.000 0.3052478E+00 ( 8.306eV) occ=1.000

527""", re.M) # noqa: E501, W291

528

529# Parser for kpoint weights for a plane wave calculation

530_kpt_weight = _define_pattern(

531 r'^[ \t]+Brillouin zone point:[ \t]+([\S]+)[ \t]*\n'

532 r'^[ \t]+weight=[ \t]+([\S]+)[ \t]*\n',

533 """\

534 Brillouin zone point: 1

535 weight= 0.074074

536""", re.M) # noqa: W291

537

538

539# Parse eigenvalues and occupancies from a plane wave calculation

540def _get_pw_kpts(chunk):

541 eval_blocks = []

542 for block in _nwpw_eval_block.findall(chunk):

543 if 'pathlength' not in block:

544 eval_blocks.append(block)

545 if not eval_blocks:

546 return []

547 if 'virtual' in eval_blocks[-1]:

548 occ_block = eval_blocks[-2]

549 virt_block = eval_blocks[-1]

550 else:

551 occ_block = eval_blocks[-1]

552 virt_block = ''

553 kpts = NWChemKpts()

554 _extract_pw_kpts(occ_block, kpts, 1.)

555 _extract_pw_kpts(virt_block, kpts, 0.)

556 for match in _kpt_weight.finditer(occ_block):

557 index, weight = match.groups()

558 kpts.set_weight(index, float(weight))

559 return kpts.to_ibz_kpts(), kpts.to_singlepointkpts()

560

561

562# Helper class for keeping track of kpoints and converting to

563# SinglePointKPoint objects.

564class NWChemKpts:

565 def __init__(self):

566 self.data = {}

567 self.ibz_kpts = {}

568 self.weights = {}

569

570 def add_ibz_kpt(self, index, raw_kpt):

571 kpt = np.array([float(x.strip('>')) for x in raw_kpt.split()[1:4]])

572 self.ibz_kpts[index] = kpt

573

574 def add_eval(self, index, spin, energy, occ):

575 if index not in self.data:

576 self.data[index] = {}

577 if spin not in self.data[index]:

578 self.data[index][spin] = []

579 self.data[index][spin].append((energy, occ))

580

581 def set_weight(self, index, weight):

582 self.weights[index] = weight

583

584 def to_ibz_kpts(self):

585 if not self.ibz_kpts:

586 return np.array([[0., 0., 0.]])

587 sorted_kpts = sorted(list(self.ibz_kpts.items()), key=lambda x: x[0])

588 return np.array(list(zip(*sorted_kpts))[1])

589

590 def to_singlepointkpts(self):

591 kpts = []

592 for i, (index, spins) in enumerate(self.data.items()):

593 weight = self.weights[index]

594 for spin, (_, data) in enumerate(spins.items()):

595 energies, occs = np.array(sorted(data, key=lambda x: x[0])).T

596 kpts.append(SinglePointKPoint(weight, spin, i, energies, occs))

597 return kpts

598

599

600# Extracts MO/band data from a pattern matched by _nwpw_eval_block above

601_kpt = _define_pattern(

602 r'^[ \t]+Brillouin zone point:[ \t]+([\S]+)[ \t]*\n'

603 r'^[ \t]+weight=[ \t]+([\S])+[ \t]*\n'

604 r'^[ \t]+k[ \t]+([ \t\S]+)\n'

605 r'(?:^[ \t\S]*\n){1,2}'

606 r'^[ \t]+(?:virtual )?orbital energies:\n'

607 r'((?:^(?:(?:[ \t]+[\S]+){3,4}){1,2}[ \t]*\n)+)',

608 """\

609 Brillouin zone point: 1

610 weight= 0.074074

611 k =< 0.333 0.333 0.333> . <b1,b2,b3>

612 =< 0.307 0.307 0.307>

613

614 orbital energies:

615 0.3919370E+00 ( 10.665eV) occ=1.000

616 0.3908827E+00 ( 10.637eV) occ=1.000 0.4155535E+00 ( 11.308eV) occ=1.000

617 0.3607689E+00 ( 9.817eV) occ=1.000 0.3827820E+00 ( 10.416eV) occ=1.000

618 0.3544000E+00 ( 9.644eV) occ=1.000 0.3782641E+00 ( 10.293eV) occ=1.000

619 0.3531137E+00 ( 9.609eV) occ=1.000 0.3778819E+00 ( 10.283eV) occ=1.000

620 0.2596367E+00 ( 7.065eV) occ=1.000 0.2820723E+00 ( 7.676eV) occ=1.000

621""", re.M) # noqa: E501, W291

622

623

624# Extracts kpoints from a plane wave calculation

625def _extract_pw_kpts(chunk, kpts, default_occ):

626 for match in _kpt.finditer(chunk):

627 point, weight, raw_kpt, orbitals = match.groups()

628 index = int(point) - 1

629 for line in orbitals.split('\n'):

630 tokens = line.strip().split()

631 if not tokens:

632 continue

633 ntokens = len(tokens)

634 a_e = float(tokens[0]) * Hartree

635 if ntokens % 3 == 0:

636 a_o = default_occ

637 else:

638 a_o = float(tokens[3].split('=')[1])

639 kpts.add_eval(index, 0, a_e, a_o)

640

641 if ntokens <= 4:

642 continue

643 if ntokens == 6:

644 b_e = float(tokens[3]) * Hartree

645 b_o = default_occ

646 elif ntokens == 8:

647 b_e = float(tokens[4]) * Hartree

648 b_o = float(tokens[7].split('=')[1])

649 kpts.add_eval(index, 1, b_e, b_o)

650 kpts.set_weight(index, float(weight))

651 kpts.add_ibz_kpt(index, raw_kpt)