Coverage for /builds/debichem-team/python-ase/ase/io/nwchem/nwreader.py: 65.44%

1import re

2from collections import OrderedDict

4import numpy as np

6from ase import Atoms

7from ase.calculators.singlepoint import (

8 SinglePointDFTCalculator,

9 SinglePointKPoint,

10)

11from ase.units import Bohr, Hartree

13from .parser import _define_pattern

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

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

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

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

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

20# the regular expressions are still working correctly.

22# Matches the beginning of a GTO calculation

23_gauss_block = _define_pattern(

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

25 " NWChem SCF Module\n",

26)

29# Matches the beginning of a plane wave calculation

30_pw_block = _define_pattern(

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

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

33 " * NWPW PSPW Calculation *\n",

34)

37# Top-level parser

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

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

40 individual single point calculations."""

41 lines = fobj.readlines()

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

44 for line in lines:

45 if _gauss_block.match(line):

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

47 if _pw_block.match(line):

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

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

50 'output file.')

52 # First, find each SCF block

53 group = []

54 atomslist = []

55 header = True

56 lastgroup = []

57 lastparser = None

58 parser = None

59 for line in lines:

60 group.append(line)

61 if _gauss_block.match(line):

62 next_parser = parse_gto_chunk

63 elif _pw_block.match(line):

64 next_parser = parse_pw_chunk

65 else:

66 continue

68 if header:

69 header = False

70 else:

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

72 if atoms is None and parser is lastparser:

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

74 if atoms is not None:

75 atomslist[-1] = atoms

76 lastgroup += group

77 else:

78 atomslist.append(atoms)

79 lastgroup = group

80 lastparser = parser

81 group = []

82 parser = next_parser

83 if not header:

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

85 if atoms is not None:

86 atomslist.append(atoms)

88 return atomslist[index]

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

92_geom = _define_pattern(

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

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

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

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

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

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

99 """\

100

101 Geometry "geometry" -> ""

102 -------------------------

103

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

105

106 No. Tag Charge X Y Z

107 ---- ---------------- ---------- -------------- -------------- --------------

108 1 C 6.0000 0.00000000 0.00000000 0.00000000

109 2 H 1.0000 0.62911800 0.62911800 0.62911800

110 3 H 1.0000 -0.62911800 -0.62911800 0.62911800

111 4 H 1.0000 0.62911800 -0.62911800 -0.62911800

112""", re.M)

113

114# Unit cell parser

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

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

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

118 """\

119 Lattice Parameters

120 ------------------

121

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

123

124 a1=< 4.000 0.000 0.000 >

125 a2=< 0.000 5.526 0.000 >

126 a3=< 0.000 0.000 4.596 >

127 a= 4.000 b= 5.526 c= 4.596

128 alpha= 90.000 beta= 90.000 gamma= 90.000

129 omega= 101.6

130""", re.M)

131

132

133# Parses the geometry and returns the corresponding Atoms object

134def _parse_geomblock(chunk):

135 geomblocks = _geom.findall(chunk)

136 if not geomblocks:

137 return None

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

139 natoms = len(geomblock)

140 symbols = []

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

142 for i, line in enumerate(geomblock):

143 line = line.strip().split()

144 symbols.append(line[1])

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

146

147 cellblocks = _cell_block.findall(chunk)

148 if cellblocks:

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

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

151 for i, line in enumerate(cellblock):

152 line = line.strip().split()

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

154 else:

155 cell = None

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

157

158

159# GTO-specific parser stuff

160

161# Matches gradient block from a GTO calculation

162_gto_grad = _define_pattern(

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

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

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

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

167 """\

168 UHF ENERGY GRADIENTS

169

170 atom coordinates gradient

171 x y z x y z

172 1 C 0.293457 -0.293457 0.293457 -0.000083 0.000083 -0.000083

173 2 H 1.125380 1.355351 1.125380 0.000086 0.000089 0.000086

174 3 H -1.355351 -1.125380 1.125380 -0.000089 -0.000086 0.000086

175 4 H 1.125380 -1.125380 -1.355351 0.000086 -0.000086 -0.000089

176

