Coverage for /builds/kinetik161/ase/ase/io/abinit.py: 85.36%

1import os

2import re

3from glob import glob

4from os.path import join

5from pathlib import Path

7import numpy as np

9from ase import Atoms

10from ase.calculators.calculator import all_properties

11from ase.calculators.singlepoint import SinglePointCalculator

12from ase.data import chemical_symbols

13from ase.units import Bohr, Hartree, fs

16def read_abinit_in(fd):

17 """Import ABINIT input file.

19 Reads cell, atom positions, etc. from abinit input file

20 """

22 tokens = []

24 for line in fd:

25 meat = line.split('#', 1)[0]

26 tokens += meat.lower().split()

28 # note that the file can not be scanned sequentially

30 index = tokens.index("acell")

31 unit = 1.0

32 if tokens[index + 4].lower()[:3] != 'ang':

33 unit = Bohr

34 acell = [unit * float(tokens[index + 1]),

35 unit * float(tokens[index + 2]),

36 unit * float(tokens[index + 3])]

38 index = tokens.index("natom")

39 natom = int(tokens[index + 1])

41 index = tokens.index("ntypat")

42 ntypat = int(tokens[index + 1])

44 index = tokens.index("typat")

45 typat = []

46 while len(typat) < natom:

47 token = tokens[index + 1]

48 if '*' in token: # e.g. typat 4*1 3*2 ...

49 nrepeat, typenum = token.split('*')

50 typat += [int(typenum)] * int(nrepeat)

51 else:

52 typat.append(int(token))

53 index += 1

54 assert natom == len(typat)

56 index = tokens.index("znucl")

57 znucl = []

58 for i in range(ntypat):

59 znucl.append(int(tokens[index + 1 + i]))

61 index = tokens.index("rprim")

62 rprim = []

63 for i in range(3):

64 rprim.append([acell[i] * float(tokens[index + 3 * i + 1]),

65 acell[i] * float(tokens[index + 3 * i + 2]),

66 acell[i] * float(tokens[index + 3 * i + 3])])

68 # create a list with the atomic numbers

69 numbers = []

70 for i in range(natom):

71 ii = typat[i] - 1

72 numbers.append(znucl[ii])

74 # now the positions of the atoms

75 if "xred" in tokens:

76 index = tokens.index("xred")

77 xred = []

78 for i in range(natom):

79 xred.append([float(tokens[index + 3 * i + 1]),

80 float(tokens[index + 3 * i + 2]),

81 float(tokens[index + 3 * i + 3])])

82 atoms = Atoms(cell=rprim, scaled_positions=xred, numbers=numbers,

83 pbc=True)

84 else:

85 if "xcart" in tokens:

86 index = tokens.index("xcart")

87 unit = Bohr

88 elif "xangst" in tokens:

89 unit = 1.0

90 index = tokens.index("xangst")

91 else:

92 raise OSError(

93 "No xred, xcart, or xangs keyword in abinit input file")

95 xangs = []

96 for i in range(natom):

97 xangs.append([unit * float(tokens[index + 3 * i + 1]),

98 unit * float(tokens[index + 3 * i + 2]),

99 unit * float(tokens[index + 3 * i + 3])])

100 atoms = Atoms(cell=rprim, positions=xangs, numbers=numbers, pbc=True)

101

102 try:

103 ii = tokens.index('nsppol')

104 except ValueError:

105 nsppol = None

106 else:

107 nsppol = int(tokens[ii + 1])

108

109 if nsppol == 2:

110 index = tokens.index('spinat')

111 magmoms = [float(tokens[index + 3 * i + 3]) for i in range(natom)]

112 atoms.set_initial_magnetic_moments(magmoms)

113

114 assert len(atoms) == natom

115 return atoms

116

117

