Coverage for /builds/kinetik161/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 else:

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

49 'output file.')

51 # First, find each SCF block

52 group = []

53 atomslist = []

54 header = True

55 lastgroup = []

56 lastparser = None

57 parser = None

58 for line in lines:

59 group.append(line)

60 if _gauss_block.match(line):

61 next_parser = parse_gto_chunk

62 elif _pw_block.match(line):

63 next_parser = parse_pw_chunk

64 else:

65 continue

67 if header:

68 header = False

69 else:

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

71 if atoms is None and parser is lastparser:

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

73 if atoms is not None:

74 atomslist[-1] = atoms

75 lastgroup += group

76 else:

77 atomslist.append(atoms)

78 lastgroup = group

79 lastparser = parser

80 group = []

81 parser = next_parser

82 else:

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 quadrupole = None

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

217 matches = pattern.findall(chunk)

218 if matches:

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

220 break

221

222 gradblocks = _gto_grad.findall(chunk)

223 if gradblocks:

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

225 natoms = len(gradblock)

226 symbols = []

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

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

229 for i, line in enumerate(gradblock):

230 line = line.strip().split()

231 symbols.append(line[1])

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

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

234 pos *= Bohr

235 forces *= Hartree / Bohr

236 atoms = Atoms(symbols, positions=pos)

237

238 dipole, quadrupole = _get_multipole(chunk)

239

240 kpts = _get_gto_kpts(chunk)

241

242 if atoms is None:

243 atoms = _parse_geomblock(chunk)

244

245 if atoms is None:

246 return

247

248 # SinglePointDFTCalculator doesn't support quadrupole moment currently

249 calc = SinglePointDFTCalculator(atoms=atoms,

250 energy=energy,

251 free_energy=energy, # XXX Is this right?

252 forces=forces,

253 dipole=dipole,

254 # quadrupole=quadrupole,

255 )

256 calc.kpts = kpts

257 atoms.calc = calc

258 return atoms

259

260

261# Extracts dipole and quadrupole moment for a GTO calculation

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

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

264# in between \n and \n

265_multipole = _define_pattern(

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

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

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

269 """\

270 Multipole analysis of the density

271 ---------------------------------

272

273 L x y z total alpha beta nuclear

274 - - - - ----- ----- ---- -------

275 0 0 0 0 -0.000000 -5.000000 -5.000000 10.000000

276

277 1 1 0 0 0.000000 0.000000 0.000000 0.000000

278 1 0 1 0 -0.000001 -0.000017 -0.000017 0.000034

279 1 0 0 1 -0.902084 -0.559881 -0.559881 0.217679

280

281 2 2 0 0 -5.142958 -2.571479 -2.571479 0.000000

282 2 1 1 0 -0.000000 -0.000000 -0.000000 0.000000

283 2 1 0 1 0.000000 0.000000 0.000000 0.000000

284 2 0 2 0 -3.153324 -3.807308 -3.807308 4.461291

285 2 0 1 1 0.000001 -0.000009 -0.000009 0.000020

286 2 0 0 2 -4.384288 -3.296205 -3.296205 2.208122

287""", re.M)

288

289

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

291def _get_multipole(chunk):

292 matches = _multipole.findall(chunk)

293 if not matches:

294 return None, None

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

296 # this column contains the actual moments.

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

298 if x and not x.isspace()]

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

300 quadrupole = np.zeros(9)

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

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

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

304

305

306# MO eigenvalue and occupancy parser for GTO calculations

