1import os 

2from os.path import join 

3import re 

4from glob import glob 

5from pathlib import Path 

6 

7import numpy as np 

8 

9from ase import Atoms 

10from ase.data import chemical_symbols 

11from ase.units import Hartree, Bohr, fs 

12 

13 

14def read_abinit_in(fd): 

15 """Import ABINIT input file. 

16 

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

18 """ 

19 

20 tokens = [] 

21 

22 for line in fd: 

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

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

25 

26 # note that the file can not be scanned sequentially 

27 

28 index = tokens.index("acell") 

29 unit = 1.0 

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

31 unit = Bohr 

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

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

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

35 

36 index = tokens.index("natom") 

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

38 

39 index = tokens.index("ntypat") 

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

41 

42 index = tokens.index("typat") 

43 typat = [] 

44 while len(typat) < natom: 

45 token = tokens[index + 1] 

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

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

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

49 else: 

50 typat.append(int(token)) 

51 index += 1 

52 assert natom == len(typat) 

53 

54 index = tokens.index("znucl") 

55 znucl = [] 

56 for i in range(ntypat): 

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

58 

59 index = tokens.index("rprim") 

60 rprim = [] 

61 for i in range(3): 

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

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

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

65 

66 # create a list with the atomic numbers 

67 numbers = [] 

68 for i in range(natom): 

69 ii = typat[i] - 1 

70 numbers.append(znucl[ii]) 

71 

72 # now the positions of the atoms 

73 if "xred" in tokens: 

74 index = tokens.index("xred") 

75 xred = [] 

76 for i in range(natom): 

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

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

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

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

81 pbc=True) 

82 else: 

83 if "xcart" in tokens: 

84 index = tokens.index("xcart") 

85 unit = Bohr 

86 elif "xangst" in tokens: 

87 unit = 1.0 

88 index = tokens.index("xangst") 

89 else: 

90 raise IOError( 

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

92 

93 xangs = [] 

94 for i in range(natom): 

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

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

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

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

99 

100 try: 

101 ii = tokens.index('nsppol') 

102 except ValueError: 

103 nsppol = None 

104 else: 

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

106 

107 if nsppol == 2: 

108 index = tokens.index('spinat') 

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

110 atoms.set_initial_magnetic_moments(magmoms) 

111 

112 assert len(atoms) == natom 

113 return atoms 

114 

115 

116keys_with_units = { 

117 'toldfe': 'eV', 

118 'tsmear': 'eV', 

119 'paoenergyshift': 'eV', 

120 'zmunitslength': 'Bohr', 

121 'zmunitsangle': 'rad', 

122 'zmforcetollength': 'eV/Ang', 

123 'zmforcetolangle': 'eV/rad', 

124 'zmmaxdispllength': 'Ang', 

125 'zmmaxdisplangle': 'rad', 

126 'ecut': 'eV', 

127 'pawecutdg': 'eV', 

128 'dmenergytolerance': 'eV', 

129 'electronictemperature': 'eV', 

130 'oneta': 'eV', 

131 'onetaalpha': 'eV', 

132 'onetabeta': 'eV', 

133 'onrclwf': 'Ang', 

134 'onchemicalpotentialrc': 'Ang', 

135 'onchemicalpotentialtemperature': 'eV', 

136 'mdmaxcgdispl': 'Ang', 

137 'mdmaxforcetol': 'eV/Ang', 

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

139 'mdlengthtimestep': 'fs', 

140 'mdinitialtemperature': 'eV', 

141 'mdtargettemperature': 'eV', 

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

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

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

145 'mdtaurelax': 'fs', 

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

147 'mdfcdispl': 'Ang', 

148 'warningminimumatomicdistance': 'Ang', 

149 'rcspatial': 'Ang', 

150 'kgridcutoff': 'Ang', 

151 'latticeconstant': 'Ang'} 

152 

153 

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

155 from ase.calculators.calculator import kpts2mp 

156 

157 if param is None: 

158 param = {} 

159 

160 if species is None: 

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

162 

163 inp = dict(param) 

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

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

166 inp.pop(key, None) 

167 

168 smearing = param.get('smearing') 

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

170 assert smearing is None 

171 

172 if smearing is not None: 

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

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

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

176 

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

178 

179 if 'nbands' in param: 

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

181 del inp['nbands'] 

182 

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

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

185 if 'ixc' not in param: 

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

187 'PBE': 11, 

188 'revPBE': 14, 

189 'RPBE': 15, 

190 'WC': 23}[xc] 