118keys_with_units = {

119 'toldfe': 'eV',

120 'tsmear': 'eV',

121 'paoenergyshift': 'eV',

122 'zmunitslength': 'Bohr',

123 'zmunitsangle': 'rad',

124 'zmforcetollength': 'eV/Ang',

125 'zmforcetolangle': 'eV/rad',

126 'zmmaxdispllength': 'Ang',

127 'zmmaxdisplangle': 'rad',

128 'ecut': 'eV',

129 'pawecutdg': 'eV',

130 'dmenergytolerance': 'eV',

131 'electronictemperature': 'eV',

132 'oneta': 'eV',

133 'onetaalpha': 'eV',

134 'onetabeta': 'eV',

135 'onrclwf': 'Ang',

136 'onchemicalpotentialrc': 'Ang',

137 'onchemicalpotentialtemperature': 'eV',

138 'mdmaxcgdispl': 'Ang',

139 'mdmaxforcetol': 'eV/Ang',

140 'mdmaxstresstol': 'eV/Ang**3',

141 'mdlengthtimestep': 'fs',

142 'mdinitialtemperature': 'eV',

143 'mdtargettemperature': 'eV',

144 'mdtargetpressure': 'eV/Ang**3',

145 'mdnosemass': 'eV*fs**2',

146 'mdparrinellorahmanmass': 'eV*fs**2',

147 'mdtaurelax': 'fs',

148 'mdbulkmodulus': 'eV/Ang**3',

149 'mdfcdispl': 'Ang',

150 'warningminimumatomicdistance': 'Ang',

151 'rcspatial': 'Ang',

152 'kgridcutoff': 'Ang',

153 'latticeconstant': 'Ang'}

154

155

156def write_abinit_in(fd, atoms, param=None, species=None, pseudos=None):

157 from ase.calculators.calculator import kpts2mp

158

159 if param is None:

160 param = {}

161

162 if species is None:

163 species = sorted(set(atoms.numbers))

164

165 inp = dict(param)

166 xc = inp.pop('xc', 'LDA')

167 for key in ['smearing', 'kpts', 'pps', 'raw']:

168 inp.pop(key, None)

169

170 smearing = param.get('smearing')

171 if 'tsmear' in param or 'occopt' in param:

172 assert smearing is None

173

174 if smearing is not None:

175 inp['occopt'] = {'fermi-dirac': 3,

176 'gaussian': 7}[smearing[0].lower()]

177 inp['tsmear'] = smearing[1]

178

179 inp['natom'] = len(atoms)

180

181 if 'nbands' in param:

182 inp['nband'] = param['nbands']

183 del inp['nbands']

184

185 # ixc is set from paw/xml file. Ignore 'xc' setting then.

186 if param.get('pps') not in ['pawxml']:

187 if 'ixc' not in param:

188 inp['ixc'] = {'LDA': 7,

189 'PBE': 11,

190 'revPBE': 14,

191 'RPBE': 15,

192 'WC': 23}[xc]

193

194 magmoms = atoms.get_initial_magnetic_moments()

195 if magmoms.any():

196 inp['nsppol'] = 2

197 fd.write('spinat\n')

198 for n, M in enumerate(magmoms):

199 fd.write(f'{0:.14f} {0:.14f} {M:.14f}\n')

200 else:

201 inp['nsppol'] = 1

202

203 if param.get('kpts') is not None:

204 mp = kpts2mp(atoms, param['kpts'])

205 fd.write('kptopt 1\n')

206 fd.write('ngkpt %d %d %d\n' % tuple(mp))

207 fd.write('nshiftk 1\n')

208 fd.write('shiftk\n')

209 fd.write('%.1f %.1f %.1f\n' % tuple((np.array(mp) + 1) % 2 * 0.5))

210

211 valid_lists = (list, np.ndarray)

212 for key in sorted(inp):

213 value = inp[key]

214 unit = keys_with_units.get(key)

215 if unit is not None:

216 if 'fs**2' in unit:

217 value /= fs**2

218 elif 'fs' in unit:

219 value /= fs

220 if isinstance(value, valid_lists):

221 if isinstance(value[0], valid_lists):

222 fd.write(f"{key}\n")

223 for dim in value:

224 write_list(fd, dim, unit)

225 else:

226 fd.write(f"{key}\n")

227 write_list(fd, value, unit)

228 else:

229 if unit is None:

230 fd.write(f"{key} {value}\n")

231 else:

232 fd.write(f"{key} {value} {unit}\n")

233

234 if param.get('raw') is not None:

235 if isinstance(param['raw'], str):

236 raise TypeError('The raw parameter is a single string; expected '

237 'a sequence of lines')

238 for line in param['raw']:

239 if isinstance(line, tuple):

240 fd.write(' '.join([f'{x}' for x in line]) + '\n')

241 else:

242 fd.write(f'{line}\n')

243

244 fd.write('#Definition of the unit cell\n')

245 fd.write('acell\n')

246 fd.write(f'{1.0:.14f} {1.0:.14f} {1.0:.14f} Angstrom\n')

247 fd.write('rprim\n')

248 if atoms.cell.rank != 3:

249 raise RuntimeError('Abinit requires a 3D cell, but cell is {}'

250 .format(atoms.cell))