177""", re.M)

178

179# Energy parsers for a variety of different GTO calculations

180_e_gto = OrderedDict()

181_e_gto['tce'] = _define_pattern(

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

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

184 " CCD total energy / hartree "

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

186)

187_e_gto['ccsd'] = _define_pattern(

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

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

190 re.M,

191)

192_e_gto['tddft'] = _define_pattern(

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

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

195 re.M,

196)

197_e_gto['mp2'] = _define_pattern(

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

199 " Total MP2 energy -75.708800087578\n",

200 re.M,

201)

202_e_gto['mf'] = _define_pattern(

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

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

205 re.M,

206)

207

208

209# GTO parser

210def parse_gto_chunk(chunk):

211 atoms = None

212 forces = None

213 energy = None

214 dipole = None

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

216 matches = pattern.findall(chunk)

217 if matches:

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

219 break

220

221 gradblocks = _gto_grad.findall(chunk)

222 if gradblocks:

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

224 natoms = len(gradblock)

225 symbols = []

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

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

228 for i, line in enumerate(gradblock):

229 line = line.strip().split()

230 symbols.append(line[1])

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

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

233 pos *= Bohr

234 forces *= Hartree / Bohr

235 atoms = Atoms(symbols, positions=pos)

236

237 dipole, _quadrupole = _get_multipole(chunk)

238

239 kpts = _get_gto_kpts(chunk)

240

241 if atoms is None:

242 atoms = _parse_geomblock(chunk)

243

244 if atoms is None:

245 return None

246

247 # SinglePointDFTCalculator doesn't support quadrupole moment currently

248 calc = SinglePointDFTCalculator(atoms=atoms,

249 energy=energy,

250 free_energy=energy, # XXX Is this right?

251 forces=forces,

252 dipole=dipole,

253 # quadrupole=quadrupole,

254 )

255 calc.kpts = kpts

256 atoms.calc = calc

257 return atoms

258

259

260# Extracts dipole and quadrupole moment for a GTO calculation

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

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

263# in between \n and \n

264_multipole = _define_pattern(

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

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

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

268 """\

269 Multipole analysis of the density

270 ---------------------------------

271

272 L x y z total alpha beta nuclear

273 - - - - ----- ----- ---- -------

274 0 0 0 0 -0.000000 -5.000000 -5.000000 10.000000

275

276 1 1 0 0 0.000000 0.000000 0.000000 0.000000

277 1 0 1 0 -0.000001 -0.000017 -0.000017 0.000034

278 1 0 0 1 -0.902084 -0.559881 -0.559881 0.217679

279

280 2 2 0 0 -5.142958 -2.571479 -2.571479 0.000000

281 2 1 1 0 -0.000000 -0.000000 -0.000000 0.000000

282 2 1 0 1 0.000000 0.000000 0.000000 0.000000

283 2 0 2 0 -3.153324 -3.807308 -3.807308 4.461291

284 2 0 1 1 0.000001 -0.000009 -0.000009 0.000020

285 2 0 0 2 -4.384288 -3.296205 -3.296205 2.208122

286""", re.M)

287

288

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

290def _get_multipole(chunk):

291 matches = _multipole.findall(chunk)

292 if not matches:

293 return None, None

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

295 # this column contains the actual moments.

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

297 if x and not x.isspace()]

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

299 quadrupole = np.zeros(9)

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

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

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

303

304

305# MO eigenvalue and occupancy parser for GTO calculations

306_eval_block = _define_pattern(

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

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

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

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

311 """\

312 ROHF Final Molecular Orbital Analysis

313 -------------------------------------

314

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

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

317 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function

318 ----- ------------ --------------- ----- ------------ ---------------

319 1 0.983233 1 O s

320

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

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

323 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function

324 ----- ------------ --------------- ----- ------------ ---------------

325 6 0.708998 1 O s 1 -0.229426 1 O s

326 2 0.217752 1 O s

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

328

329

330# Parses the eigenvalues and occupations from a GTO calculation

331def _get_gto_kpts(chunk):

332 eval_blocks = _eval_block.findall(chunk)

333 if not eval_blocks:

334 return []

335 kpts = []

336 kpt = _get_gto_evals(eval_blocks[-1])

337 if kpt.s == 1:

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

339 kpts.append(kpt)

340 return kpts

341

342

343# Extracts MO eigenvalue and occupancy for a GTO calculation

344_extract_vector = _define_pattern(

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

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

347)

348

349

350# Extracts the eigenvalues and occupations from a GTO calculation

351def _get_gto_evals(chunk):

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

353 data = []

354 for vector in _extract_vector.finditer(chunk):

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

356 data = np.array(data)

357 occ = data[:, 0]

358 energies = data[:, 1] * Hartree

359

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

361

362

363# Plane wave specific parsing stuff

364

365# Matches the gradient block from a plane wave calculation

366_nwpw_grad = _define_pattern(

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

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

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

370 """\

371 ============= Ion Gradients =================

372 Ion Forces:

373 1 O ( -0.000012 0.000027 -0.005199 )

374 2 H ( 0.000047 -0.013082 0.020790 )

375 3 H ( 0.000047 0.012863 0.020786 )

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

377 ===============================================

378""", re.M)

379

380# Matches the gradient block from a PAW calculation

381_paw_grad = _define_pattern(

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

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

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

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

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

387 """\

388 ============= Ion Gradients =================

389 Ion Positions:

390 1 O ( -3.77945 -5.22176 -3.77945 )

391 2 H ( -3.77945 -3.77945 3.77945 )

392 3 H ( -3.77945 3.77945 3.77945 )

393 Ion Forces:

394 1 O ( -0.00001 -0.00000 0.00081 )

395 2 H ( 0.00005 -0.00026 -0.00322 )

396 3 H ( 0.00005 0.00030 -0.00322 )

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

398 ===============================================

399""", re.M)

400

401# Energy parser for plane wave calculations

402_nwpw_energy = _define_pattern(

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

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

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

406 re.M,

407)

408

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