191 

192 magmoms = atoms.get_initial_magnetic_moments() 

193 if magmoms.any(): 

194 inp['nsppol'] = 2 

195 fd.write('spinat\n') 

196 for n, M in enumerate(magmoms): 

197 fd.write('%.14f %.14f %.14f\n' % (0, 0, M)) 

198 else: 

199 inp['nsppol'] = 1 

200 

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

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

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

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

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

206 fd.write('shiftk\n') 

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

208 

209 valid_lists = (list, np.ndarray) 

210 for key in sorted(inp): 

211 value = inp[key] 

212 unit = keys_with_units.get(key) 

213 if unit is not None: 

214 if 'fs**2' in unit: 

215 value /= fs**2 

216 elif 'fs' in unit: 

217 value /= fs 

218 if isinstance(value, valid_lists): 

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

220 fd.write("{}\n".format(key)) 

221 for dim in value: 

222 write_list(fd, dim, unit) 

223 else: 

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

225 write_list(fd, value, unit) 

226 else: 

227 if unit is None: 

228 fd.write("{} {}\n".format(key, value)) 

229 else: 

230 fd.write("{} {} {}\n".format(key, value, unit)) 

231 

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

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

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

235 'a sequence of lines') 

236 for line in param['raw']: 

237 if isinstance(line, tuple): 

238 fd.write(' '.join(['%s' % x for x in line]) + '\n') 

239 else: 

240 fd.write('%s\n' % line) 

241 

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

243 fd.write('acell\n') 

244 fd.write('%.14f %.14f %.14f Angstrom\n' % (1.0, 1.0, 1.0)) 

245 fd.write('rprim\n') 

246 if atoms.cell.rank != 3: 

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

248 .format(atoms.cell)) 

249 for v in atoms.cell: 

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

251 

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

253 

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

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

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

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

258 fd.write('typat') 

259 fd.write('\n') 

260 

261 types = [] 

262 for Z in atoms.numbers: 

263 for n, Zs in enumerate(species): 

264 if Z == Zs: 

265 types.append(n + 1) 

266 n_entries_int = 20 # integer entries per line 

267 for n, type in enumerate(types): 

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

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

270 fd.write('\n') 

271 fd.write('\n') 

272 

273 if pseudos is not None: 

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

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

276 fd.write(line) 

277 

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

279 fd.write('xcart\n') 

280 for pos in atoms.positions / Bohr: 

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

282 

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

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

285 

286 

287def write_list(fd, value, unit): 

288 for element in value: 

289 fd.write("{} ".format(element)) 

290 if unit is not None: 

291 fd.write("{}".format(unit)) 

292 fd.write("\n") 

293 

294 

295def read_stress(fd): 

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

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

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

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

300 stress = np.empty(6) 

301 for i in range(3): 

302 line = next(fd) 

303 m = pat.match(line) 

304 if m is None: 

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

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

307 .format(line, pat)) 

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

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

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

311 unit = Hartree / Bohr**3 

312 return stress * unit 

313 

314 

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

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

317 

318 Example: 

319 

320 typat 1 1 1 1 1 

321 1 1 1 1 