251 for v in atoms.cell:

252 fd.write('%.14f %.14f %.14f\n' % tuple(v))

253

254 fd.write('chkprim 0 # Allow non-primitive cells\n')

255

256 fd.write('#Definition of the atom types\n')

257 fd.write('ntypat %d\n' % (len(species)))

258 fd.write('znucl {}\n'.format(' '.join(str(Z) for Z in species)))

259 fd.write('#Enumerate different atomic species\n')

260 fd.write('typat')

261 fd.write('\n')

262

263 types = []

264 for Z in atoms.numbers:

265 for n, Zs in enumerate(species):

266 if Z == Zs:

267 types.append(n + 1)

268 n_entries_int = 20 # integer entries per line

269 for n, type in enumerate(types):

270 fd.write(' %d' % (type))

271 if n > 1 and ((n % n_entries_int) == 1):

272 fd.write('\n')

273 fd.write('\n')

274

275 if pseudos is not None:

276 listing = ',\n'.join(pseudos)

277 line = f'pseudos "{listing}"\n'

278 fd.write(line)

279

280 fd.write('#Definition of the atoms\n')

281 fd.write('xcart\n')

282 for pos in atoms.positions / Bohr:

283 fd.write('%.14f %.14f %.14f\n' % tuple(pos))

284

285 fd.write('chkexit 1 # abinit.exit file in the running '

286 'directory terminates after the current SCF\n')

287

288

289def write_list(fd, value, unit):

290 for element in value:

291 fd.write(f"{element} ")

292 if unit is not None:

293 fd.write(f"{unit}")

294 fd.write("\n")

295

296

297def read_stress(fd):

298 # sigma(1 1)= 4.02063464E-04 sigma(3 2)= 0.00000000E+00

299 # sigma(2 2)= 4.02063464E-04 sigma(3 1)= 0.00000000E+00

300 # sigma(3 3)= 4.02063464E-04 sigma(2 1)= 0.00000000E+00

301 pat = re.compile(r'\s*sigma\(\d \d\)=\s*(\S+)\s*sigma\(\d \d\)=\s*(\S+)')

302 stress = np.empty(6)

303 for i in range(3):

304 line = next(fd)

305 m = pat.match(line)

306 if m is None:

307 # Not a real value error. What should we raise?

308 raise ValueError('Line {!r} does not match stress pattern {!r}'

309 .format(line, pat))

310 s1, s2 = m.group(1, 2)

311 stress[i] = float(m.group(1))

312 stress[i + 3] = float(m.group(2))

313 unit = Hartree / Bohr**3

314 return stress * unit

315

316

317def consume_multiline(fd, headerline, nvalues, dtype):

318 """Parse abinit-formatted "header + values" sections.

319

320 Example:

321

322 typat 1 1 1 1 1

323 1 1 1 1

324 """

325 tokens = headerline.split()

326 assert tokens[0].isalpha()

327

328 values = tokens[1:]

329 while len(values) < nvalues:

330 line = next(fd)

331 values.extend(line.split())

332 assert len(values) == nvalues

333 return np.array(values).astype(dtype)

334

335

336def read_abinit_out(fd):

337 results = {}

338

339 def skipto(string):

340 for line in fd:

341 if string in line:

342 return line

343 raise RuntimeError(f'Not found: {string}')

344

345 line = skipto('Version')

346 m = re.match(r'\.*?Version\s+(\S+)\s+of ABINIT', line)

347 assert m is not None

348 version = m.group(1)

349 results['version'] = version

350

351 use_v9_format = int(version.split('.', 1)[0]) >= 9

352

353 shape_vars = {}

354

355 skipto('echo values of preprocessed input variables')

356

357 for line in fd:

358 if '===============' in line:

359 break

360

361 tokens = line.split()

362 if not tokens:

363 continue

364

365 for key in ['natom', 'nkpt', 'nband', 'ntypat']:

366 if tokens[0] == key:

367 shape_vars[key] = int(tokens[1])

368

369 if line.lstrip().startswith('typat'): # Avoid matching ntypat

370 types = consume_multiline(fd, line, shape_vars['natom'], int)

371

372 if 'znucl' in line:

373 znucl = consume_multiline(fd, line, shape_vars['ntypat'], float)

374

375 if 'rprim' in line:

376 cell = consume_multiline(fd, line, 9, float)

377 cell = cell.reshape(3, 3)

378

379 natoms = shape_vars['natom']

380

381 # Skip ahead to results:

382 for line in fd:

383 if 'was not enough scf cycles to converge' in line:

384 raise RuntimeError(line)

385 if 'iterations are completed or convergence reached' in line:

386 break

387 else:

388 raise RuntimeError('Cannot find results section')

389

390 def read_array(fd, nlines):

391 arr = []

392 for i in range(nlines):

393 line = next(fd)

394 arr.append(line.split()[1:])

395 arr = np.array(arr).astype(float)

396 return arr

397

398 if use_v9_format:

399 energy_header = '--- !EnergyTerms'

400 total_energy_name = 'total_energy_eV'

401

402 def parse_energy(line):

403 return float(line.split(':')[1].strip())

404 else:

405 energy_header = 'Components of total free energy (in Hartree) :'

406 total_energy_name = '>>>>>>>>> Etotal'

407

408 def parse_energy(line):

409 return float(line.rsplit('=', 2)[1]) * Hartree

410

411 for line in fd:

412 if 'cartesian coordinates (angstrom) at end' in line:

413 positions = read_array(fd, natoms)

414 if 'cartesian forces (eV/Angstrom) at end' in line:

415 results['forces'] = read_array(fd, natoms)

416 if 'Cartesian components of stress tensor (hartree/bohr^3)' in line:

417 results['stress'] = read_stress(fd)

418

419 if line.strip() == energy_header:

420 # Header not to be confused with EnergyTermsDC,

421 # therefore we don't use .startswith()

422 energy = None

423 for line in fd:

424 # Which of the listed energies should we include?

425 if total_energy_name in line:

426 energy = parse_energy(line)

427 break

428 if energy is None:

429 raise RuntimeError('No energy found in output')

430 results['energy'] = results['free_energy'] = energy

431

432 if 'END DATASET(S)' in line:

433 break

434

435 znucl_int = znucl.astype(int)

436 znucl_int[znucl_int != znucl] = 0 # (Fractional Z)

437 numbers = znucl_int[types - 1]

438

439 atoms = Atoms(numbers=numbers,

440 positions=positions,

441 cell=cell,

442 pbc=True)

443

444 calc_results = {name: results[name] for name in

445 set(all_properties) & set(results)}

446 atoms.calc = SinglePointCalculator(atoms,

447 **calc_results)

448 atoms.calc.name = "abinit"

449

450 results['atoms'] = atoms

451 return results

452

453

454def match_kpt_header(line):

455 headerpattern = (r'\s*kpt#\s*\S+\s*'

456 r'nband=\s*(\d+),\s*'

457 r'wtk=\s*(\S+?),\s*'

458 r'kpt=\s*(\S+)+\s*(\S+)\s*(\S+)')

459 m = re.match(headerpattern, line)

460 assert m is not None, line

461 nbands = int(m.group(1))

462 weight = float(m.group(2))

463 kvector = np.array(m.group(3, 4, 5)).astype(float)

464 return nbands, weight, kvector

465

466

467def read_eigenvalues_for_one_spin(fd, nkpts):

468

469 kpoint_weights = []

470 kpoint_coords = []

471

472 eig_kn = []

473 for ikpt in range(nkpts):

474 header = next(fd)

475 nbands, weight, kvector = match_kpt_header(header)

476 kpoint_coords.append(kvector)

477 kpoint_weights.append(weight)

478

479 eig_n = []

480 while len(eig_n) < nbands:

481 line = next(fd)

482 tokens = line.split()

483 values = np.array(tokens).astype(float) * Hartree

484 eig_n.extend(values)

485 assert len(eig_n) == nbands

486 eig_kn.append(eig_n)

487 assert nbands == len(eig_kn[0])

488

489 kpoint_weights = np.array(kpoint_weights)

490 kpoint_coords = np.array(kpoint_coords)

491 eig_kn = np.array(eig_kn)

492 return kpoint_coords, kpoint_weights, eig_kn

493

494

495def read_eig(fd):

496 line = next(fd)

497 results = {}

498 m = re.match(r'\s*Fermi \(or HOMO\) energy \(hartree\)\s*=\s*(\S+)', line)

499 if m is not None:

500 results['fermilevel'] = float(m.group(1)) * Hartree

501 line = next(fd)

502

503 nspins = 1

504

505 m = re.match(r'\s*Magnetization \(Bohr magneton\)=\s*(\S+)', line)

506 if m is not None:

507 nspins = 2

508 magmom = float(m.group(1))

509 results['magmom'] = magmom

510 line = next(fd)

511

512 if 'Total spin up' in line:

513 assert nspins == 2

514 line = next(fd)

515

516 m = re.match(r'\s*Eigenvalues \(hartree\) for nkpt\s*='

517 r'\s*(\S+)\s*k\s*points', line)

518 if 'SPIN' in line or 'spin' in line:

519 # If using spinpol with fixed magmoms, we don't get the magmoms

520 # listed before now.

521 nspins = 2

522 assert m is not None

523 nkpts = int(m.group(1))

524

525 eig_skn = []

526

527 kpts, weights, eig_kn = read_eigenvalues_for_one_spin(fd, nkpts)

528 nbands = eig_kn.shape[1]

529

530 eig_skn.append(eig_kn)

531 if nspins == 2:

532 line = next(fd)

533 assert 'SPIN DOWN' in line

534 _, _, eig_kn = read_eigenvalues_for_one_spin(fd, nkpts)

535 assert eig_kn.shape == (nkpts, nbands)

536 eig_skn.append(eig_kn)

537 eig_skn = np.array(eig_skn)

538

539 eigshape = (nspins, nkpts, nbands)

540 assert eig_skn.shape == eigshape, (eig_skn.shape, eigshape)

541

542 results['ibz_kpoints'] = kpts

543 results['kpoint_weights'] = weights

544 results['eigenvalues'] = eig_skn

545 return results

546

547

548def get_default_abinit_pp_paths():

549 return os.environ.get('ABINIT_PP_PATH', '.').split(':')

550

551

552def prepare_abinit_input(directory, atoms, properties, parameters,

553 pp_paths=None,

554 raise_exception=True):

555 directory = Path(directory)

556 species = sorted(set(atoms.numbers))

557 if pp_paths is None:

558 pp_paths = get_default_abinit_pp_paths()

559 ppp = get_ppp_list(atoms, species,

560 raise_exception=raise_exception,

561 xc=parameters['xc'],

562 pps=parameters['pps'],

563 search_paths=pp_paths)

564

565 inputfile = directory / 'abinit.in'

566

567 # XXX inappropriate knowledge about choice of outputfile

568 outputfile = directory / 'abinit.abo'

569

570 # Abinit will write to label.txtA if label.txt already exists,

571 # so we remove it if it's there:

572 if outputfile.exists():

573 outputfile.unlink()

574

575 with open(inputfile, 'w') as fd:

576 write_abinit_in(fd, atoms, param=parameters, species=species,

577 pseudos=ppp)

578

579

580def read_abinit_outputs(directory, label):

581 directory = Path(directory)

582 textfilename = directory / f'{label}.abo'

583 results = {}

584 with open(textfilename) as fd:

585 dct = read_abinit_out(fd)

586 results.update(dct)

587

588 # The eigenvalues section in the main file is shortened to

589 # a limited number of kpoints. We read the complete one from

590 # the EIG file then:

591 with open(directory / f'{label}o_EIG') as fd:

592 dct = read_eig(fd)

593 results.update(dct)

594 return results

595

596

597def read_abinit_gsr(filename):

598 import netCDF4

599 data = netCDF4.Dataset(filename)

600 data.set_auto_mask(False)

601 version = data.abinit_version

602

603 typat = data.variables['atom_species'][:]

604 cell = data.variables['primitive_vectors'][:] * Bohr

605 znucl = data.variables['atomic_numbers'][:]

606 xred = data.variables['reduced_atom_positions'][:]

607

608 znucl_int = znucl.astype(int)

609 znucl_int[znucl_int != znucl] = 0 # (Fractional Z)

610 numbers = znucl_int[typat - 1]

611

612 atoms = Atoms(numbers=numbers,

613 scaled_positions=xred,

614 cell=cell,

615 pbc=True)

616

617 results = {}

618

619 def addresult(name, abinit_name, unit=1):

620 if abinit_name not in data.variables:

621 return

622 values = data.variables[abinit_name][:]

623 # Within the netCDF4 dataset, the float variables return a array(float)

624 # The tolist() is here to ensure that the result is of type float

625 if not values.shape:

626 values = values.tolist()

627 results[name] = values * unit

628

629 addresult('energy', 'etotal', Hartree)

630 addresult('free_energy', 'etotal', Hartree)

631 addresult('forces', 'cartesian_forces', Hartree / Bohr)

632 addresult('stress', 'cartesian_stress_tensor', Hartree / Bohr**3)

633

634 atoms.calc = SinglePointCalculator(atoms, **results)

635 atoms.calc.name = 'abinit'

