Coverage for /builds/ase/ase/ase/io/abinit.py: 85.32%

1# fmt: off

3import os

4import re

5from glob import glob

6from os.path import join

7from pathlib import Path

9import numpy as np

11from ase import Atoms

12from ase.calculators.calculator import all_properties

13from ase.calculators.singlepoint import SinglePointCalculator

14from ase.data import chemical_symbols

15from ase.units import Bohr, Hartree, fs

18def read_abinit_in(fd):

19 """Import ABINIT input file.

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

22 """

24 tokens = []

26 for line in fd:

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

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

30 # note that the file can not be scanned sequentially

32 index = tokens.index("acell")

33 unit = 1.0

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

35 unit = Bohr

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

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

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

40 index = tokens.index("natom")

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

43 index = tokens.index("ntypat")

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

46 index = tokens.index("typat")

47 typat = []

48 while len(typat) < natom:

49 token = tokens[index + 1]

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

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

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

53 else:

54 typat.append(int(token))

55 index += 1

56 assert natom == len(typat)

58 index = tokens.index("znucl")

59 znucl = []

60 for i in range(ntypat):

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

63 index = tokens.index("rprim")

64 rprim = []

65 for i in range(3):

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

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

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

70 # create a list with the atomic numbers

71 numbers = []

72 for i in range(natom):

73 ii = typat[i] - 1

74 numbers.append(znucl[ii])

76 # now the positions of the atoms

77 if "xred" in tokens:

78 index = tokens.index("xred")

79 xred = []

80 for i in range(natom):

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

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

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

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

85 pbc=True)

86 else:

87 if "xcart" in tokens:

88 index = tokens.index("xcart")

89 unit = Bohr

90 elif "xangst" in tokens:

91 unit = 1.0

92 index = tokens.index("xangst")

93 else:

94 raise OSError(

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

97 xangs = []

98 for i in range(natom):

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

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

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

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

103

104 try:

105 ii = tokens.index('nsppol')

106 except ValueError:

107 nsppol = None

108 else:

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

110

111 if nsppol == 2:

112 index = tokens.index('spinat')

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

114 atoms.set_initial_magnetic_moments(magmoms)

115

116 assert len(atoms) == natom

117 return atoms

118

119

120keys_with_units = {

121 'toldfe': 'eV',

122 'tsmear': 'eV',

123 'paoenergyshift': 'eV',

124 'zmunitslength': 'Bohr',

125 'zmunitsangle': 'rad',

126 'zmforcetollength': 'eV/Ang',

127 'zmforcetolangle': 'eV/rad',

128 'zmmaxdispllength': 'Ang',

129 'zmmaxdisplangle': 'rad',

130 'ecut': 'eV',

131 'pawecutdg': 'eV',

132 'dmenergytolerance': 'eV',

133 'electronictemperature': 'eV',

134 'oneta': 'eV',

135 'onetaalpha': 'eV',

136 'onetabeta': 'eV',

137 'onrclwf': 'Ang',

138 'onchemicalpotentialrc': 'Ang',

139 'onchemicalpotentialtemperature': 'eV',

140 'mdmaxcgdispl': 'Ang',

141 'mdmaxforcetol': 'eV/Ang',

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

143 'mdlengthtimestep': 'fs',

144 'mdinitialtemperature': 'eV',

145 'mdtargettemperature': 'eV',

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

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

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

149 'mdtaurelax': 'fs',

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

151 'mdfcdispl': 'Ang',

152 'warningminimumatomicdistance': 'Ang',

153 'rcspatial': 'Ang',

154 'kgridcutoff': 'Ang',

155 'latticeconstant': 'Ang'}

156

157

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

159 from ase.calculators.calculator import kpts2mp

160

161 if param is None:

162 param = {}

163

164 if species is None:

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

166

167 inp = dict(param)

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

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

170 inp.pop(key, None)

171

172 smearing = param.get('smearing')

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

174 assert smearing is None

175

176 if smearing is not None:

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

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

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

180

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

182

183 if 'nbands' in param:

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

185 del inp['nbands']

186

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

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

189 if 'ixc' not in param:

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

191 'PBE': 11,

192 'revPBE': 14,

193 'RPBE': 15,

194 'WC': 23}[xc]

195

196 magmoms = atoms.get_initial_magnetic_moments()

197 if magmoms.any():

198 inp['nsppol'] = 2

199 fd.write('spinat\n')

200 for n, M in enumerate(magmoms):

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

202 else:

203 inp['nsppol'] = 1

204

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

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

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

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

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

210 fd.write('shiftk\n')

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

212

213 valid_lists = (list, np.ndarray)

214 for key in sorted(inp):

215 value = inp[key]

216 unit = keys_with_units.get(key)

217 if unit is not None:

218 if 'fs**2' in unit:

219 value /= fs**2

220 elif 'fs' in unit:

221 value /= fs

222 if isinstance(value, valid_lists):

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

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

225 for dim in value:

226 write_list(fd, dim, unit)