322 """ 

323 tokens = headerline.split() 

324 assert tokens[0].isalpha() 

325 

326 values = tokens[1:] 

327 while len(values) < nvalues: 

328 line = next(fd) 

329 values.extend(line.split()) 

330 assert len(values) == nvalues 

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

332 

333 

334def read_abinit_out(fd): 

335 results = {} 

336 

337 def skipto(string): 

338 for line in fd: 

339 if string in line: 

340 return line 

341 raise RuntimeError('Not found: {}'.format(string)) 

342 

343 line = skipto('Version') 

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

345 assert m is not None 

346 version = m.group(1) 

347 results['version'] = version 

348 

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

350 

351 shape_vars = {} 

352 

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

354 

355 for line in fd: 

356 if '===============' in line: 

357 break 

358 

359 tokens = line.split() 

360 if not tokens: 

361 continue 

362 

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

364 if tokens[0] == key: 

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

366 

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

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

369 

370 if 'znucl' in line: 

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

372 

373 if 'rprim' in line: 

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

375 cell = cell.reshape(3, 3) 

376 

377 natoms = shape_vars['natom'] 

378 

379 # Skip ahead to results: 

380 for line in fd: 

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

382 raise RuntimeError(line) 

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

384 break 

385 else: 

386 raise RuntimeError('Cannot find results section') 

387 

388 def read_array(fd, nlines): 

389 arr = [] 

390 for i in range(nlines): 

391 line = next(fd) 

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

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

394 return arr 

395 

396 if use_v9_format: 

397 energy_header = '--- !EnergyTerms' 

398 total_energy_name = 'total_energy_eV' 

399 

400 def parse_energy(line): 

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

402 else: 

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

404 total_energy_name = '>>>>>>>>> Etotal' 

405 

406 def parse_energy(line): 

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

408 

409 for line in fd: 

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

411 positions = read_array(fd, natoms) 

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

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

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

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

416 

417 if line.strip() == energy_header: 

418 # Header not to be confused with EnergyTermsDC, 

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

420 energy = None 

421 for line in fd: 

422 # Which of the listed energies should we include? 

423 if total_energy_name in line: 

424 energy = parse_energy(line) 

425 break 

426 if energy is None: 

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

428 results['energy'] = results['free_energy'] = energy 

429 

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

431 break 

432 

433 znucl_int = znucl.astype(int) 

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

435 numbers = znucl_int[types - 1] 

436 

437 atoms = Atoms(numbers=numbers, 

438 positions=positions, 

439 cell=cell, 

440 pbc=True) 

441 

442 results['atoms'] = atoms 

443 return results 

444 

445 

446def match_kpt_header(line): 

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

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

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

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

451 m = re.match(headerpattern, line) 

452 assert m is not None, line 

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

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

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

456 return nbands, weight, kvector 

457 

458 

459def read_eigenvalues_for_one_spin(fd, nkpts): 

460 

461 kpoint_weights = [] 

462 kpoint_coords = [] 

463 

464 eig_kn = [] 

465 for ikpt in range(nkpts): 

466 header = next(fd) 

467 nbands, weight, kvector = match_kpt_header(header) 

468 kpoint_coords.append(kvector) 

469 kpoint_weights.append(weight) 

470 

471 eig_n = [] 

472 while len(eig_n) < nbands: 

473 line = next(fd) 

474 tokens = line.split() 

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

476 eig_n.extend(values) 

477 assert len(eig_n) == nbands 

478 eig_kn.append(eig_n) 

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

480 

481 kpoint_weights = np.array(kpoint_weights) 

482 kpoint_coords = np.array(kpoint_coords) 

483 eig_kn = np.array(eig_kn) 

484 return kpoint_coords, kpoint_weights, eig_kn 

485 

486 

487def read_eig(fd): 

488 line = next(fd) 

489 results = {} 

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

491 if m is not None: 

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

493 line = next(fd) 

494 

495 nspins = 1 

496 

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

498 if m is not None: 

499 nspins = 2 

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

501 results['magmom'] = magmom 

502 line = next(fd) 

503 

504 if 'Total spin up' in line: 

505 assert nspins == 2 

506 line = next(fd) 

507 

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

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

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

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

512 # listed before now. 

513 nspins = 2 

514 assert m is not None 

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

516 

517 eig_skn = [] 

518 

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

520 nbands = eig_kn.shape[1] 

521 

522 eig_skn.append(eig_kn) 

523 if nspins == 2: 

524 line = next(fd) 

525 assert 'SPIN DOWN' in line 

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

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

528 eig_skn.append(eig_kn) 

529 eig_skn = np.array(eig_skn) 

530 

531 eigshape = (nspins, nkpts, nbands) 

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

533 

534 results['ibz_kpoints'] = kpts 

535 results['kpoint_weights'] = weights 

536 results['eigenvalues'] = eig_skn 

537 return results 

538 

539 

540def get_default_abinit_pp_paths(): 

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

542 

543 

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

545 pp_paths=None, 

546 raise_exception=True): 

547 directory = Path(directory) 

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

549 if pp_paths is None: 

550 pp_paths = get_default_abinit_pp_paths() 

551 ppp = get_ppp_list(atoms, species, 

552 raise_exception=raise_exception, 

553 xc=parameters['xc'], 

554 pps=parameters['pps'], 

555 search_paths=pp_paths) 

556 

557 inputfile = directory / 'abinit.in' 

558 

559 # XXX inappropriate knowledge about choice of outputfile 

560 outputfile = directory / 'abinit.abo' 

561 

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

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

564 if outputfile.exists(): 

565 outputfile.unlink() 

566 

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

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

569 pseudos=ppp) 

570 

571 

572def read_abinit_outputs(directory, label): 

573 directory = Path(directory) 

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

575 results = {} 

576 with open(textfilename) as fd: 

577 dct = read_abinit_out(fd) 

578 results.update(dct) 

579 

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

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

582 # the EIG file then: 

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

584 dct = read_eig(fd) 

585 results.update(dct) 

586 return results 

587 

588 

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

590 search_paths): 

591 ppp_list = [] 

592 

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

594 for Z in species: 

595 number = abs(Z) 

596 symbol = chemical_symbols[number] 

597 

598 names = [] 

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

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

601 if pps in ['paw']: 

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

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

604 names.append('%s[.-_]*.paw' % s) 

605 elif pps in ['pawxml']: 

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

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

608 names.append('%s[.-_]*.xml' % s) 

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

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

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

612 names.append('%s[.-_]*.hgh.k' % s) 

613 names.append('%s[.-_]*.hgh' % s) 

614 elif pps in ['tm']: 

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

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

617 names.append('%s[.-_]*.pspnc' % s) 

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

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

620 # There might be multiple files with different valence 

621 # electron counts, so we must choose between 

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

623 # 

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

625 # then pick the correct one afterwards. 

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

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

628 names.append('%s[.-_]*.hgh' % s) 

629 else: # default extension 

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

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

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

633 names.append('%s[.-_]*.%s' % (s, pps)) 

634 

635 found = False 

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

637 for path in search_paths: # in all available directories 

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

639 if not filenames: 

640 continue 

641 if pps == 'paw': 

642 # warning: see download.sh in 

643 # abinit-pseudopotentials*tar.gz for additional 

644 # information! 

645 # 

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

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

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

649 # just be inconsistent. --askhl 

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

651 elif pps == 'hgh': 

652 # Lowest valence electron count 

653 filenames[0] = min(filenames) 

654 elif pps == 'hgh.k': 

655 # Semicore - highest electron count 

656 filenames[0] = max(filenames) 

657 elif pps == 'tm': 

658 # Semicore - highest electron count 

659 filenames[0] = max(filenames) 

660 elif pps == 'hgh.sc': 

661 # Semicore - highest electron count 

662 filenames[0] = max(filenames) 

663 

664 if filenames: 

665 found = True 

666 ppp_list.append(filenames[0]) 

667 break 

668 if found: 

669 break 

670 

671 if not found: 

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

673 if raise_exception: 

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

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

676 raise RuntimeError(msg) 

677 

678 return ppp_list