636 results['atoms'] = atoms

637

638 addresult('fermilevel', 'fermie', Hartree)

639 addresult('ibz_kpoints', 'reduced_coordinates_of_kpoints')

640 addresult('eigenvalues', 'eigenvalues', Hartree)

641 addresult('occupations', 'occupations')

642 addresult('kpoint_weights', 'kpoint_weights')

643 results['version'] = version

644

645 return results

646

647

648def get_ppp_list(atoms, species, raise_exception, xc, pps,

649 search_paths):

650 ppp_list = []

651

652 xcname = 'GGA' if xc != 'LDA' else 'LDA'

653 for Z in species:

654 number = abs(Z)

655 symbol = chemical_symbols[number]

656

657 names = []

658 for s in [symbol, symbol.lower()]:

659 for xcn in [xcname, xcname.lower()]:

660 if pps in ['paw']:

661 hghtemplate = '%s-%s-%s.paw' # E.g. "H-GGA-hard-uspp.paw"

662 names.append(hghtemplate % (s, xcn, '*'))

663 names.append(f'{s}[.-_]*.paw')

664 elif pps in ['pawxml']:

665 hghtemplate = '%s.%s%s.xml' # E.g. "H.GGA_PBE-JTH.xml"

666 names.append(hghtemplate % (s, xcn, '*'))

667 names.append(f'{s}[.-_]*.xml')

668 elif pps in ['hgh.k']:

669 hghtemplate = '%s-q%s.hgh.k' # E.g. "Co-q17.hgh.k"

670 names.append(hghtemplate % (s, '*'))

671 names.append(f'{s}[.-_]*.hgh.k')

672 names.append(f'{s}[.-_]*.hgh')

673 elif pps in ['tm']:

674 hghtemplate = '%d%s%s.pspnc' # E.g. "44ru.pspnc"

675 names.append(hghtemplate % (number, s, '*'))

676 names.append(f'{s}[.-_]*.pspnc')

677 elif pps in ['psp8']:

678 hghtemplate = '%s.psp8' # E.g. "Si.psp8"

679 names.append(hghtemplate % (s))

680 elif pps in ['hgh', 'hgh.sc']:

681 hghtemplate = '%d%s.%s.hgh' # E.g. "42mo.6.hgh"

682 # There might be multiple files with different valence

683 # electron counts, so we must choose between

684 # the ordinary and the semicore versions for some elements.

685 #

686 # Therefore we first use glob to get all relevant files,

687 # then pick the correct one afterwards.

688 names.append(hghtemplate % (number, s, '*'))

689 names.append('%d%s%s.hgh' % (number, s, '*'))

690 names.append(f'{s}[.-_]*.hgh')

691 else: # default extension

692 names.append('%02d-%s.%s.%s' % (number, s, xcn, pps))

693 names.append('%02d[.-_]%s*.%s' % (number, s, pps))

694 names.append('%02d%s*.%s' % (number, s, pps))

695 names.append(f'{s}[.-_]*.{pps}')

696

697 found = False

698 for name in names: # search for file names possibilities

699 for path in search_paths: # in all available directories

700 filenames = glob(join(path, name))

701 if not filenames:

702 continue

703 if pps == 'paw':

704 # warning: see download.sh in

705 # abinit-pseudopotentials*tar.gz for additional

706 # information!

707 #

708 # XXXX This is probably buggy, max(filenames) uses

709 # an lexicographic order so 14 < 8, and it's

710 # untested so if I change it I'm sure things will

711 # just be inconsistent. --askhl

712 filenames[0] = max(filenames) # Semicore or hard

713 elif pps == 'hgh':

714 # Lowest valence electron count

715 filenames[0] = min(filenames)

716 elif pps == 'hgh.k':

717 # Semicore - highest electron count

718 filenames[0] = max(filenames)

719 elif pps == 'tm':

720 # Semicore - highest electron count

721 filenames[0] = max(filenames)

722 elif pps == 'hgh.sc':

723 # Semicore - highest electron count

724 filenames[0] = max(filenames)

725

726 if filenames:

727 found = True

728 ppp_list.append(filenames[0])

729 break

730 if found:

731 break

732

733 if not found:

734 ppp_list.append("Provide {}.{}.{}?".format(symbol, '*', pps))

735 if raise_exception:

736 msg = ('Could not find {} pseudopotential {} for {} in {}'

737 .format(xcname.lower(), pps, symbol, search_paths))

738 raise RuntimeError(msg)

739

740 return ppp_list