227 else:

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

229 write_list(fd, value, unit)

230 else:

231 if unit is None:

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

233 else:

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

235

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

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

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

239 'a sequence of lines')

240 for line in param['raw']:

241 if isinstance(line, tuple):

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

243 else:

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

245

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

247 fd.write('acell\n')

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

249 fd.write('rprim\n')

250 if atoms.cell.rank != 3:

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

252 .format(atoms.cell))

253 for v in atoms.cell:

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

255

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

257

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

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

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

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

262 fd.write('typat')

263 fd.write('\n')

264

265 types = []

266 for Z in atoms.numbers:

267 for n, Zs in enumerate(species):

268 if Z == Zs:

269 types.append(n + 1)

270 n_entries_int = 20 # integer entries per line

271 for n, type in enumerate(types):

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

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

274 fd.write('\n')

275 fd.write('\n')

276

277 if pseudos is not None:

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

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

280 fd.write(line)

281

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

283 fd.write('xcart\n')

284 for pos in atoms.positions / Bohr:

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

286

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

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

289

290

291def write_list(fd, value, unit):

292 for element in value:

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

294 if unit is not None:

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

296 fd.write("\n")

297

298

299def read_stress(fd):

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

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

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

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

304 stress = np.empty(6)

305 for i in range(3):

306 line = next(fd)

307 m = pat.match(line)

308 if m is None:

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

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

311 .format(line, pat))

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

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

314 unit = Hartree / Bohr**3

315 return stress * unit

316

317

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

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

320

321 Example:

322

323 typat 1 1 1 1 1

324 1 1 1 1