410_fermi_energy = _define_pattern(

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

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

413)

414

415

416# Plane wave parser

417def parse_pw_chunk(chunk):

418 atoms = _parse_geomblock(chunk)

419 if atoms is None:

420 return None

421

422 energy = None

423 efermi = None

424 forces = None

425 stress = None

426

427 matches = _nwpw_energy.findall(chunk)

428 if matches:

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

430

431 matches = _fermi_energy.findall(chunk)

432 if matches:

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

434

435 gradblocks = _nwpw_grad.findall(chunk)

436 if not gradblocks:

437 gradblocks = _paw_grad.findall(chunk)

438 if gradblocks:

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

440 natoms = len(gradblock)

441 symbols = []

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

443 for i, line in enumerate(gradblock):

444 line = line.strip().split()

445 symbols.append(line[1])

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

447 forces *= Hartree / Bohr

448

449 if atoms.cell:

450 stress = _get_stress(chunk, atoms.cell)

451

452 ibz_kpts, kpts = _get_pw_kpts(chunk)

453

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

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

456 calc = SinglePointDFTCalculator(atoms=atoms,

457 energy=energy,

458 efermi=efermi,

459 free_energy=energy,

460 forces=forces,

461 stress=stress,

462 ibzkpts=ibz_kpts)

463 calc.kpts = kpts

464 atoms.calc = calc

465 return atoms

466

467

468# Extracts stress tensor from a plane wave calculation

469_stress = _define_pattern(

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

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

472 """\

473 ============= total gradient ==============

474 S = ( -0.22668 0.27174 0.19134 )

475 ( 0.23150 -0.26760 0.23226 )

476 ( 0.19090 0.27206 -0.22700 )

477 ===================================================

478""", re.M)

479

480

481# Extract stress tensor from a plane wave calculation

482def _get_stress(chunk, cell):

483 stress_blocks = _stress.findall(chunk)

484 if not stress_blocks:

485 return None

486 stress_block = stress_blocks[-1]

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

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

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

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

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

492 # convert from 3x3 array to Voigt form

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

494

495

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

497_nwpw_eval_block = _define_pattern(

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

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

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

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

502 """\

503 Brillouin zone point: 1

504 weight= 0.074074

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

506 =< 0.307 0.307 0.307>

507

508 orbital energies:

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

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

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

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

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

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

515

516 Brillouin zone point: 2

517 weight= 0.074074

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

519 =< 0.614 0.000 0.000>

520

521 orbital energies:

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

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

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

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

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

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

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

529

530# Parser for kpoint weights for a plane wave calculation

531_kpt_weight = _define_pattern(

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

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

534 """\

535 Brillouin zone point: 1

536 weight= 0.074074

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

538

539

540# Parse eigenvalues and occupancies from a plane wave calculation

541def _get_pw_kpts(chunk):

542 eval_blocks = []

543 for block in _nwpw_eval_block.findall(chunk):

544 if 'pathlength' not in block:

545 eval_blocks.append(block)

546 if not eval_blocks:

547 return []

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

549 occ_block = eval_blocks[-2]

550 virt_block = eval_blocks[-1]

551 else:

552 occ_block = eval_blocks[-1]

553 virt_block = ''

554 kpts = NWChemKpts()

555 _extract_pw_kpts(occ_block, kpts, 1.)

556 _extract_pw_kpts(virt_block, kpts, 0.)

557 for match in _kpt_weight.finditer(occ_block):

558 index, weight = match.groups()

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

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

561

562

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

564# SinglePointKPoint objects.

565class NWChemKpts:

566 def __init__(self):

567 self.data = {}

568 self.ibz_kpts = {}

569 self.weights = {}

570

571 def add_ibz_kpt(self, index, raw_kpt):

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

573 self.ibz_kpts[index] = kpt

574

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

576 if index not in self.data:

577 self.data[index] = {}

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

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

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

581

582 def set_weight(self, index, weight):

583 self.weights[index] = weight

584

585 def to_ibz_kpts(self):

586 if not self.ibz_kpts:

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

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

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

590

591 def to_singlepointkpts(self):

592 kpts = []

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

594 weight = self.weights[index]

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

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

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

598 return kpts

599

600

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

602_kpt = _define_pattern(

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

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

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

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

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

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

609 """\

610 Brillouin zone point: 1

611 weight= 0.074074

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

613 =< 0.307 0.307 0.307>

614

615 orbital energies:

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

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

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

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

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

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

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

623

624

625# Extracts kpoints from a plane wave calculation

626def _extract_pw_kpts(chunk, kpts, default_occ):

627 for match in _kpt.finditer(chunk):

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

629 index = int(point) - 1

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

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

632 if not tokens:

633 continue

634 ntokens = len(tokens)

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

636 if ntokens % 3 == 0:

637 a_o = default_occ

638 else:

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

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

641

642 if ntokens <= 4:

643 continue

644 if ntokens == 6:

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

646 b_o = default_occ

647 elif ntokens == 8:

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

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

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

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

652 kpts.add_ibz_kpt(index, raw_kpt)