Coverage for /builds/debichem-team/python-ase/ase/io/abinit.py: 85.32%

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 stress[i] = float(m.group(1))

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

312 unit = Hartree / Bohr**3

313 return stress * unit

314

315

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

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

318

319 Example:

320

321 typat 1 1 1 1 1

322 1 1 1 1

323 """

324 tokens = headerline.split()

325 assert tokens[0].isalpha()

326

327 values = tokens[1:]

328 while len(values) < nvalues:

329 line = next(fd)

330 values.extend(line.split())

331 assert len(values) == nvalues

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

333

334

335def read_abinit_out(fd):

336 results = {}

337

338 def skipto(string):

339 for line in fd:

340 if string in line:

341 return line

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

343

344 line = skipto('Version')

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

346 assert m is not None

347 version = m.group(1)

348 results['version'] = version

349

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

351

352 shape_vars = {}

353

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

355

356 for line in fd:

357 if '===============' in line:

358 break

359

360 tokens = line.split()

361 if not tokens:

362 continue

363

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

365 if tokens[0] == key:

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

367

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

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

370

371 if 'znucl' in line:

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

373

374 if 'rprim' in line:

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

376 cell = cell.reshape(3, 3)

377

378 natoms = shape_vars['natom']

379

380 # Skip ahead to results:

381 for line in fd:

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

383 raise RuntimeError(line)

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

385 break

386 else:

387 raise RuntimeError('Cannot find results section')

388

389 def read_array(fd, nlines):

390 arr = []

391 for i in range(nlines):

392 line = next(fd)

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

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

395 return arr

396

397 if use_v9_format:

398 energy_header = '--- !EnergyTerms'

399 total_energy_name = 'total_energy_eV'

400

401 def parse_energy(line):

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

403 else:

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

405 total_energy_name = '>>>>>>>>> Etotal'

406

407 def parse_energy(line):

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

409

410 for line in fd:

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

412 positions = read_array(fd, natoms)

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

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

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

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

417

418 if line.strip() == energy_header:

419 # Header not to be confused with EnergyTermsDC,

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

421 energy = None

422 for line in fd:

423 # Which of the listed energies should we include?

424 if total_energy_name in line:

425 energy = parse_energy(line)

426 break

427 if energy is None:

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

429 results['energy'] = results['free_energy'] = energy

430

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

432 break

433

434 znucl_int = znucl.astype(int)

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

436 numbers = znucl_int[types - 1]

437

438 atoms = Atoms(numbers=numbers,

439 positions=positions,

440 cell=cell,

441 pbc=True)

442

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

444 set(all_properties) & set(results)}

445 atoms.calc = SinglePointCalculator(atoms,

446 **calc_results)

447 atoms.calc.name = "abinit"

448

449 results['atoms'] = atoms

450 return results

451

452

453def match_kpt_header(line):

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

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

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

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

458 m = re.match(headerpattern, line)

459 assert m is not None, line

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

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

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

463 return nbands, weight, kvector

464

465

466def read_eigenvalues_for_one_spin(fd, nkpts):

467

468 kpoint_weights = []

469 kpoint_coords = []

470

471 eig_kn = []

472 for _ in range(nkpts):

473 header = next(fd)

474 nbands, weight, kvector = match_kpt_header(header)

475 kpoint_coords.append(kvector)

476 kpoint_weights.append(weight)

477

478 eig_n = []

479 while len(eig_n) < nbands:

480 line = next(fd)

481 tokens = line.split()

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

483 eig_n.extend(values)

484 assert len(eig_n) == nbands

485 eig_kn.append(eig_n)

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

487

488 kpoint_weights = np.array(kpoint_weights)

489 kpoint_coords = np.array(kpoint_coords)

490 eig_kn = np.array(eig_kn)

491 return kpoint_coords, kpoint_weights, eig_kn

492

493

494def read_eig(fd):

495 line = next(fd)

496 results = {}

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

498 if m is not None:

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

500 line = next(fd)

501

502 nspins = 1

503

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

505 if m is not None:

506 nspins = 2

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

508 results['magmom'] = magmom

509 line = next(fd)

510

511 if 'Total spin up' in line:

512 assert nspins == 2

513 line = next(fd)

514

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

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

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

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

519 # listed before now.

520 nspins = 2

521 assert m is not None

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

523

524 eig_skn = []

525

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

527 nbands = eig_kn.shape[1]

528

529 eig_skn.append(eig_kn)

530 if nspins == 2:

531 line = next(fd)

532 assert 'SPIN DOWN' in line

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

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

535 eig_skn.append(eig_kn)

536 eig_skn = np.array(eig_skn)

537

538 eigshape = (nspins, nkpts, nbands)

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

540

541 results['ibz_kpoints'] = kpts

542 results['kpoint_weights'] = weights

543 results['eigenvalues'] = eig_skn

544 return results

545

546

547def get_default_abinit_pp_paths():

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

549

550

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

552 pp_paths=None,

553 raise_exception=True):

554 directory = Path(directory)

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

556 if pp_paths is None:

557 pp_paths = get_default_abinit_pp_paths()

558 ppp = get_ppp_list(atoms, species,

559 raise_exception=raise_exception,

560 xc=parameters['xc'],

561 pps=parameters['pps'],

562 search_paths=pp_paths)

563

564 inputfile = directory / 'abinit.in'

565

566 # XXX inappropriate knowledge about choice of outputfile

567 outputfile = directory / 'abinit.abo'

568

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

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

571 if outputfile.exists():

572 outputfile.unlink()

573

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

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

576 pseudos=ppp)

577

578

579def read_abinit_outputs(directory, label):

580 directory = Path(directory)

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

582 results = {}

583 with open(textfilename) as fd:

584 dct = read_abinit_out(fd)

585 results.update(dct)

586

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

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

589 # the EIG file then:

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

591 dct = read_eig(fd)

592 results.update(dct)

593 return results

594

595

596def read_abinit_gsr(filename):

597 import netCDF4

598 data = netCDF4.Dataset(filename)

599 data.set_auto_mask(False)

600 version = data.abinit_version

601

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

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

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

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

606

607 znucl_int = znucl.astype(int)

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

609 numbers = znucl_int[typat - 1]

610

611 atoms = Atoms(numbers=numbers,

612 scaled_positions=xred,

613 cell=cell,

614 pbc=True)

615

616 results = {}

617

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

619 if abinit_name not in data.variables:

620 return

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

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

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

624 if not values.shape:

625 values = values.tolist()

626 results[name] = values * unit

627

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

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

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

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

632

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

634 atoms.calc.name = 'abinit'

635 results['atoms'] = atoms

636

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

638 addresult('ibz_kpoints', 'reduced_coordinates_of_kpoints')

639 addresult('eigenvalues', 'eigenvalues', Hartree)

640 addresult('occupations', 'occupations')

641 addresult('kpoint_weights', 'kpoint_weights')

642 results['version'] = version

643

644 return results

645

646

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

648 search_paths):

649 ppp_list = []

650

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

652 for Z in species:

653 number = abs(Z)

654 symbol = chemical_symbols[number]

655

656 names = []

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

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

659 if pps in ['paw']:

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

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

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

663 elif pps in ['pawxml']:

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

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

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

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

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

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

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

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

672 elif pps in ['tm']:

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

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

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

676 elif pps in ['psp8']:

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

678 names.append(hghtemplate % (s))

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

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

681 # There might be multiple files with different valence

682 # electron counts, so we must choose between

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

684 #

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

686 # then pick the correct one afterwards.

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

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

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

690 else: # default extension

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

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

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

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

695

696 found = False

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

698 for path in search_paths: # in all available directories

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

700 if not filenames:

701 continue

702 if pps == 'paw':

703 # warning: see download.sh in

704 # abinit-pseudopotentials*tar.gz for additional

705 # information!

706 #

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

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

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

710 # just be inconsistent. --askhl

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

712 elif pps == 'hgh':

713 # Lowest valence electron count

714 filenames[0] = min(filenames)

715 elif pps == 'hgh.k':

716 # Semicore - highest electron count

717 filenames[0] = max(filenames)

718 elif pps == 'tm':

719 # Semicore - highest electron count

720 filenames[0] = max(filenames)

721 elif pps == 'hgh.sc':

722 # Semicore - highest electron count

723 filenames[0] = max(filenames)

724

725 if filenames:

726 found = True

727 ppp_list.append(filenames[0])

728 break

729 if found:

730 break

731

732 if not found:

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

734 if raise_exception:

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

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

737 raise RuntimeError(msg)

738

739 return ppp_list