325 """

326 tokens = headerline.split()

327 assert tokens[0].isalpha()

328

329 values = tokens[1:]

330 while len(values) < nvalues:

331 line = next(fd)

332 values.extend(line.split())

333 assert len(values) == nvalues

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

335

336

337def read_abinit_out(fd):

338 results = {}

339

340 def skipto(string):

341 for line in fd:

342 if string in line:

343 return line

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

345

346 line = skipto('Version')

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

348 assert m is not None

349 version = m.group(1)

350 results['version'] = version

351

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

353

354 shape_vars = {}

355

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

357

358 for line in fd:

359 if '===============' in line:

360 break

361

362 tokens = line.split()

363 if not tokens:

364 continue

365

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

367 if tokens[0] == key:

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

369

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

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

372

373 if 'znucl' in line:

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

375

376 if 'rprim' in line:

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

378 cell = cell.reshape(3, 3)

379

380 natoms = shape_vars['natom']

381

382 # Skip ahead to results:

383 for line in fd:

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

385 raise RuntimeError(line)

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

387 break

388 else:

389 raise RuntimeError('Cannot find results section')

390

391 def read_array(fd, nlines):

392 arr = []

393 for i in range(nlines):

394 line = next(fd)

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

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

397 return arr

398

399 if use_v9_format:

400 energy_header = '--- !EnergyTerms'

401 total_energy_name = 'total_energy_eV'

402

403 def parse_energy(line):

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

405 else:

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

407 total_energy_name = '>>>>>>>>> Etotal'

408

409 def parse_energy(line):

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

411

412 for line in fd:

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

414 positions = read_array(fd, natoms)

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

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

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

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

419

420 if line.strip() == energy_header:

421 # Header not to be confused with EnergyTermsDC,

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

423 energy = None

424 for line in fd:

425 # Which of the listed energies should we include?

426 if total_energy_name in line:

427 energy = parse_energy(line)

428 break

429 if energy is None:

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

431 results['energy'] = results['free_energy'] = energy

432

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

434 break

435

436 znucl_int = znucl.astype(int)

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

438 numbers = znucl_int[types - 1]

439

440 atoms = Atoms(numbers=numbers,

441 positions=positions,

442 cell=cell,

443 pbc=True)

444

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

446 set(all_properties) & set(results)}

447 atoms.calc = SinglePointCalculator(atoms,

448 **calc_results)

449 atoms.calc.name = "abinit"

450

451 results['atoms'] = atoms

452 return results

453

454

455def match_kpt_header(line):

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

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

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

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

460 m = re.match(headerpattern, line)

461 assert m is not None, line

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

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

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

465 return nbands, weight, kvector

466

467

468def read_eigenvalues_for_one_spin(fd, nkpts):

469

470 kpoint_weights = []

471 kpoint_coords = []

472

473 eig_kn = []

474 for _ in range(nkpts):

475 header = next(fd)

476 nbands, weight, kvector = match_kpt_header(header)

477 kpoint_coords.append(kvector)

478 kpoint_weights.append(weight)

479

480 eig_n = []

481 while len(eig_n) < nbands:

482 line = next(fd)

483 tokens = line.split()

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

485 eig_n.extend(values)

486 assert len(eig_n) == nbands

487 eig_kn.append(eig_n)

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

489

490 kpoint_weights = np.array(kpoint_weights)

491 kpoint_coords = np.array(kpoint_coords)

492 eig_kn = np.array(eig_kn)

493 return kpoint_coords, kpoint_weights, eig_kn

494

495

496def read_eig(fd):

497 line = next(fd)

498 results = {}

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

500 if m is not None:

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

502 line = next(fd)

503

504 nspins = 1

505

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

507 if m is not None:

508 nspins = 2

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

510 results['magmom'] = magmom

511 line = next(fd)

512

513 if 'Total spin up' in line:

514 assert nspins == 2

515 line = next(fd)

516

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

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

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

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

521 # listed before now.

522 nspins = 2

523 assert m is not None

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

525

526 eig_skn = []

527

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

529 nbands = eig_kn.shape[1]

530

531 eig_skn.append(eig_kn)

532 if nspins == 2:

533 line = next(fd)

534 assert 'SPIN DOWN' in line

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

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

537 eig_skn.append(eig_kn)

538 eig_skn = np.array(eig_skn)

539

540 eigshape = (nspins, nkpts, nbands)

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

542

543 results['ibz_kpoints'] = kpts

544 results['kpoint_weights'] = weights

545 results['eigenvalues'] = eig_skn

546 return results

547

548

549def get_default_abinit_pp_paths():

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

551

552

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

554 pp_paths=None,

555 raise_exception=True):

556 directory = Path(directory)

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

558 if pp_paths is None:

559 pp_paths = get_default_abinit_pp_paths()

560 ppp = get_ppp_list(atoms, species,

561 raise_exception=raise_exception,

562 xc=parameters['xc'],

563 pps=parameters['pps'],

564 search_paths=pp_paths)

565

566 inputfile = directory / 'abinit.in'

567

568 # XXX inappropriate knowledge about choice of outputfile

569 outputfile = directory / 'abinit.abo'

570

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

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

573 if outputfile.exists():

574 outputfile.unlink()

575

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

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

578 pseudos=ppp)

579

580

581def read_abinit_outputs(directory, label):

582 directory = Path(directory)

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

584 results = {}

585 with open(textfilename) as fd:

586 dct = read_abinit_out(fd)

587 results.update(dct)

588

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

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

591 # the EIG file then:

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

593 dct = read_eig(fd)

594 results.update(dct)

595 return results

596

597

598def read_abinit_gsr(filename):

599 import netCDF4

600 data = netCDF4.Dataset(filename)

601 data.set_auto_mask(False)

602 version = data.abinit_version

603

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

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

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

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

608

609 znucl_int = znucl.astype(int)

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

611 numbers = znucl_int[typat - 1]

612

613 atoms = Atoms(numbers=numbers,

614 scaled_positions=xred,

615 cell=cell,

616 pbc=True)

617

618 results = {}

619

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

621 if abinit_name not in data.variables:

622 return

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

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

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

626 if not values.shape:

627 values = values.tolist()

628 results[name] = values * unit

629

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

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

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

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

634

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

636 atoms.calc.name = 'abinit'

637 results['atoms'] = atoms

638

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

640 addresult('ibz_kpoints', 'reduced_coordinates_of_kpoints')

641 addresult('eigenvalues', 'eigenvalues', Hartree)

642 addresult('occupations', 'occupations')

643 addresult('kpoint_weights', 'kpoint_weights')

644 results['version'] = version

645

646 return results

647

648

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

650 search_paths):

651 ppp_list = []

652

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

654 for Z in species:

655 number = abs(Z)

656 symbol = chemical_symbols[number]

657

658 names = []

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

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

661 if pps in ['paw']:

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

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

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

665 elif pps in ['pawxml']:

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

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

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

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

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

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

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

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

674 elif pps in ['tm']:

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

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

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

678 elif pps in ['psp8']:

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

680 names.append(hghtemplate % (s))

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

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

683 # There might be multiple files with different valence

684 # electron counts, so we must choose between

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

686 #

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

688 # then pick the correct one afterwards.

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

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

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

692 else: # default extension

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

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

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

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

697

698 found = False

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

700 for path in search_paths: # in all available directories

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

702 if not filenames:

703 continue

704 if pps == 'paw':

705 # warning: see download.sh in

706 # abinit-pseudopotentials*tar.gz for additional

707 # information!

708 #

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

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

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

712 # just be inconsistent. --askhl

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

714 elif pps == 'hgh':

715 # Lowest valence electron count

716 filenames[0] = min(filenames)

717 elif pps == 'hgh.k':

718 # Semicore - highest electron count

719 filenames[0] = max(filenames)

720 elif pps == 'tm':

721 # Semicore - highest electron count

722 filenames[0] = max(filenames)

723 elif pps == 'hgh.sc':

724 # Semicore - highest electron count

725 filenames[0] = max(filenames)

726

727 if filenames:

728 found = True

729 ppp_list.append(filenames[0])

730 break

731 if found:

732 break

733

734 if not found:

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

736 if raise_exception:

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

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

739 raise RuntimeError(msg)

740

741 return ppp_list