307_eval_block = _define_pattern(

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

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

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

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

312 """\

313 ROHF Final Molecular Orbital Analysis

314 -------------------------------------

315

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

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

318 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function

319 ----- ------------ --------------- ----- ------------ ---------------

320 1 0.983233 1 O s

321

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

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

324 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function

325 ----- ------------ --------------- ----- ------------ ---------------

326 6 0.708998 1 O s 1 -0.229426 1 O s

327 2 0.217752 1 O s

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

329

330

331# Parses the eigenvalues and occupations from a GTO calculation

332def _get_gto_kpts(chunk):

333 eval_blocks = _eval_block.findall(chunk)

334 if not eval_blocks:

335 return []

336 kpts = []

337 kpt = _get_gto_evals(eval_blocks[-1])

338 if kpt.s == 1:

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

340 kpts.append(kpt)

341 return kpts

342

343

344# Extracts MO eigenvalue and occupancy for a GTO calculation

345_extract_vector = _define_pattern(

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

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

348)

349

350

351# Extracts the eigenvalues and occupations from a GTO calculation

352def _get_gto_evals(chunk):

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

354 data = []

355 for vector in _extract_vector.finditer(chunk):

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

357 data = np.array(data)

358 occ = data[:, 0]

359 energies = data[:, 1] * Hartree

360

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

362

363

364# Plane wave specific parsing stuff

365

366# Matches the gradient block from a plane wave calculation

367_nwpw_grad = _define_pattern(

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

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

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

371 """\

372 ============= Ion Gradients =================

373 Ion Forces:

374 1 O ( -0.000012 0.000027 -0.005199 )

375 2 H ( 0.000047 -0.013082 0.020790 )

376 3 H ( 0.000047 0.012863 0.020786 )

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

378 ===============================================

379""", re.M)

380

381# Matches the gradient block from a PAW calculation

382_paw_grad = _define_pattern(

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

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

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

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

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

388 """\

389 ============= Ion Gradients =================

390 Ion Positions:

391 1 O ( -3.77945 -5.22176 -3.77945 )

392 2 H ( -3.77945 -3.77945 3.77945 )

393 3 H ( -3.77945 3.77945 3.77945 )

394 Ion Forces:

395 1 O ( -0.00001 -0.00000 0.00081 )

396 2 H ( 0.00005 -0.00026 -0.00322 )

397 3 H ( 0.00005 0.00030 -0.00322 )

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

399 ===============================================

400""", re.M)

401

402# Energy parser for plane wave calculations

403_nwpw_energy = _define_pattern(

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

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

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

407 re.M,

408)

409

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

411_fermi_energy = _define_pattern(

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

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

414)

415

416

417# Plane wave parser

418def parse_pw_chunk(chunk):

419 atoms = _parse_geomblock(chunk)

420 if atoms is None:

421 return

422

423 energy = None

424 efermi = None

425 forces = None

426 stress = None

427

428 matches = _nwpw_energy.findall(chunk)

429 if matches:

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

431

432 matches = _fermi_energy.findall(chunk)

433 if matches:

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

435

436 gradblocks = _nwpw_grad.findall(chunk)

437 if not gradblocks:

438 gradblocks = _paw_grad.findall(chunk)

439 if gradblocks:

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

441 natoms = len(gradblock)

442 symbols = []

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

444 for i, line in enumerate(gradblock):

445 line = line.strip().split()

446 symbols.append(line[1])

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

448 forces *= Hartree / Bohr

449

450 if atoms.cell:

451 stress = _get_stress(chunk, atoms.cell)

452

453 ibz_kpts, kpts = _get_pw_kpts(chunk)

454

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

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

457 calc = SinglePointDFTCalculator(atoms=atoms,

458 energy=energy,

459 efermi=efermi,

460 free_energy=energy,

461 forces=forces,

462 stress=stress,

463 ibzkpts=ibz_kpts)

464 calc.kpts = kpts

465 atoms.calc = calc

466 return atoms

467

468

469# Extracts stress tensor from a plane wave calculation

470_stress = _define_pattern(

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

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

473 """\

474 ============= total gradient ==============

475 S = ( -0.22668 0.27174 0.19134 )

476 ( 0.23150 -0.26760 0.23226 )

477 ( 0.19090 0.27206 -0.22700 )

478 ===================================================

479""", re.M)

480

481

482# Extract stress tensor from a plane wave calculation

483def _get_stress(chunk, cell):

484 stress_blocks = _stress.findall(chunk)

485 if not stress_blocks:

486 return None

487 stress_block = stress_blocks[-1]

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

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

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

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

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

493 # convert from 3x3 array to Voigt form

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

495

496

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

498_nwpw_eval_block = _define_pattern(

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

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

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

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

503 """\

504 Brillouin zone point: 1

505 weight= 0.074074

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

507 =< 0.307 0.307 0.307>

508

509 orbital energies:

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

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

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

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

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

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

516

517 Brillouin zone point: 2

518 weight= 0.074074

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

520 =< 0.614 0.000 0.000>

521

522 orbital energies:

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

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

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

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

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

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

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

530

531# Parser for kpoint weights for a plane wave calculation

532_kpt_weight = _define_pattern(

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

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

535 """\

536 Brillouin zone point: 1

537 weight= 0.074074

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

539

540

541# Parse eigenvalues and occupancies from a plane wave calculation

542def _get_pw_kpts(chunk):

543 eval_blocks = []

544 for block in _nwpw_eval_block.findall(chunk):

545 if 'pathlength' not in block:

546 eval_blocks.append(block)

547 if not eval_blocks:

548 return []

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

550 occ_block = eval_blocks[-2]

551 virt_block = eval_blocks[-1]

552 else:

553 occ_block = eval_blocks[-1]

554 virt_block = ''

555 kpts = NWChemKpts()

556 _extract_pw_kpts(occ_block, kpts, 1.)

557 _extract_pw_kpts(virt_block, kpts, 0.)

558 for match in _kpt_weight.finditer(occ_block):

559 index, weight = match.groups()

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

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

562

563

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

565# SinglePointKPoint objects.

566class NWChemKpts:

567 def __init__(self):

568 self.data = {}

569 self.ibz_kpts = {}

570 self.weights = {}

571

572 def add_ibz_kpt(self, index, raw_kpt):

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

574 self.ibz_kpts[index] = kpt

575

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

577 if index not in self.data:

578 self.data[index] = {}

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

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

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

582

583 def set_weight(self, index, weight):

584 self.weights[index] = weight

585

586 def to_ibz_kpts(self):

587 if not self.ibz_kpts:

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

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

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

591

592 def to_singlepointkpts(self):

593 kpts = []

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

595 weight = self.weights[index]

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

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

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

599 return kpts

600

601

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

603_kpt = _define_pattern(

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

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

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

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

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

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

610 """\

611 Brillouin zone point: 1

612 weight= 0.074074

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

614 =< 0.307 0.307 0.307>

615

616 orbital energies:

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

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

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

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

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

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

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

624

625

626# Extracts kpoints from a plane wave calculation

627def _extract_pw_kpts(chunk, kpts, default_occ):

628 for match in _kpt.finditer(chunk):

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

630 index = int(point) - 1

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

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

633 if not tokens:

634 continue

635 ntokens = len(tokens)

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

637 if ntokens % 3 == 0:

638 a_o = default_occ

639 else:

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

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

642

643 if ntokens <= 4:

644 continue

645 if ntokens == 6:

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

647 b_o = default_occ

648 elif ntokens == 8:

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

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

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

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

653 kpts.add_ibz_kpt(index, raw_kpt)