1"""This module defines an interface to CASTEP for 

2 use by the ASE (Webpage: http://wiki.fysik.dtu.dk/ase) 

3 

4Authors: 

5 Max Hoffmann, max.hoffmann@ch.tum.de 

6 Joerg Meyer, joerg.meyer@ch.tum.de 

7 Simon P. Rittmeyer, simon.rittmeyer@tum.de 

8 

9Contributors: 

10 Juan M. Lorenzi, juan.lorenzi@tum.de 

11 Georg S. Michelitsch, georg.michelitsch@tch.tum.de 

12 Reinhard J. Maurer, reinhard.maurer@yale.edu 

13 Simone Sturniolo, simone.sturniolo@stfc.ac.uk 

14""" 

15 

16import difflib 

17import numpy as np 

18import os 

19import re 

20import glob 

21import shutil 

22import sys 

23import json 

24import time 

25import tempfile 

26import warnings 

27import subprocess 

28from copy import deepcopy 

29from collections import namedtuple 

30from itertools import product 

31from pathlib import Path 

32from typing import List, Set 

33 

34import ase 

35import ase.units as units 

36from ase.calculators.general import Calculator 

37from ase.calculators.calculator import compare_atoms 

38from ase.calculators.calculator import PropertyNotImplementedError 

39from ase.calculators.calculator import kpts2sizeandoffsets 

40from ase.dft.kpoints import BandPath 

41from ase.parallel import paropen 

42from ase.io.castep import read_param 

43from ase.io.castep import read_bands 

44from ase.constraints import FixCartesian 

45 

46__all__ = [ 

47 'Castep', 

48 'CastepCell', 

49 'CastepParam', 

50 'create_castep_keywords'] 

51 

52contact_email = 'simon.rittmeyer@tum.de' 

53 

54# A convenient table to avoid the previously used "eval" 

55_tf_table = { 

56 '': True, # Just the keyword is equivalent to True 

57 'True': True, 

58 'False': False} 

59 

60 

61def _self_getter(getf): 

62 # A decorator that makes it so that if no 'atoms' argument is passed to a 

63 # getter function, self.atoms is used instead 

64 

65 def decor_getf(self, atoms=None, *args, **kwargs): 

66 

67 if atoms is None: 

68 atoms = self.atoms 

69 

70 return getf(self, atoms, *args, **kwargs) 

71 

72 return decor_getf 

73 

74 

75def _parse_tss_block(value, scaled=False): 

76 # Parse the assigned value for a Transition State Search structure block 

77 is_atoms = isinstance(value, ase.atoms.Atoms) 

78 try: 

79 is_strlist = all(map(lambda x: isinstance(x, str), value)) 

80 except TypeError: 

81 is_strlist = False 

82 

83 if not is_atoms: 

84 if not is_strlist: 

85 # Invalid! 

86 raise TypeError('castep.cell.positions_abs/frac_intermediate/' 

87 'product expects Atoms object or list of strings') 

88 

89 # First line must be Angstroms, or nothing 

90 has_units = len(value[0].strip().split()) == 1 

91 if (not scaled) and has_units and value[0].strip() != 'ang': 

92 raise RuntimeError('Only ang units currently supported in castep.' 

93 'cell.positions_abs_intermediate/product') 

94 return '\n'.join(map(str.strip, value)) 

95 else: 

96 text_block = '' if scaled else 'ang\n' 

97 positions = (value.get_scaled_positions() if scaled else 

98 value.get_positions()) 

99 symbols = value.get_chemical_symbols() 

100 for s, p in zip(symbols, positions): 

101 text_block += ' {0} {1:.3f} {2:.3f} {3:.3f}\n'.format(s, *p) 

102 

103 return text_block 

104 

105 

106class Castep(Calculator): 

107 r""" 

108CASTEP Interface Documentation 

109 

110 

111Introduction 

112============ 

113 

114CASTEP_ [1]_ W_ is a software package which uses density functional theory to 

115provide a good atomic-level description of all manner of materials and 

116molecules. CASTEP can give information about total energies, forces and 

117stresses on an atomic system, as well as calculating optimum geometries, band 

118structures, optical spectra, phonon spectra and much more. It can also perform 

119molecular dynamics simulations. 

120 

121The CASTEP calculator interface class offers intuitive access to all CASTEP 

122settings and most results. All CASTEP specific settings are accessible via 

123attribute access (*i.e*. ``calc.param.keyword = ...`` or 

124``calc.cell.keyword = ...``) 

125 

126 

127Getting Started: 

128================ 

129 

130Set the environment variables appropriately for your system. 

131 

132>>> export CASTEP_COMMAND=' ... ' 

133>>> export CASTEP_PP_PATH=' ... ' 

134 

135Note: alternatively to CASTEP_PP_PATH one can set PSPOT_DIR 

136as CASTEP consults this by default, i.e. 

137 

138>>> export PSPOT_DIR=' ... ' 

139 

140 

141Running the Calculator 

142====================== 

143 

144The default initialization command for the CASTEP calculator is 

145 

146.. class:: Castep(directory='CASTEP', label='castep') 

147 

148To do a minimal run one only needs to set atoms, this will use all 

149default settings of CASTEP, meaning LDA, singlepoint, etc.. 

150 

151With a generated *castep_keywords.json* in place all options are accessible 

152by inspection, *i.e.* tab-completion. This works best when using ``ipython``. 

153All options can be accessed via ``calc.param.<TAB>`` or ``calc.cell.<TAB>`` 

154and documentation is printed with ``calc.param.<keyword> ?`` or 

155``calc.cell.<keyword> ?``. All options can also be set directly 

156using ``calc.keyword = ...`` or ``calc.KEYWORD = ...`` or even 

157``ialc.KeYwOrD`` or directly as named arguments in the call to the constructor 

158(*e.g.* ``Castep(task='GeometryOptimization')``). 

159If using this calculator on a machine without CASTEP, one might choose to copy 

160a *castep_keywords.json* file generated elsewhere in order to access this 

161feature: the file will be used if located in the working directory, 

162*$HOME/.ase/* or *ase/ase/calculators/* within the ASE library. The file should 

163be generated the first time it is needed, but you can generate a new keywords 

164file in the currect directory with ``python -m ase.calculators.castep``. 

165 

166All options that go into the ``.param`` file are held in an ``CastepParam`` 

167instance, while all options that go into the ``.cell`` file and don't belong 

168to the atoms object are held in an ``CastepCell`` instance. Each instance can 

169be created individually and can be added to calculators by attribute 

170assignment, *i.e.* ``calc.param = param`` or ``calc.cell = cell``. 

171 

172All internal variables of the calculator start with an underscore (_). 

173All cell attributes that clearly belong into the atoms object are blocked. 

174Setting ``calc.atoms_attribute`` (*e.g.* ``= positions``) is sent directly to 

175the atoms object. 

176 

177 

178Arguments: 

179========== 

180 

181========================= ==================================================== 

182Keyword Description 

183========================= ==================================================== 

184``directory`` The relative path where all input and output files 

185 will be placed. If this does not exist, it will be 

186 created. Existing directories will be moved to 

187 directory-TIMESTAMP unless self._rename_existing_dir 

188 is set to false. 

189 

190``label`` The prefix of .param, .cell, .castep, etc. files. 

191 

192``castep_command`` Command to run castep. Can also be set via the bash 

193 environment variable ``CASTEP_COMMAND``. If none is 

194 given or found, will default to ``castep`` 

195 

196``check_castep_version`` Boolean whether to check if the installed castep 

197 version matches the version from which the available 

198 options were deduced. Defaults to ``False``. 

199 

200``castep_pp_path`` The path where the pseudopotentials are stored. Can 

201 also be set via the bash environment variables 

202 ``PSPOT_DIR`` (preferred) and ``CASTEP_PP_PATH``. 

203 Will default to the current working directory if 

204 none is given or found. Note that pseudopotentials 

205 may be generated on-the-fly if they are not found. 

206 

207``find_pspots`` Boolean whether to search for pseudopotentials in 

208 ``<castep_pp_path>`` or not. If activated, files in 

209 this directory will be checked for typical names. If 

210 files are not found, they will be generated on the 

211 fly, depending on the ``_build_missing_pspots`` 

212 value. A RuntimeError will be raised in case 

213 multiple files per element are found. Defaults to 

214 ``False``. 

215``keyword_tolerance`` Integer to indicate the level of tolerance to apply 

216 validation of any parameters set in the CastepCell 

217 or CastepParam objects against the ones found in 

218 castep_keywords. Levels are as following: 

219 

220 0 = no tolerance, keywords not found in 

221 castep_keywords will raise an exception 

222 

223 1 = keywords not found will be accepted but produce 

224 a warning (default) 

225 

226 2 = keywords not found will be accepted silently 

227 

228 3 = no attempt is made to look for 

229 castep_keywords.json at all 

230``castep_keywords`` Can be used to pass a CastepKeywords object that is 

231 then used with no attempt to actually load a 

232 castep_keywords.json file. Most useful for debugging 

233 and testing purposes. 

234 

235========================= ==================================================== 

236 

237 

238Additional Settings 

239=================== 

240 

241========================= ==================================================== 

242Internal Setting Description 

243========================= ==================================================== 

244``_castep_command`` (``=castep``): the actual shell command used to 

245 call CASTEP. 

246 

247``_check_checkfile`` (``=True``): this makes write_param() only 

248 write a continue or reuse statement if the 

249 addressed .check or .castep_bin file exists in the 

250 directory. 

251 

252``_copy_pspots`` (``=False``): if set to True the calculator will 

253 actually copy the needed pseudo-potential (\*.usp) 

254 file, usually it will only create symlinks. 

255 

256``_link_pspots`` (``=True``): if set to True the calculator will 

257 actually will create symlinks to the needed pseudo 

258 potentials. Set this option (and ``_copy_pspots``) 

259 to False if you rather want to access your pseudo 

260 potentials using the PSPOT_DIR environment variable 

261 that is read by CASTEP. 

262 *Note:* This option has no effect if ``copy_pspots`` 

263 is True.. 

264 

265``_build_missing_pspots`` (``=True``): if set to True, castep will generate 

266 missing pseudopotentials on the fly. If not, a 

267 RuntimeError will be raised if not all files were 

268 found. 

269 

270``_export_settings`` (``=True``): if this is set to 

271 True, all calculator internal settings shown here 

272 will be included in the .param in a comment line (#) 

273 and can be read again by merge_param. merge_param 

274 can be forced to ignore this directive using the 

275 optional argument ``ignore_internal_keys=True``. 

276 

277``_force_write`` (``=True``): this controls wether the \*cell and 

278 \*param will be overwritten. 

279 

280``_prepare_input_only`` (``=False``): If set to True, the calculator will 

281 create \*cell und \*param file but not 

282 start the calculation itself. 

283 If this is used to prepare jobs locally 

284 and run on a remote cluster it is recommended 

285 to set ``_copy_pspots = True``. 

286 

287``_castep_pp_path`` (``='.'``) : the place where the calculator 

288 will look for pseudo-potential files. 

289 

290``_find_pspots`` (``=False``): if set to True, the calculator will 

291 try to find the respective pseudopotentials from 

292 <_castep_pp_path>. As long as there are no multiple 

293 files per element in this directory, the auto-detect 

294 feature should be very robust. Raises a RuntimeError 

295 if required files are not unique (multiple files per 

296 element). Non existing pseudopotentials will be 

297 generated, though this could be dangerous. 

298 

299``_rename_existing_dir`` (``=True``) : when using a new instance 

300 of the calculator, this will move directories out of 

301 the way that would be overwritten otherwise, 

302 appending a date string. 

303 

304``_set_atoms`` (``=False``) : setting this to True will overwrite 

305 any atoms object previously attached to the 

306 calculator when reading a \.castep file. By de- 

307 fault, the read() function will only create a new 

308 atoms object if none has been attached and other- 

309 wise try to assign forces etc. based on the atom's 

310 positions. ``_set_atoms=True`` could be necessary 

311 if one uses CASTEP's internal geometry optimization 

312 (``calc.param.task='GeometryOptimization'``) 

313 because then the positions get out of sync. 

314 *Warning*: this option is generally not recommended 

315 unless one knows one really needs it. There should 

316 never be any need, if CASTEP is used as a 

317 single-point calculator. 

318 

319``_track_output`` (``=False``) : if set to true, the interface 

320 will append a number to the label on all input 

321 and output files, where n is the number of calls 

322 to this instance. *Warning*: this setting may con- 

323 sume a lot more disk space because of the additio- 

324 nal \*check files. 

325 

326``_try_reuse`` (``=_track_output``) : when setting this, the in- 

327 terface will try to fetch the reuse file from the 

328 previous run even if _track_output is True. By de- 

329 fault it is equal to _track_output, but may be 

330 overridden. 

331 

332 Since this behavior may not always be desirable for 

333 single-point calculations. Regular reuse for *e.g.* 

334 a geometry-optimization can be achieved by setting 

335 ``calc.param.reuse = True``. 

336 

337``_pedantic`` (``=False``) if set to true, the calculator will 

338 inform about settings probably wasting a lot of CPU 

339 time or causing numerical inconsistencies. 

340 

341========================= ==================================================== 

342 

343Special features: 

344================= 

345 

346 

347``.dryrun_ok()`` 

348 Runs ``castep_command seed -dryrun`` in a temporary directory return True if 

349 all variables initialized ok. This is a fast way to catch errors in the 

350 input. Afterwards _kpoints_used is set. 

351 

352``.merge_param()`` 

353 Takes a filename or filehandler of a .param file or CastepParam instance and 

354 merges it into the current calculator instance, overwriting current settings 

355 

356``.keyword.clear()`` 

357 Can be used on any option like ``calc.param.keyword.clear()`` or 

358 ``calc.cell.keyword.clear()`` to return to the CASTEP default. 

359 

360``.initialize()`` 

361 Creates all needed input in the ``_directory``. This can then copied to and 

362 run in a place without ASE or even python. 

363 

364``.set_pspot('<library>')`` 

365 This automatically sets the pseudo-potential for all present species to 

366 ``<Species>_<library>.usp``. Make sure that ``_castep_pp_path`` is set 

367 correctly. Note that there is no check, if the file actually exists. If it 

368 doesn't castep will crash! You may want to use ``find_pspots()`` instead. 

369 

370``.find_pspots(pspot=<library>, suffix=<suffix>)`` 

371 This automatically searches for pseudopotentials of type 

372 ``<Species>_<library>.<suffix>`` or ``<Species>-<library>.<suffix>`` in 

373 ``castep_pp_path` (make sure this is set correctly). Note that ``<Species>`` 

374 will be searched for case insensitive. Regular expressions are accepted, and 

375 arguments ``'*'`` will be regarded as bash-like wildcards. Defaults are any 

376 ``<library>`` and any ``<suffix>`` from ``['usp', 'UPF', 'recpot']``. If you 

377 have well-organized folders with pseudopotentials of one kind, this should 

378 work with the defaults. 

379 

380``print(calc)`` 

381 Prints a short summary of the calculator settings and atoms. 

382 

383``ase.io.castep.read_seed('path-to/seed')`` 

384 Given you have a combination of seed.{param,cell,castep} this will return an 

385 atoms object with the last ionic positions in the .castep file and all other 

386 settings parsed from the .cell and .param file. If no .castep file is found 

387 the positions are taken from the .cell file. The output directory will be 

388 set to the same directory, only the label is preceded by 'copy_of\_' to 

389 avoid overwriting. 

390 

391``.set_kpts(kpoints)`` 

392 This is equivalent to initialising the calculator with 

393 ``calc = Castep(kpts=kpoints)``. ``kpoints`` can be specified in many 

394 convenient forms: simple Monkhorst-Pack grids can be specified e.g. 

395 ``(2, 2, 3)`` or ``'2 2 3'``; lists of specific weighted k-points can be 

396 given in reciprocal lattice coordinates e.g. 

397 ``[[0, 0, 0, 0.25], [0.25, 0.25, 0.25, 0.75]]``; a dictionary syntax is 

398 available for more complex requirements e.g. 

399 ``{'size': (2, 2, 2), 'gamma': True}`` will give a Gamma-centered 2x2x2 M-P 

400 grid, ``{'density': 10, 'gamma': False, 'even': False}`` will give a mesh 

401 with density of at least 10 Ang (based on the unit cell of currently-attached 

402 atoms) with an odd number of points in each direction and avoiding the Gamma 

403 point. 

404 

405``.set_bandpath(bandpath)`` 

406 This is equivalent to initialialising the calculator with 

407 ``calc=Castep(bandpath=bandpath)`` and may be set simultaneously with *kpts*. 

408 It allows an electronic band structure path to be set up using ASE BandPath 

409 objects. This enables a band structure calculation to be set up conveniently 

410 using e.g. calc.set_bandpath(atoms.cell.bandpath().interpolate(npoints=200)) 

411 

412``.band_structure(bandfile=None)`` 

413 Read a band structure from _seedname.bands_ file. This returns an ase 

414 BandStructure object which may be plotted with e.g. 

415 ``calc.band_structure().plot()`` 

416 

417Notes/Issues: 

418============== 

419 

420* Currently *only* the FixAtoms *constraint* is fully supported for reading and 

421 writing. There is some experimental support for the FixCartesian constraint. 

422 

423* There is no support for the CASTEP *unit system*. Units of eV and Angstrom 

424 are used throughout. In particular when converting total energies from 

425 different calculators, one should check that the same CODATA_ version is 

426 used for constants and conversion factors, respectively. 

427 

428.. _CASTEP: http://www.castep.org/ 

429 

430.. _W: https://en.wikipedia.org/wiki/CASTEP 

431 

432.. _CODATA: https://physics.nist.gov/cuu/Constants/index.html 

433 

434.. [1] S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J. Probert, 

435 K. Refson, M. C. Payne Zeitschrift für Kristallographie 220(5-6) 

436 pp.567- 570 (2005) PDF_. 

437 

438.. _PDF: http://www.tcm.phy.cam.ac.uk/castep/papers/ZKristallogr_2005.pdf 

439 

440 

441End CASTEP Interface Documentation 

442 """ 

443 

444 # Class attributes ! 

445 # keys set through atoms object 

446 atoms_keys = [ 

447 'charges', 

448 'ionic_constraints', 

449 'lattice_abs', 

450 'lattice_cart', 

451 'positions_abs', 

452 'positions_abs_final', 

453 'positions_abs_intermediate', 

454 'positions_frac', 

455 'positions_frac_final', 

456 'positions_frac_intermediate'] 

457 

458 atoms_obj_keys = [ 

459 'dipole', 

460 'energy_free', 

461 'energy_zero', 

462 'fermi', 

463 'forces', 

464 'nbands', 

465 'positions', 

466 'stress', 

467 'pressure'] 

468 

469 internal_keys = [ 

470 '_castep_command', 

471 '_check_checkfile', 

472 '_copy_pspots', 

473 '_link_pspots', 

474 '_find_pspots', 

475 '_build_missing_pspots', 

476 '_directory', 

477 '_export_settings', 

478 '_force_write', 

479 '_label', 

480 '_prepare_input_only', 

481 '_castep_pp_path', 

482 '_rename_existing_dir', 

483 '_set_atoms', 

484 '_track_output', 

485 '_try_reuse', 

486 '_pedantic'] 

487 

488 def __init__(self, directory='CASTEP', label='castep', 

489 castep_command=None, check_castep_version=False, 

490 castep_pp_path=None, find_pspots=False, keyword_tolerance=1, 

491 castep_keywords=None, **kwargs): 

492 

493 self.__name__ = 'Castep' 

494 

495 # initialize the ase.calculators.general calculator 

496 Calculator.__init__(self) 

497 

498 from ase.io.castep import write_cell 

499 self._write_cell = write_cell 

500 

501 if castep_keywords is None: 

502 castep_keywords = CastepKeywords(make_param_dict(), 

503 make_cell_dict(), 

504 [], 

505 [], 

506 0) 

507 if keyword_tolerance < 3: 

508 try: 

509 castep_keywords = import_castep_keywords(castep_command) 

510 except CastepVersionError as e: 

511 if keyword_tolerance == 0: 

512 raise e 

513 else: 

514 warnings.warn(str(e)) 

515 

516 self._kw_tol = keyword_tolerance 

517 keyword_tolerance = max(keyword_tolerance, 2) # 3 not accepted below 

518 self.param = CastepParam(castep_keywords, 

519 keyword_tolerance=keyword_tolerance) 

520 self.cell = CastepCell(castep_keywords, 

521 keyword_tolerance=keyword_tolerance) 

522 

523 ################################### 

524 # Calculator state variables # 

525 ################################### 

526 self._calls = 0 

527 self._castep_version = castep_keywords.castep_version 

528 

529 # collects warning from .castep files 

530 self._warnings = [] 

531 # collects content from *.err file 

532 self._error = None 

533 # warnings raised by the ASE interface 

534 self._interface_warnings = [] 

535 

536 # store to check if recalculation is necessary 

537 self._old_atoms = None 

538 self._old_cell = None 

539 self._old_param = None 

540 

541 ################################### 

542 # Internal keys # 

543 # Allow to tweak the behavior # 

544 ################################### 

545 self._opt = {} 

546 self._castep_command = get_castep_command(castep_command) 

547 self._castep_pp_path = get_castep_pp_path(castep_pp_path) 

548 self._check_checkfile = True 

549 self._copy_pspots = False 

550 self._link_pspots = True 

551 self._find_pspots = find_pspots 

552 self._build_missing_pspots = True 

553 self._directory = os.path.abspath(directory) 

554 self._export_settings = True 

555 self._force_write = True 

556 self._label = label 

557 self._prepare_input_only = False 

558 self._rename_existing_dir = True 

559 self._set_atoms = False 

560 self._track_output = False 

561 self._try_reuse = False 

562 

563 # turn off the pedantic user warnings 

564 self._pedantic = False 

565 

566 # will be set on during runtime 

567 self._seed = None 

568 

569 ################################### 

570 # (Physical) result variables # 

571 ################################### 

572 self.atoms = None 

573 # initialize result variables 

574 self._forces = None 

575 self._energy_total = None 

576 self._energy_free = None 

577 self._energy_0K = None 

578 self._energy_total_corr = None 

579 self._eigenvalues = None 

580 self._efermi = None 

581 self._ibz_kpts = None 

582 self._ibz_weights = None 

583 self._band_structure = None 

584 

585 # dispersion corrections 

586 self._dispcorr_energy_total = None 

587 self._dispcorr_energy_free = None 

588 self._dispcorr_energy_0K = None 

589 

590 # spins and hirshfeld volumes 

591 self._spins = None 

592 self._hirsh_volrat = None 

593 

594 # Mulliken charges 

595 self._mulliken_charges = None 

596 # Hirshfeld charges 

597 self._hirshfeld_charges = None 

598 

599 self._number_of_cell_constraints = None 

600 self._output_verbosity = None 

601 self._stress = None 

602 self._pressure = None 

603 self._unit_cell = None 

604 self._kpoints = None 

605 

606 # pointers to other files used at runtime 

607 self._check_file = None 

608 self._castep_bin_file = None 

609 

610 # plane wave cutoff energy (may be derived during PP generation) 

611 self._cut_off_energy = None 

612 

613 # runtime information 

614 self._total_time = None 

615 self._peak_memory = None 

616 

617 # check version of CASTEP options module against current one 

618 if check_castep_version: 

619 local_castep_version = get_castep_version(self._castep_command) 

620 if not hasattr(self, '_castep_version'): 

621 warnings.warn('No castep version found') 

622 return 

623 if not local_castep_version == self._castep_version: 

624 warnings.warn( 

625 'The options module was generated from version %s ' 

626 'while your are currently using CASTEP version %s' % 

627 (self._castep_version, 

628 get_castep_version(self._castep_command))) 

629 self._castep_version = local_castep_version 

630 

631 # processes optional arguments in kw style 

632 for keyword, value in kwargs.items(): 

633 # first fetch special keywords issued by ASE CLI 

634 if keyword == 'kpts': 

635 self.set_kpts(value) 

636 elif keyword == 'bandpath': 

637 self.set_bandpath(value) 

638 elif keyword == 'xc': 

639 self.xc_functional = value 

640 elif keyword == 'ecut': 

641 self.cut_off_energy = value 

642 else: # the general case 

643 self.__setattr__(keyword, value) 

644 

645 def band_structure(self, bandfile=None): 

646 from ase.spectrum.band_structure import BandStructure 

647 

648 if bandfile is None: 

649 bandfile = os.path.join(self._directory, self._seed) + '.bands' 

650 

651 if not os.path.exists(bandfile): 

652 raise ValueError('Cannot find band file "{}".'.format(bandfile)) 

653 

654 kpts, weights, eigenvalues, efermi = read_bands(bandfile) 

655 

656 # Get definitions of high-symmetry points 

657 special_points = self.atoms.cell.bandpath(npoints=0).special_points 

658 bandpath = BandPath(self.atoms.cell, 

659 kpts=kpts, 

660 special_points=special_points) 

661 return BandStructure(bandpath, eigenvalues, reference=efermi) 

662 

663 def set_bandpath(self, bandpath): 

664 """Set a band structure path from ase.dft.kpoints.BandPath object 

665 

666 This will set the bs_kpoint_list block with a set of specific points 

667 determined in ASE. bs_kpoint_spacing will not be used; to modify the 

668 number of points, consider using e.g. bandpath.resample(density=20) to 

669 obtain a new dense path. 

670 

671 Args: 

672 bandpath (:obj:`ase.dft.kpoints.BandPath` or None): 

673 Set to None to remove list of band structure points. Otherwise, 

674 sampling will follow BandPath parameters. 

675 

676 """ 

677 

678 def clear_bs_keywords(): 

679 bs_keywords = product({'bs_kpoint', 'bs_kpoints'}, 

680 {'path', 'path_spacing', 

681 'list', 

682 'mp_grid', 'mp_spacing', 'mp_offset'}) 

683 for bs_tag in bs_keywords: 

684 setattr(self.cell, '_'.join(bs_tag), None) 

685 

686 if bandpath is None: 

687 clear_bs_keywords() 

688 elif isinstance(bandpath, BandPath): 

689 clear_bs_keywords() 

690 self.cell.bs_kpoint_list = [' '.join(map(str, row)) 

691 for row in bandpath.kpts] 

692 else: 

693 raise TypeError('Band structure path must be an ' 

694 'ase.dft.kpoint.BandPath object') 

695 

696 def set_kpts(self, kpts): 

697 """Set k-point mesh/path using a str, tuple or ASE features 

698 

699 Args: 

700 kpts (None, tuple, str, dict): 

701 

702 This method will set the CASTEP parameters kpoints_mp_grid, 

703 kpoints_mp_offset and kpoints_mp_spacing as appropriate. Unused 

704 parameters will be set to None (i.e. not included in input files.) 

705 

706 If kpts=None, all these parameters are set as unused. 

707 

708 The simplest useful case is to give a 3-tuple of integers specifying 

709 a Monkhorst-Pack grid. This may also be formatted as a string separated 

710 by spaces; this is the format used internally before writing to the 

711 input files. 

712 

713 A more powerful set of features is available when using a python 

714 dictionary with the following allowed keys: 

715 

716 - 'size' (3-tuple of int) mesh of mesh dimensions 

717 - 'density' (float) for BZ sampling density in points per recip. Ang 

718 ( kpoint_mp_spacing = 1 / (2pi * density) ). An explicit MP mesh will 

719 be set to allow for rounding/centering. 

720 - 'spacing' (float) for BZ sampling density for maximum space between 

721 sample points in reciprocal space. This is numerically equivalent to 

722 the inbuilt ``calc.cell.kpoint_mp_spacing``, but will be converted to 

723 'density' to allow for rounding/centering. 

724 - 'even' (bool) to round each direction up to the nearest even number; 

725 set False for odd numbers, leave as None for no odd/even rounding. 

726 - 'gamma' (bool) to offset the Monkhorst-Pack grid to include 

727 (0, 0, 0); set False to offset each direction avoiding 0. 

728 """ 

729 

730 def clear_mp_keywords(): 

731 mp_keywords = product({'kpoint', 'kpoints'}, 

732 {'mp_grid', 'mp_offset', 

733 'mp_spacing', 'list'}) 

734 for kp_tag in mp_keywords: 

735 setattr(self.cell, '_'.join(kp_tag), None) 

736 

737 # Case 1: Clear parameters with set_kpts(None) 

738 if kpts is None: 

739 clear_mp_keywords() 

740 pass 

741 

742 # Case 2: list of explicit k-points with weights 

743 # e.g. [[ 0, 0, 0, 0.125], 

744 # [ 0, -0.5, 0, 0.375], 

745 # [-0.5, 0, -0.5, 0.375], 

746 # [-0.5, -0.5, -0.5, 0.125]] 

747 

748 elif (isinstance(kpts, (tuple, list)) 

749 and isinstance(kpts[0], (tuple, list))): 

750 

751 if not all(map((lambda row: len(row) == 4), kpts)): 

752 raise ValueError( 

753 'In explicit kpt list each row should have 4 elements') 

754 

755 clear_mp_keywords() 

756 self.cell.kpoint_list = [' '.join(map(str, row)) for row in kpts] 

757 

758 # Case 3: list of explicit kpts formatted as list of str 

759 # i.e. the internal format of calc.kpoint_list split on \n 

760 # e.g. ['0 0 0 0.125', '0 -0.5 0 0.375', '-0.5 0 -0.5 0.375'] 

761 elif isinstance(kpts, (tuple, list)) and isinstance(kpts[0], str): 

762 

763 if not all(map((lambda row: len(row.split()) == 4), kpts)): 

764 raise ValueError( 

765 'In explicit kpt list each row should have 4 elements') 

766 

767 clear_mp_keywords() 

768 self.cell.kpoint_list = kpts 

769 

770 # Case 4: list or tuple of MP samples e.g. [3, 3, 2] 

771 elif isinstance(kpts, (tuple, list)) and isinstance(kpts[0], int): 

772 if len(kpts) != 3: 

773 raise ValueError('Monkhorst-pack grid should have 3 values') 

774 clear_mp_keywords() 

775 self.cell.kpoint_mp_grid = '%d %d %d' % tuple(kpts) 

776 

777 # Case 5: str representation of Case 3 e.g. '3 3 2' 

778 elif isinstance(kpts, str): 

779 self.set_kpts([int(x) for x in kpts.split()]) 

780 

781 # Case 6: dict of options e.g. {'size': (3, 3, 2), 'gamma': True} 

782 # 'spacing' is allowed but transformed to 'density' to get mesh/offset 

783 elif isinstance(kpts, dict): 

784 kpts = kpts.copy() 

785 

786 if (kpts.get('spacing') is not None 

787 and kpts.get('density') is not None): 

788 raise ValueError( 

789 'Cannot set kpts spacing and density simultaneously.') 

790 else: 

791 if kpts.get('spacing') is not None: 

792 kpts = kpts.copy() 

793 spacing = kpts.pop('spacing') 

794 kpts['density'] = 1 / (np.pi * spacing) 

795 

796 clear_mp_keywords() 

797 size, offsets = kpts2sizeandoffsets(atoms=self.atoms, **kpts) 

798 self.cell.kpoint_mp_grid = '%d %d %d' % tuple(size) 

799 self.cell.kpoint_mp_offset = '%f %f %f' % tuple(offsets) 

800 

801 # Case 7: some other iterator. Try treating as a list: 

802 elif hasattr(kpts, '__iter__'): 

803 self.set_kpts(list(kpts)) 

804 

805 # Otherwise, give up 

806 else: 

807 raise TypeError('Cannot interpret kpts of this type') 

808 

809 def todict(self, skip_default=True): 

810 """Create dict with settings of .param and .cell""" 

811 dct = {} 

812 dct['param'] = self.param.get_attr_dict() 

813 dct['cell'] = self.cell.get_attr_dict() 

814 

815 return dct 

816 

817 def check_state(self, atoms, tol=1e-15): 

818 """Check for system changes since last calculation.""" 

819 return compare_atoms(self._old_atoms, atoms) 

820 

821 def _castep_find_last_record(self, castep_file): 

822 """Checks wether a given castep file has a regular 

823 ending message following the last banner message. If this 

824 is the case, the line number of the last banner is message 

825 is return, otherwise False. 

826 

827 returns (record_start, record_end, end_found, last_record_complete) 

828 """ 

829 if isinstance(castep_file, str): 

830 castep_file = paropen(castep_file, 'r') 

831 file_opened = True 

832 else: 

833 file_opened = False 

834 record_starts = [] 

835 while True: 

836 line = castep_file.readline() 

837 if (('Welcome' in line or 'Materials Studio' in line) 

838 and 'CASTEP' in line): 

839 record_starts = [castep_file.tell()] + record_starts 

840 if not line: 

841 break 

842 

843 if record_starts == []: 

844 warnings.warn('Could not find CASTEP label in result file: %s.' 

845 ' Are you sure this is a .castep file?' % castep_file) 

846 return 

847 

848 # search for regular end of file 

849 end_found = False 

850 # start to search from record beginning from the back 

851 # and see if 

852 record_end = -1 

853 for record_nr, record_start in enumerate(record_starts): 

854 castep_file.seek(record_start) 

855 while True: 

856 line = castep_file.readline() 

857 if not line: 

858 break 

859 if 'warn' in line.lower(): 

860 self._warnings.append(line) 

861 

862 if 'Finalisation time =' in line: 

863 end_found = True 

864 record_end = castep_file.tell() 

865 break 

866 

867 if end_found: 

868 break 

869 

870 if file_opened: 

871 castep_file.close() 

872 

873 if end_found: 

874 # record_nr == 0 corresponds to the last record here 

875 if record_nr == 0: 

876 return (record_start, record_end, True, True) 

877 else: 

878 return (record_start, record_end, True, False) 

879 else: 

880 return (0, record_end, False, False) 

881 

882 def read(self, castep_file=None): 

883 """Read a castep file into the current instance.""" 

884 

885 _close = True 

886 

887 if castep_file is None: 

888 if self._castep_file: 

889 castep_file = self._castep_file 

890 out = paropen(castep_file, 'r') 

891 else: 

892 warnings.warn('No CASTEP file specified') 

893 return 

894 if not os.path.exists(castep_file): 

895 warnings.warn('No CASTEP file found') 

896 

897 elif isinstance(castep_file, str): 

898 out = paropen(castep_file, 'r') 

899 

900 else: 

901 # in this case we assume that we have a fileobj already, but check 

902 # for attributes in order to avoid extended EAFP blocks. 

903 out = castep_file 

904 

905 # look before you leap... 

906 attributes = ['name', 

907 'seek', 

908 'close', 

909 'readline', 

910 'tell'] 

911 

912 for attr in attributes: 

913 if not hasattr(out, attr): 

914 raise TypeError( 

915 '"castep_file" is neither str nor valid fileobj') 

916 

917 castep_file = out.name 

918 _close = False 

919 

920 if self._seed is None: 

921 self._seed = os.path.splitext(os.path.basename(castep_file))[0] 

922 

923 err_file = '%s.0001.err' % self._seed 

924 if os.path.exists(err_file): 

925 err_file = paropen(err_file) 

926 self._error = err_file.read() 

927 err_file.close() 

928 # we return right-away because it might 

929 # just be here from a previous run 

930 # look for last result, if several CASTEP 

931 # run are appended 

932 

933 record_start, record_end, end_found, _\ 

934 = self._castep_find_last_record(out) 

935 if not end_found: 

936 warnings.warn('No regular end found in %s file. %s' % 

937 (castep_file, self._error)) 

938 if _close: 

939 out.close() 

940 return 

941 # we return here, because the file has no a regular end 

942 

943 # now iterate over last CASTEP output in file to extract information 

944 # could be generalized as well to extract trajectory from file 

945 # holding several outputs 

946 n_cell_const = 0 

947 forces = [] 

948 

949 # HOTFIX: 

950 # we have to initialize the _stress variable as a zero array 

951 # otherwise the calculator crashes upon pickling trajectories 

952 # Alternative would be to raise a NotImplementedError() which 

953 # is also kind of not true, since we can extract stresses if 

954 # the user configures CASTEP to print them in the outfile 

955 # stress = [] 

956 stress = np.zeros([3, 3]) 

957 hirsh_volrat = [] 

958 

959 # Two flags to check whether spin-polarized or not, and whether 

960 # Hirshfeld volumes are calculated 

961 spin_polarized = False 

962 calculate_hirshfeld = False 

963 mulliken_analysis = False 

964 hirshfeld_analysis = False 

965 kpoints = None 

966 

967 positions_frac_list = [] 

968 mulliken_charges = [] 

969 spins = [] 

970 

971 out.seek(record_start) 

972 while True: 

973 # TODO: add a switch if we have a geometry optimization: record 

974 # atoms objects for intermediate steps. 

975 try: 

976 line = out.readline() 

977 if not line or out.tell() > record_end: 

978 break 

979 elif 'Hirshfeld Analysis' in line: 

980 hirshfeld_charges = [] 

981 

982 hirshfeld_analysis = True 

983 # skip the separating line 

984 line = out.readline() 

985 # this is the headline 

986 line = out.readline() 

987 

988 if 'Charge' in line: 

989 # skip the next separator line 

990 line = out.readline() 

991 while True: 

992 line = out.readline() 

993 fields = line.split() 

994 if len(fields) == 1: 

995 break 

996 else: 

997 hirshfeld_charges.append(float(fields[-1])) 

998 elif 'stress calculation' in line: 

999 if line.split()[-1].strip() == 'on': 

1000 self.param.calculate_stress = True 

1001 elif 'basis set accuracy' in line: 

1002 self.param.basis_precision = line.split()[-1] 

1003 elif 'plane wave basis set cut-off' in line: 

1004 # NB this is set as a private "result" attribute to avoid 

1005 # conflict with input option basis_precision 

1006 cutoff = float(line.split()[-2]) 

1007 self._cut_off_energy = cutoff 

1008 if self.param.basis_precision.value is None: 

1009 self.param.cut_off_energy = cutoff 

1010 elif 'total energy / atom convergence tol.' in line: 

1011 elec_energy_tol = float(line.split()[-2]) 

1012 self.param.elec_energy_tol = elec_energy_tol 

1013 elif 'convergence tolerance window' in line: 

1014 elec_convergence_win = int(line.split()[-2]) 

1015 self.param.elec_convergence_win = elec_convergence_win 

1016 elif re.match(r'\sfinite basis set correction\s*:', line): 

1017 finite_basis_corr = line.split()[-1] 

1018 fbc_possibilities = {'none': 0, 

1019 'manual': 1, 'automatic': 2} 

1020 fbc = fbc_possibilities[finite_basis_corr] 

1021 self.param.finite_basis_corr = fbc 

1022 elif 'Treating system as non-metallic' in line: 

1023 self.param.fix_occupancy = True 

1024 elif 'max. number of SCF cycles:' in line: 

1025 max_no_scf = float(line.split()[-1]) 

1026 self.param.max_scf_cycles = max_no_scf 

1027 elif 'density-mixing scheme' in line: 

1028 mixing_scheme = line.split()[-1] 

1029 self.param.mixing_scheme = mixing_scheme 

1030 elif 'dump wavefunctions every' in line: 

1031 no_dump_cycles = float(line.split()[-3]) 

1032 self.param.num_dump_cycles = no_dump_cycles 

1033 elif 'optimization strategy' in line: 

1034 lspl = line.split(":") 

1035 if lspl[0].strip() != 'optimization strategy': 

1036 # This can happen in iprint: 3 calculations 

1037 continue 

1038 if 'memory' in line: 

1039 self.param.opt_strategy = 'Memory' 

1040 if 'speed' in line: 

1041 self.param.opt_strategy = 'Speed' 

1042 elif 'calculation limited to maximum' in line: 

1043 calc_limit = float(line.split()[-2]) 

1044 self.param.run_time = calc_limit 

1045 elif 'type of calculation' in line: 

1046 lspl = line.split(":") 

1047 if lspl[0].strip() != 'type of calculation': 

1048 # This can happen in iprint: 3 calculations 

1049 continue 

1050 calc_type = lspl[-1] 

1051 calc_type = re.sub(r'\s+', ' ', calc_type) 

1052 calc_type = calc_type.strip() 

1053 if calc_type != 'single point energy': 

1054 calc_type_possibilities = { 

1055 'geometry optimization': 'GeometryOptimization', 

1056 'band structure': 'BandStructure', 

1057 'molecular dynamics': 'MolecularDynamics', 

1058 'optical properties': 'Optics', 

1059 'phonon calculation': 'Phonon', 

1060 'E-field calculation': 'Efield', 

1061 'Phonon followed by E-field': 'Phonon+Efield', 

1062 'transition state search': 'TransitionStateSearch', 

1063 'Magnetic Resonance': 'MagRes', 

1064 'Core level spectra': 'Elnes', 

1065 'Electronic Spectroscopy': 'ElectronicSpectroscopy' 

1066 } 

1067 ctype = calc_type_possibilities[calc_type] 

1068 self.param.task = ctype 

1069 elif 'using functional' in line: 

1070 used_functional = line.split(":")[-1] 

1071 used_functional = re.sub(r'\s+', ' ', used_functional) 

1072 used_functional = used_functional.strip() 

1073 if used_functional != 'Local Density Approximation': 

1074 used_functional_possibilities = { 

1075 'Perdew Wang (1991)': 'PW91', 

1076 'Perdew Burke Ernzerhof': 'PBE', 

1077 'revised Perdew Burke Ernzerhof': 'RPBE', 

1078 'PBE with Wu-Cohen exchange': 'WC', 

1079 'PBE for solids (2008)': 'PBESOL', 

1080 'Hartree-Fock': 'HF', 

1081 'Hartree-Fock +': 'HF-LDA', 

1082 'Screened Hartree-Fock': 'sX', 

1083 'Screened Hartree-Fock + ': 'sX-LDA', 

1084 'hybrid PBE0': 'PBE0', 

1085 'hybrid B3LYP': 'B3LYP', 

1086 'hybrid HSE03': 'HSE03', 

1087 'hybrid HSE06': 'HSE06' 

1088 } 

1089 used_func = used_functional_possibilities[ 

1090 used_functional] 

1091 self.param.xc_functional = used_func 

1092 elif 'output verbosity' in line: 

1093 iprint = int(line.split()[-1][1]) 

1094 if int(iprint) != 1: 

1095 self.param.iprint = iprint 

1096 elif 'treating system as spin-polarized' in line: 

1097 spin_polarized = True 

1098 self.param.spin_polarized = spin_polarized 

1099 elif 'treating system as non-spin-polarized' in line: 

1100 spin_polarized = False 

1101 elif 'Number of kpoints used' in line: 

1102 kpoints = int(line.split('=')[-1].strip()) 

1103 elif 'Unit Cell' in line: 

1104 lattice_real = [] 

1105 lattice_reci = [] 

1106 while True: 

1107 line = out.readline() 

1108 fields = line.split() 

1109 if len(fields) == 6: 

1110 break 

1111 for i in range(3): 

1112 lattice_real.append([float(f) for f in fields[0:3]]) 

1113 lattice_reci.append([float(f) for f in fields[3:7]]) 

1114 line = out.readline() 

1115 fields = line.split() 

1116 elif 'Cell Contents' in line: 

1117 while True: 

1118 line = out.readline() 

1119 if 'Total number of ions in cell' in line: 

1120 n_atoms = int(line.split()[7]) 

1121 if 'Total number of species in cell' in line: 

1122 int(line.split()[7]) 

1123 fields = line.split() 

1124 if len(fields) == 0: 

1125 break 

1126 elif 'Fractional coordinates of atoms' in line: 

1127 species = [] 

1128 custom_species = None # A CASTEP special thing 

1129 positions_frac = [] 

1130 # positions_cart = [] 

1131 while True: 

1132 line = out.readline() 

1133 fields = line.split() 

1134 if len(fields) == 7: 

1135 break 

1136 for n in range(n_atoms): 

1137 spec_custom = fields[1].split(':', 1) 

1138 elem = spec_custom[0] 

1139 if len(spec_custom) > 1 and custom_species is None: 

1140 # Add it to the custom info! 

1141 custom_species = list(species) 

1142 species.append(elem) 

1143 if custom_species is not None: 

1144 custom_species.append(fields[1]) 

1145 positions_frac.append([float(s) for s in fields[3:6]]) 

1146 line = out.readline() 

1147 fields = line.split() 

1148 positions_frac_list.append(positions_frac) 

1149 elif 'Files used for pseudopotentials' in line: 

1150 while True: 

1151 line = out.readline() 

1152 if 'Pseudopotential generated on-the-fly' in line: 

1153 continue 

1154 fields = line.split() 

1155 if (len(fields) >= 2): 

1156 elem, pp_file = fields 

1157 self.cell.species_pot = (elem, pp_file) 

1158 else: 

1159 break 

1160 elif 'k-Points For BZ Sampling' in line: 

1161 # TODO: generalize for non-Monkhorst Pack case 

1162 # (i.e. kpoint lists) - 

1163 # kpoints_offset cannot be read this way and 

1164 # is hence always set to None 

1165 while True: 

1166 line = out.readline() 

1167 if not line.strip(): 

1168 break 

1169 if 'MP grid size for SCF calculation' in line: 

1170 # kpoints = ' '.join(line.split()[-3:]) 

1171 # self.kpoints_mp_grid = kpoints 

1172 # self.kpoints_mp_offset = '0. 0. 0.' 

1173 # not set here anymore because otherwise 

1174 # two calculator objects go out of sync 

1175 # after each calculation triggering unnecessary 

1176 # recalculation 

1177 break 

1178 elif 'Number of cell constraints' in line: 

1179 n_cell_const = int(line.split()[4]) 

1180 elif 'Final energy' in line: 

1181 self._energy_total = float(line.split()[-2]) 

1182 elif 'Final free energy' in line: 

1183 self._energy_free = float(line.split()[-2]) 

1184 elif 'NB est. 0K energy' in line: 

1185 self._energy_0K = float(line.split()[-2]) 

1186 # check if we had a finite basis set correction 

1187 elif 'Total energy corrected for finite basis set' in line: 

1188 self._energy_total_corr = float(line.split()[-2]) 

1189 

1190 # Add support for dispersion correction 

1191 # filtering due to SEDC is done in get_potential_energy 

1192 elif 'Dispersion corrected final energy' in line: 

1193 self._dispcorr_energy_total = float(line.split()[-2]) 

1194 elif 'Dispersion corrected final free energy' in line: 

1195 self._dispcorr_energy_free = float(line.split()[-2]) 

1196 elif 'dispersion corrected est. 0K energy' in line: 

1197 self._dispcorr_energy_0K = float(line.split()[-2]) 

1198 

1199 # remember to remove constraint labels in force components 

1200 # (lacking a space behind the actual floating point number in 

1201 # the CASTEP output) 

1202 elif '******************** Forces *********************'\ 

1203 in line or\ 

1204 '************** Symmetrised Forces ***************'\ 

1205 in line or\ 

1206 '************** Constrained Symmetrised Forces ****'\ 

1207 '**********'\ 

1208 in line or\ 

1209 '******************** Constrained Forces **********'\ 

1210 '**********'\ 

1211 in line or\ 

1212 '******************* Unconstrained Forces *********'\ 

1213 '**********'\ 

1214 in line: 

1215 fix = [] 

1216 fix_cart = [] 

1217 forces = [] 

1218 while True: 

1219 line = out.readline() 

1220 fields = line.split() 

1221 if len(fields) == 7: 

1222 break 

1223 for n in range(n_atoms): 

1224 consd = np.array([0, 0, 0]) 

1225 fxyz = [0, 0, 0] 

1226 for (i, force_component) in enumerate(fields[-4:-1]): 

1227 if force_component.count("(cons'd)") > 0: 

1228 consd[i] = 1 

1229 fxyz[i] = float(force_component.replace( 

1230 "(cons'd)", '')) 

1231 if consd.all(): 

1232 fix.append(n) 

1233 elif consd.any(): 

1234 fix_cart.append(FixCartesian(n, consd)) 

1235 forces.append(fxyz) 

1236 line = out.readline() 

1237 fields = line.split() 

1238 

1239 # add support for Hirshfeld analysis 

1240 elif 'Hirshfeld / free atomic volume :' in line: 

1241 # if we are here, then params must be able to cope with 

1242 # Hirshfeld flag (if castep_keywords.py matches employed 

1243 # castep version) 

1244 calculate_hirshfeld = True 

1245 hirsh_volrat = [] 

1246 while True: 

1247 line = out.readline() 

1248 fields = line.split() 

1249 if len(fields) == 1: 

1250 break 

1251 for n in range(n_atoms): 

1252 hirsh_atom = float(fields[0]) 

1253 hirsh_volrat.append(hirsh_atom) 

1254 while True: 

1255 line = out.readline() 

1256 if 'Hirshfeld / free atomic volume :' in line or\ 

1257 'Hirshfeld Analysis' in line: 

1258 break 

1259 line = out.readline() 

1260 fields = line.split() 

1261 

1262 elif '***************** Stress Tensor *****************'\ 

1263 in line or\ 

1264 '*********** Symmetrised Stress Tensor ***********'\ 

1265 in line: 

1266 stress = [] 

1267 while True: 

1268 line = out.readline() 

1269 fields = line.split() 

1270 if len(fields) == 6: 

1271 break 

1272 for n in range(3): 

1273 stress.append([float(s) for s in fields[2:5]]) 

1274 line = out.readline() 

1275 fields = line.split() 

1276 line = out.readline() 

1277 if "Pressure:" in line: 

1278 self._pressure = float(line.split()[-2]) * units.GPa 

1279 elif ('BFGS: starting iteration' in line 

1280 or 'BFGS: improving iteration' in line): 

1281 if n_cell_const < 6: 

1282 lattice_real = [] 

1283 lattice_reci = [] 

1284 # backup previous configuration first: 

1285 # for highly symmetric systems (where essentially only the 

1286 # stress is optimized, but the atomic positions) positions 

1287 # are only printed once. 

1288 if species: 

1289 prev_species = deepcopy(species) 

1290 if positions_frac: 

1291 prev_positions_frac = deepcopy(positions_frac) 

1292 species = [] 

1293 positions_frac = [] 

1294 forces = [] 

1295 

1296 # HOTFIX: 

1297 # Same reason for the stress initialization as before 

1298 # stress = [] 

1299 stress = np.zeros([3, 3]) 

1300 

1301 # extract info from the Mulliken analysis 

1302 elif 'Atomic Populations' in line: 

1303 # sometimes this appears twice in a castep file 

1304 

1305 mulliken_analysis = True 

1306 # skip the separating line 

1307 line = out.readline() 

1308 # this is the headline 

1309 line = out.readline() 

1310 

1311 if 'Charge' in line: 

1312 # skip the next separator line 

1313 line = out.readline() 

1314 while True: 

1315 line = out.readline() 

1316 fields = line.split() 

1317 if len(fields) == 1: 

1318 break 

1319 

1320 # the check for len==7 is due to CASTEP 18 

1321 # outformat changes 

1322 if spin_polarized: 

1323 if len(fields) != 7: 

1324 spins.append(float(fields[-1])) 

1325 mulliken_charges.append(float(fields[-2])) 

1326 else: 

1327 mulliken_charges.append(float(fields[-1])) 

1328 

1329 # There is actually no good reason to get out of the loop 

1330 # already at this point... or do I miss something? 

1331 # elif 'BFGS: Final Configuration:' in line: 

1332 # break 

1333 elif 'warn' in line.lower(): 

1334 self._warnings.append(line) 

1335 

1336 # fetch some last info 

1337 elif 'Total time' in line: 

1338 pattern = r'.*=\s*([\d\.]+) s' 

1339 self._total_time = float(re.search(pattern, line).group(1)) 

1340 

1341 elif 'Peak Memory Use' in line: 

1342 pattern = r'.*=\s*([\d]+) kB' 

1343 self._peak_memory = int(re.search(pattern, line).group(1)) 

1344 

1345 except Exception as exception: 

1346 sys.stderr.write(line + '|-> line triggered exception: ' 

1347 + str(exception)) 

1348 raise 

1349 

1350 if _close: 

1351 out.close() 

1352 

1353 # in highly summetric crystals, positions and symmetry are only printed 

1354 # upon init, hence we here restore these original values 

1355 if not positions_frac: 

1356 positions_frac = prev_positions_frac 

1357 if not species: 

1358 species = prev_species 

1359 

1360 if not spin_polarized: 

1361 # set to zero spin if non-spin polarized calculation 

1362 spins = np.zeros(len(positions_frac)) 

1363 elif len(spins) != len(positions_frac): 

1364 warnings.warn('Spins could not be read for the atoms despite' 

1365 ' spin-polarized calculation; spins will be ignored') 

1366 spins = np.zeros(len(positions_frac)) 

1367 

1368 positions_frac_atoms = np.array(positions_frac) 

1369 forces_atoms = np.array(forces) 

1370 spins_atoms = np.array(spins) 

1371 

1372 if mulliken_analysis: 

1373 if len(mulliken_charges) != len(positions_frac): 

1374 warnings.warn( 

1375 'Mulliken charges could not be read for the atoms;' 

1376 ' charges will be ignored') 

1377 mulliken_charges = np.zeros(len(positions_frac)) 

1378 mulliken_charges_atoms = np.array(mulliken_charges) 

1379 else: 

1380 mulliken_charges_atoms = np.zeros(len(positions_frac)) 

1381 

1382 if hirshfeld_analysis: 

1383 hirshfeld_charges_atoms = np.array(hirshfeld_charges) 

1384 else: 

1385 hirshfeld_charges_atoms = None 

1386 

1387 if calculate_hirshfeld: 

1388 hirsh_atoms = np.array(hirsh_volrat) 

1389 else: 

1390 hirsh_atoms = np.zeros_like(spins) 

1391 

1392 if self.atoms and not self._set_atoms: 

1393 # compensate for internal reordering of atoms by CASTEP 

1394 # using the fact that the order is kept within each species 

1395 

1396 # positions_frac_ase = self.atoms.get_scaled_positions(wrap=False) 

1397 atoms_assigned = [False] * len(self.atoms) 

1398 

1399 # positions_frac_castep_init = np.array(positions_frac_list[0]) 

1400 positions_frac_castep = np.array(positions_frac_list[-1]) 

1401 

1402 # species_castep = list(species) 

1403 forces_castep = np.array(forces) 

1404 hirsh_castep = np.array(hirsh_volrat) 

1405 spins_castep = np.array(spins) 

1406 mulliken_charges_castep = np.array(mulliken_charges_atoms) 

1407 

1408 # go through the atoms position list and replace 

1409 # with the corresponding one from the 

1410 # castep file corresponding atomic number 

1411 for iase in range(n_atoms): 

1412 for icastep in range(n_atoms): 

1413 if (species[icastep] == self.atoms[iase].symbol 

1414 and not atoms_assigned[icastep]): 

1415 positions_frac_atoms[iase] = \ 

1416 positions_frac_castep[icastep] 

1417 forces_atoms[iase] = np.array(forces_castep[icastep]) 

1418 if iprint > 1 and calculate_hirshfeld: 

1419 hirsh_atoms[iase] = np.array(hirsh_castep[icastep]) 

1420 if spin_polarized: 

1421 # reordering not necessary in case all spins == 0 

1422 spins_atoms[iase] = np.array(spins_castep[icastep]) 

1423 mulliken_charges_atoms[iase] = np.array( 

1424 mulliken_charges_castep[icastep]) 

1425 atoms_assigned[icastep] = True 

1426 break 

1427 

1428 if not all(atoms_assigned): 

1429 not_assigned = [i for (i, assigned) 

1430 in zip(range(len(atoms_assigned)), 

1431 atoms_assigned) if not assigned] 

1432 warnings.warn( 

1433 '%s atoms not assigned. ' 

1434 ' DEBUGINFO: The following atoms where not assigned: %s' % 

1435 (atoms_assigned.count(False), not_assigned)) 

1436 else: 

1437 self.atoms.set_scaled_positions(positions_frac_atoms) 

1438 

1439 else: 

1440 # If no atoms, object has been previously defined 

1441 # we define it here and set the Castep() instance as calculator. 

1442 # This covers the case that we simply want to open a .castep file. 

1443 

1444 # The next time around we will have an atoms object, since 

1445 # set_calculator also set atoms in the calculator. 

1446 if self.atoms: 

1447 constraints = self.atoms.constraints 

1448 else: 

1449 constraints = [] 

1450 atoms = ase.atoms.Atoms(species, 

1451 cell=lattice_real, 

1452 constraint=constraints, 

1453 pbc=True, 

1454 scaled_positions=positions_frac, 

1455 ) 

1456 if custom_species is not None: 

1457 atoms.new_array('castep_custom_species', 

1458 np.array(custom_species)) 

1459 

1460 if self.param.spin_polarized: 

1461 # only set magnetic moments if this was a spin polarized 

1462 # calculation 

1463 # this one fails as is 

1464 atoms.set_initial_magnetic_moments(magmoms=spins_atoms) 

1465 

1466 if mulliken_analysis: 

1467 atoms.set_initial_charges(charges=mulliken_charges_atoms) 

1468 atoms.calc = self 

1469 

1470 self._kpoints = kpoints 

1471 self._forces = forces_atoms 

1472 # stress in .castep file is given in GPa: 

1473 self._stress = np.array(stress) * units.GPa 

1474 self._hirsh_volrat = hirsh_atoms 

1475 self._spins = spins_atoms 

1476 self._mulliken_charges = mulliken_charges_atoms 

1477 self._hirshfeld_charges = hirshfeld_charges_atoms 

1478 

1479 if self._warnings: 

1480 warnings.warn('WARNING: %s contains warnings' % castep_file) 

1481 for warning in self._warnings: 

1482 warnings.warn(warning) 

1483 # reset 

1484 self._warnings = [] 

1485 

1486 # Read in eigenvalues from bands file 

1487 bands_file = castep_file[:-7] + '.bands' 

1488 if (self.param.task.value is not None 

1489 and self.param.task.value.lower() == 'bandstructure'): 

1490 self._band_structure = self.band_structure(bandfile=bands_file) 

1491 else: 

1492 try: 

1493 (self._ibz_kpts, 

1494 self._ibz_weights, 

1495 self._eigenvalues, 

1496 self._efermi) = read_bands(filename=bands_file) 

1497 except FileNotFoundError: 

1498 warnings.warn('Could not load .bands file, eigenvalues and ' 

1499 'Fermi energy are unknown') 

1500 

1501 def get_hirsh_volrat(self): 

1502 """ 

1503 Return the Hirshfeld volumes. 

1504 """ 

1505 return self._hirsh_volrat 

1506 

1507 def get_spins(self): 

1508 """ 

1509 Return the spins from a plane-wave Mulliken analysis. 

1510 """ 

1511 return self._spins 

1512 

1513 def get_mulliken_charges(self): 

1514 """ 

1515 Return the charges from a plane-wave Mulliken analysis. 

1516 """ 

1517 return self._mulliken_charges 

1518 

1519 def get_hirshfeld_charges(self): 

1520 """ 

1521 Return the charges from a Hirshfeld analysis. 

1522 """ 

1523 return self._hirshfeld_charges 

1524 

1525 def get_total_time(self): 

1526 """ 

1527 Return the total runtime 

1528 """ 

1529 return self._total_time 

1530 

1531 def get_peak_memory(self): 

1532 """ 

1533 Return the peak memory usage 

1534 """ 

1535 return self._peak_memory 

1536 

1537 def set_label(self, label): 

1538 """The label is part of each seed, which in turn is a prefix 

1539 in each CASTEP related file. 

1540 """ 

1541 # we may think about changing this in future to set `self._directory` 

1542 # and `self._label`, as one would expect 

1543 self._label = label 

1544 

1545 def set_pspot(self, pspot, elems=None, 

1546 notelems=None, 

1547 clear=True, 

1548 suffix='usp'): 

1549 """Quickly set all pseudo-potentials: Usually CASTEP psp are named 

1550 like <Elem>_<pspot>.<suffix> so this function function only expects 

1551 the <LibraryName>. It then clears any previous pseudopotential 

1552 settings apply the one with <LibraryName> for each element in the 

1553 atoms object. The optional elems and notelems arguments can be used 

1554 to exclusively assign to some species, or to exclude with notelemens. 

1555 

1556 Parameters :: 

1557 

1558 - elems (None) : set only these elements 

1559 - notelems (None): do not set the elements 

1560 - clear (True): clear previous settings 

1561 - suffix (usp): PP file suffix 

1562 """ 

1563 if self._find_pspots: 

1564 if self._pedantic: 

1565 warnings.warn('Warning: <_find_pspots> = True. ' 

1566 'Do you really want to use `set_pspots()`? ' 

1567 'This does not check whether the PP files exist. ' 

1568 'You may rather want to use `find_pspots()` with ' 

1569 'the same <pspot>.') 

1570 

1571 if clear and not elems and not notelems: 

1572 self.cell.species_pot.clear() 

1573 for elem in set(self.atoms.get_chemical_symbols()): 

1574 if elems is not None and elem not in elems: 

1575 continue 

1576 if notelems is not None and elem in notelems: 

1577 continue 

1578 self.cell.species_pot = (elem, '%s_%s.%s' % (elem, pspot, suffix)) 

1579 

1580 def find_pspots(self, pspot='.+', elems=None, 

1581 notelems=None, clear=True, suffix='(usp|UPF|recpot)'): 

1582 r"""Quickly find and set all pseudo-potentials by searching in 

1583 castep_pp_path: 

1584 

1585 This one is more flexible than set_pspots, and also checks if the files 

1586 are actually available from the castep_pp_path. 

1587 

1588 Essentially, the function parses the filenames in <castep_pp_path> and 

1589 does a regex matching. The respective pattern is: 

1590 

1591 r"^(<elem>|<elem.upper()>|elem.lower()>(_|-)<pspot>\.<suffix>$" 

1592 

1593 In most cases, it will be sufficient to not specify anything, if you 

1594 use standard CASTEP USPPs with only one file per element in the 

1595 <castep_pp_path>. 

1596 

1597 The function raises a `RuntimeError` if there is some ambiguity 

1598 (multiple files per element). 

1599 

1600 Parameters :: 

1601 

1602 - pspots ('.+') : as defined above, will be a wildcard if not 

1603 specified. 

1604 - elems (None) : set only these elements 

1605 - notelems (None): do not set the elements 

1606 - clear (True): clear previous settings 

1607 - suffix (usp|UPF|recpot): PP file suffix 

1608 """ 

1609 if clear and not elems and not notelems: 

1610 self.cell.species_pot.clear() 

1611 

1612 if not os.path.isdir(self._castep_pp_path): 

1613 if self._pedantic: 

1614 warnings.warn( 

1615 'Cannot search directory: {} Folder does not exist' 

1616 .format(self._castep_pp_path)) 

1617 return 

1618 

1619 # translate the bash wildcard syntax to regex 

1620 if pspot == '*': 

1621 pspot = '.*' 

1622 if suffix == '*': 

1623 suffix = '.*' 

1624 if pspot == '*': 

1625 pspot = '.*' 

1626 

1627 # GBRV USPPs have a strnage naming schme 

1628 pattern = r'^({elem}|{elem_upper}|{elem_lower})(_|-){pspot}\.{suffix}$' 

1629 

1630 for elem in set(self.atoms.get_chemical_symbols()): 

1631 if elems is not None and elem not in elems: 

1632 continue 

1633 if notelems is not None and elem in notelems: 

1634 continue 

1635 p = pattern.format(elem=elem, 

1636 elem_upper=elem.upper(), 

1637 elem_lower=elem.lower(), 

1638 pspot=pspot, 

1639 suffix=suffix) 

1640 pps = [] 

1641 for f in os.listdir(self._castep_pp_path): 

1642 if re.match(p, f): 

1643 pps.append(f) 

1644 if not pps: 

1645 if self._pedantic: 

1646 warnings.warn('Pseudopotential for species {} not found!' 

1647 .format(elem)) 

1648 elif not len(pps) == 1: 

1649 raise RuntimeError( 

1650 'Pseudopotential for species ''{} not unique!\n' 

1651 .format(elem) 

1652 + 'Found the following files in {}\n' 

1653 .format(self._castep_pp_path) 

1654 + '\n'.join([' {}'.format(pp) for pp in pps]) + 

1655 '\nConsider a stricter search pattern in `find_pspots()`.') 

1656 else: 

1657 self.cell.species_pot = (elem, pps[0]) 

1658 

1659 @property 

1660 def name(self): 

1661 """Return the name of the calculator (string). """ 

1662 return self.__name__ 

1663 

1664 def get_property(self, name, atoms=None, allow_calculation=True): 

1665 # High-level getter for compliance with the database module... 

1666 # in principle this would not be necessary any longer if we properly 

1667 # based this class on `Calculator` 

1668 if name == 'forces': 

1669 return self.get_forces(atoms) 

1670 elif name == 'energy': 

1671 return self.get_potential_energy(atoms) 

1672 elif name == 'stress': 

1673 return self.get_stress(atoms) 

1674 elif name == 'charges': 

1675 return self.get_charges(atoms) 

1676 else: 

1677 raise PropertyNotImplementedError 

1678 

1679 @_self_getter 

1680 def get_forces(self, atoms): 

1681 """Run CASTEP calculation if needed and return forces.""" 

1682 self.update(atoms) 

1683 return np.array(self._forces) 

1684 

1685 @_self_getter 

1686 def get_total_energy(self, atoms): 

1687 """Run CASTEP calculation if needed and return total energy.""" 

1688 self.update(atoms) 

1689 return self._energy_total 

1690 

1691 @_self_getter 

1692 def get_total_energy_corrected(self, atoms): 

1693 """Run CASTEP calculation if needed and return total energy.""" 

1694 self.update(atoms) 

1695 return self._energy_total_corr 

1696 

1697 @_self_getter 

1698 def get_free_energy(self, atoms): 

1699 """Run CASTEP calculation if needed and return free energy. 

1700 Only defined with smearing.""" 

1701 self.update(atoms) 

1702 return self._energy_free 

1703 

1704 @_self_getter 

1705 def get_0K_energy(self, atoms): 

1706 """Run CASTEP calculation if needed and return 0K energy. 

1707 Only defined with smearing.""" 

1708 self.update(atoms) 

1709 return self._energy_0K 

1710 

1711 @_self_getter 

1712 def get_potential_energy(self, atoms, force_consistent=False): 

1713 # here for compatibility with ase/calculators/general.py 

1714 # but accessing only _name variables 

1715 """Return the total potential energy.""" 

1716 self.update(atoms) 

1717 if force_consistent: 

1718 # Assumption: If no dispersion correction is applied, then the 

1719 # respective value will default to None as initialized. 

1720 if self._dispcorr_energy_free is not None: 

1721 return self._dispcorr_energy_free 

1722 else: 

1723 return self._energy_free 

1724 else: 

1725 if self._energy_0K is not None: 

1726 if self._dispcorr_energy_0K is not None: 

1727 return self._dispcorr_energy_0K 

1728 else: 

1729 return self._energy_0K 

1730 else: 

1731 if self._dispcorr_energy_total is not None: 

1732 return self._dispcorr_energy_total 

1733 else: 

1734 if self._energy_total_corr is not None: 

1735 return self._energy_total_corr 

1736 else: 

1737 return self._energy_total 

1738 

1739 @_self_getter 

1740 def get_stress(self, atoms): 

1741 """Return the stress.""" 

1742 self.update(atoms) 

1743 # modification: we return the Voigt form directly to get rid of the 

1744 # annoying user warnings 

1745 stress = np.array( 

1746 [self._stress[0, 0], self._stress[1, 1], self._stress[2, 2], 

1747 self._stress[1, 2], self._stress[0, 2], self._stress[0, 1]]) 

1748 # return self._stress 

1749 return stress 

1750 

1751 @_self_getter 

1752 def get_pressure(self, atoms): 

1753 """Return the pressure.""" 

1754 self.update(atoms) 

1755 return self._pressure 

1756 

1757 @_self_getter 

1758 def get_unit_cell(self, atoms): 

1759 """Return the unit cell.""" 

1760 self.update(atoms) 

1761 return self._unit_cell 

1762 

1763 @_self_getter 

1764 def get_kpoints(self, atoms): 

1765 """Return the kpoints.""" 

1766 self.update(atoms) 

1767 return self._kpoints 

1768 

1769 @_self_getter 

1770 def get_number_cell_constraints(self, atoms): 

1771 """Return the number of cell constraints.""" 

1772 self.update(atoms) 

1773 return self._number_of_cell_constraints 

1774 

1775 @_self_getter 

1776 def get_charges(self, atoms): 

1777 """Run CASTEP calculation if needed and return Mulliken charges.""" 

1778 self.update(atoms) 

1779 return np.array(self._mulliken_charges) 

1780 

1781 @_self_getter 

1782 def get_magnetic_moments(self, atoms): 

1783 """Run CASTEP calculation if needed and return Mulliken charges.""" 

1784 self.update(atoms) 

1785 return np.array(self._spins) 

1786 

1787 def set_atoms(self, atoms): 

1788 """Sets the atoms for the calculator and vice versa.""" 

1789 atoms.pbc = [True, True, True] 

1790 self.__dict__['atoms'] = atoms.copy() 

1791 self.atoms._calc = self 

1792 

1793 def update(self, atoms): 

1794 """Checks if atoms object or calculator changed and 

1795 runs calculation if so. 

1796 """ 

1797 if self.calculation_required(atoms): 

1798 self.calculate(atoms) 

1799 

1800 def calculation_required(self, atoms, _=None): 

1801 """Checks wether anything changed in the atoms object or CASTEP 

1802 settings since the last calculation using this instance. 

1803 """ 

1804 # SPR: what happens with the atoms parameter here? Why don't we use it? 

1805 # from all that I can tell we need to compare against atoms instead of 

1806 # self.atoms 

1807 # if not self.atoms == self._old_atoms: 

1808 if not atoms == self._old_atoms: 

1809 return True 

1810 if self._old_param is None or self._old_cell is None: 

1811 return True 

1812 if not self.param._options == self._old_param._options: 

1813 return True 

1814 if not self.cell._options == self._old_cell._options: 

1815 return True 

1816 return False 

1817 

1818 def calculate(self, atoms): 

1819 """Write all necessary input file and call CASTEP.""" 

1820 self.prepare_input_files(atoms, force_write=self._force_write) 

1821 if not self._prepare_input_only: 

1822 self.run() 

1823 self.read() 

1824 

1825 # we need to push the old state here! 

1826 # although run() pushes it, read() may change the atoms object 

1827 # again. 

1828 # yet, the old state is supposed to be the one AFTER read() 

1829 self.push_oldstate() 

1830 

1831 def push_oldstate(self): 

1832 """This function pushes the current state of the (CASTEP) Atoms object 

1833 onto the previous state. Or in other words after calling this function, 

1834 calculation_required will return False and enquiry functions just 

1835 report the current value, e.g. get_forces(), get_potential_energy(). 

1836 """ 

1837 # make a snapshot of all current input 

1838 # to be able to test if recalculation 

1839 # is necessary 

1840 self._old_atoms = self.atoms.copy() 

1841 self._old_param = deepcopy(self.param) 

1842 self._old_cell = deepcopy(self.cell) 

1843 

1844 def initialize(self, *args, **kwargs): 

1845 """Just an alias for prepar_input_files to comply with standard 

1846 function names in ASE. 

1847 """ 

1848 self.prepare_input_files(*args, **kwargs) 

1849 

1850 def prepare_input_files(self, atoms=None, force_write=None): 

1851 """Only writes the input .cell and .param files and return 

1852 This can be useful if one quickly needs to prepare input files 

1853 for a cluster where no python or ASE is available. One can than 

1854 upload the file manually and read out the results using 

1855 Castep().read(). 

1856 """ 

1857 

1858 if self.param.reuse.value is None: 

1859 if self._pedantic: 

1860 warnings.warn( 

1861 'You have not set e.g. calc.param.reuse = True. ' 

1862 'Reusing a previous calculation may save CPU time! ' 

1863 'The interface will make sure by default, .check exists. ' 

1864 'file before adding this statement to the .param file.') 

1865 if self.param.num_dump_cycles.value is None: 

1866 if self._pedantic: 

1867 warnings.warn( 

1868 'You have not set e.g. calc.param.num_dump_cycles = 0. ' 

1869 'This can save you a lot of disk space. One only needs ' 

1870 '*wvfn* if electronic convergence is not achieved.') 

1871 from ase.io.castep import write_param 

1872 

1873 if atoms is None: 

1874 atoms = self.atoms 

1875 else: 

1876 self.atoms = atoms 

1877 

1878 if force_write is None: 

1879 force_write = self._force_write 

1880 

1881 # if we have new instance of the calculator, 

1882 # move existing results out of the way, first 

1883 if (os.path.isdir(self._directory) 

1884 and self._calls == 0 

1885 and self._rename_existing_dir): 

1886 if os.listdir(self._directory) == []: 

1887 os.rmdir(self._directory) 

1888 else: 

1889 # rename appending creation date of the directory 

1890 ctime = time.localtime(os.lstat(self._directory).st_ctime) 

1891 os.rename(self._directory, '%s.bak-%s' % 

1892 (self._directory, 

1893 time.strftime('%Y%m%d-%H%M%S', ctime))) 

1894 

1895 # create work directory 

1896 if not os.path.isdir(self._directory): 

1897 os.makedirs(self._directory, 0o775) 

1898 

1899 # we do this every time, not only upon first call 

1900 # if self._calls == 0: 

1901 self._fetch_pspots() 

1902 

1903 # if _try_reuse is requested and this 

1904 # is not the first run, we try to find 

1905 # the .check file from the previous run 

1906 # this is only necessary if _track_output 

1907 # is set to true 

1908 if self._try_reuse and self._calls > 0: 

1909 if os.path.exists(self._abs_path(self._check_file)): 

1910 self.param.reuse = self._check_file 

1911 elif os.path.exists(self._abs_path(self._castep_bin_file)): 

1912 self.param.reuse = self._castep_bin_file 

1913 self._seed = self._build_castep_seed() 

1914 self._check_file = '%s.check' % self._seed 

1915 self._castep_bin_file = '%s.castep_bin' % self._seed 

1916 self._castep_file = self._abs_path('%s.castep' % self._seed) 

1917 

1918 # write out the input file 

1919 self._write_cell(self._abs_path('%s.cell' % self._seed), 

1920 self.atoms, castep_cell=self.cell, 

1921 force_write=force_write) 

1922 

1923 if self._export_settings: 

1924 interface_options = self._opt 

1925 else: 

1926 interface_options = None 

1927 write_param(self._abs_path('%s.param' % self._seed), self.param, 

1928 check_checkfile=self._check_checkfile, 

1929 force_write=force_write, 

1930 interface_options=interface_options,) 

1931 

1932 def _build_castep_seed(self): 

1933 """Abstracts to construction of the final castep <seed> 

1934 with and without _tracking_output. 

1935 """ 

1936 if self._track_output: 

1937 return '%s-%06d' % (self._label, self._calls) 

1938 else: 

1939 return '%s' % (self._label) 

1940 

1941 def _abs_path(self, path): 

1942 # Create an absolute path for a file to put in the working directory 

1943 return os.path.join(self._directory, path) 

1944 

1945 def run(self): 

1946 """Simply call castep. If the first .err file 

1947 contains text, this will be printed to the screen. 

1948 """ 

1949 # change to target directory 

1950 self._calls += 1 

1951 

1952 # run castep itself 

1953 stdout, stderr = shell_stdouterr('%s %s' % (self._castep_command, 

1954 self._seed), 

1955 cwd=self._directory) 

1956 if stdout: 

1957 print('castep call stdout:\n%s' % stdout) 

1958 if stderr: 

1959 print('castep call stderr:\n%s' % stderr) 

1960 

1961 # shouldn't it be called after read()??? 

1962 # self.push_oldstate() 

1963 

1964 # check for non-empty error files 

1965 err_file = self._abs_path('%s.0001.err' % self._seed) 

1966 if os.path.exists(err_file): 

1967 err_file = open(err_file) 

1968 self._error = err_file.read() 

1969 err_file.close() 

1970 if self._error: 

1971 raise RuntimeError(self._error) 

1972 

1973 def __repr__(self): 

1974 """Returns generic, fast to capture representation of 

1975 CASTEP settings along with atoms object. 

1976 """ 

1977 expr = '' 

1978 expr += '-----------------Atoms--------------------\n' 

1979 if self.atoms is not None: 

1980 expr += str('%20s\n' % self.atoms) 

1981 else: 

1982 expr += 'None\n' 

1983 

1984 expr += '-----------------Param keywords-----------\n' 

1985 expr += str(self.param) 

1986 expr += '-----------------Cell keywords------------\n' 

1987 expr += str(self.cell) 

1988 expr += '-----------------Internal keys------------\n' 

1989 for key in self.internal_keys: 

1990 expr += '%20s : %s\n' % (key, self._opt[key]) 

1991 

1992 return expr 

1993 

1994 def __getattr__(self, attr): 

1995 """___getattr___ gets overloaded to reroute the internal keys 

1996 and to be able to easily store them in in the param so that 

1997 they can be read in again in subsequent calls. 

1998 """ 

1999 if attr in self.internal_keys: 

2000 return self._opt[attr] 

2001 if attr in ['__repr__', '__str__']: 

2002 raise AttributeError 

2003 elif attr not in self.__dict__: 

2004 raise AttributeError('Attribute {0} not found'.format(attr)) 

2005 else: 

2006 return self.__dict__[attr] 

2007 

2008 def __setattr__(self, attr, value): 

2009 """We overload the settattr method to make value assignment 

2010 as pythonic as possible. Internal values all start with _. 

2011 Value assigment is case insensitive! 

2012 """ 

2013 

2014 if attr.startswith('_'): 

2015 # internal variables all start with _ 

2016 # let's check first if they are close but not identical 

2017 # to one of the switches, that the user accesses directly 

2018 similars = difflib.get_close_matches(attr, self.internal_keys, 

2019 cutoff=0.9) 

2020 if attr not in self.internal_keys and similars: 

2021 warnings.warn( 

2022 'Warning: You probably tried one of: %s but typed %s' % 

2023 (similars, attr)) 

2024 if attr in self.internal_keys: 

2025 self._opt[attr] = value 

2026 if attr == '_track_output': 

2027 if value: 

2028 self._try_reuse = True 

2029 if self._pedantic: 

2030 warnings.warn( 

2031 'You switched _track_output on. This will ' 

2032 'consume a lot of disk-space. The interface ' 

2033 'also switched _try_reuse on, which will ' 

2034 'try to find the last check file. Set ' 

2035 '_try_reuse = False, if you need ' 

2036 'really separate calculations') 

2037 elif '_try_reuse' in self._opt and self._try_reuse: 

2038 self._try_reuse = False 

2039 if self._pedantic: 

2040 warnings.warn('_try_reuse is set to False, too') 

2041 else: 

2042 self.__dict__[attr] = value 

2043 return 

2044 elif attr in ['atoms', 'cell', 'param']: 

2045 if value is not None: 

2046 if attr == 'atoms' and not isinstance(value, ase.atoms.Atoms): 

2047 raise TypeError( 

2048 '%s is not an instance of ase.atoms.Atoms.' % value) 

2049 elif attr == 'cell' and not isinstance(value, CastepCell): 

2050 raise TypeError('%s is not an instance of CastepCell.' % 

2051 value) 

2052 elif attr == 'param' and not isinstance(value, CastepParam): 

2053 raise TypeError('%s is not an instance of CastepParam.' % 

2054 value) 

2055 # These 3 are accepted right-away, no matter what 

2056 self.__dict__[attr] = value 

2057 return 

2058 elif attr in self.atoms_obj_keys: 

2059 # keywords which clearly belong to the atoms object are 

2060 # rerouted to go there 

2061 self.atoms.__dict__[attr] = value 

2062 return 

2063 elif attr in self.atoms_keys: 

2064 # CASTEP keywords that should go into the atoms object 

2065 # itself are blocked 

2066 warnings.warn('Ignoring setings of "%s", since this has to be set ' 

2067 'through the atoms object' % attr) 

2068 return 

2069 

2070 attr = attr.lower() 

2071 if attr not in (list(self.cell._options.keys()) 

2072 + list(self.param._options.keys())): 

2073 # what is left now should be meant to be a castep keyword 

2074 # so we first check if it defined, and if not offer some error 

2075 # correction 

2076 if self._kw_tol == 0: 

2077 similars = difflib.get_close_matches( 

2078 attr, 

2079 self.cell._options.keys() + self.param._options.keys()) 

2080 if similars: 

2081 raise UserWarning('Option "%s" not known! You mean "%s"?' % 

2082 (attr, similars[0])) 

2083 else: 

2084 raise UserWarning('Option "%s" is not known!' % attr) 

2085 else: 

2086 warnings.warn('Option "%s" is not known - please set any new' 

2087 ' options directly in the .cell or .param ' 

2088 'objects' % attr) 

2089 return 

2090 

2091 # here we know it must go into one of the component param or cell 

2092 # so we first determine which one 

2093 if attr in self.param._options.keys(): 

2094 comp = 'param' 

2095 elif attr in self.cell._options.keys(): 

2096 comp = 'cell' 

2097 else: 

2098 raise UserWarning('Programming error: could not attach ' 

2099 'the keyword to an input file') 

2100 

2101 self.__dict__[comp].__setattr__(attr, value) 

2102 

2103 def merge_param(self, param, overwrite=True, ignore_internal_keys=False): 

2104 """Parse a param file and merge it into the current parameters.""" 

2105 if isinstance(param, CastepParam): 

2106 for key, option in param._options.items(): 

2107 if option.value is not None: 

2108 self.param.__setattr__(key, option.value) 

2109 return 

2110 

2111 elif isinstance(param, str): 

2112 param_file = open(param, 'r') 

2113 _close = True 

2114 

2115 else: 

2116 # in this case we assume that we have a fileobj already, but check 

2117 # for attributes in order to avoid extended EAFP blocks. 

2118 param_file = param 

2119 

2120 # look before you leap... 

2121 attributes = ['name', 

2122 'close' 

2123 'readlines'] 

2124 

2125 for attr in attributes: 

2126 if not hasattr(param_file, attr): 

2127 raise TypeError('"param" is neither CastepParam nor str ' 

2128 'nor valid fileobj') 

2129 

2130 param = param_file.name 

2131 _close = False 

2132 

2133 self, int_opts = read_param(fd=param_file, calc=self, 

2134 get_interface_options=True) 

2135 

2136 # Add the interface options 

2137 for k, val in int_opts.items(): 

2138 if (k in self.internal_keys and not ignore_internal_keys): 

2139 if val in _tf_table: 

2140 val = _tf_table[val] 

2141 self._opt[k] = val 

2142 

2143 if _close: 

2144 param_file.close() 

2145 

2146 def dryrun_ok(self, dryrun_flag='-dryrun'): 

2147 """Starts a CASTEP run with the -dryrun flag [default] 

2148 in a temporary and check wether all variables are initialized 

2149 correctly. This is recommended for every bigger simulation. 

2150 """ 

2151 from ase.io.castep import write_param 

2152 

2153 temp_dir = tempfile.mkdtemp() 

2154 self._fetch_pspots(temp_dir) 

2155 seed = 'dryrun' 

2156 

2157 self._write_cell(os.path.join(temp_dir, '%s.cell' % seed), 

2158 self.atoms, castep_cell=self.cell) 

2159 # This part needs to be modified now that we rely on the new formats.py 

2160 # interface 

2161 if not os.path.isfile(os.path.join(temp_dir, '%s.cell' % seed)): 

2162 warnings.warn('%s.cell not written - aborting dryrun' % seed) 

2163 return 

2164 write_param(os.path.join(temp_dir, '%s.param' % seed), self.param, ) 

2165 

2166 stdout, stderr = shell_stdouterr(('%s %s %s' % (self._castep_command, 

2167 seed, 

2168 dryrun_flag)), 

2169 cwd=temp_dir) 

2170 

2171 if stdout: 

2172 print(stdout) 

2173 if stderr: 

2174 print(stderr) 

2175 result_file = open(os.path.join(temp_dir, '%s.castep' % seed)) 

2176 

2177 txt = result_file.read() 

2178 ok_string = r'.*DRYRUN finished.*No problems found with input files.*' 

2179 match = re.match(ok_string, txt, re.DOTALL) 

2180 

2181 m = re.search(r'Number of kpoints used =\s*([0-9]+)', txt) 

2182 if m: 

2183 self._kpoints = int(m.group(1)) 

2184 else: 

2185 warnings.warn( 

2186 'Couldn\'t fetch number of kpoints from dryrun CASTEP file') 

2187 

2188 err_file = os.path.join(temp_dir, '%s.0001.err' % seed) 

2189 if match is None and os.path.exists(err_file): 

2190 err_file = open(err_file) 

2191 self._error = err_file.read() 

2192 err_file.close() 

2193 

2194 result_file.close() 

2195 shutil.rmtree(temp_dir) 

2196 

2197 # re.match return None is the string does not match 

2198 return match is not None 

2199 

2200 def _fetch_pspots(self, directory=None): 

2201 """Put all specified pseudo-potentials into the working directory. 

2202 """ 

2203 # should be a '==' right? Otherwise setting _castep_pp_path is not 

2204 # honored. 

2205 if (not os.environ.get('PSPOT_DIR', None) 

2206 and self._castep_pp_path == os.path.abspath('.')): 

2207 # By default CASTEP consults the environment variable 

2208 # PSPOT_DIR. If this contains a list of colon separated 

2209 # directories it will check those directories for pseudo- 

2210 # potential files if not in the current directory. 

2211 # Thus if PSPOT_DIR is set there is nothing left to do. 

2212 # If however PSPOT_DIR was been accidentally set 

2213 # (e.g. with regards to a different program) 

2214 # setting CASTEP_PP_PATH to an explicit value will 

2215 # still be honored. 

2216 return 

2217 

2218 if directory is None: 

2219 directory = self._directory 

2220 if not os.path.isdir(self._castep_pp_path): 

2221 warnings.warn('PSPs directory %s not found' % self._castep_pp_path) 

2222 pspots = {} 

2223 if self._find_pspots: 

2224 self.find_pspots() 

2225 if self.cell.species_pot.value is not None: 

2226 for line in self.cell.species_pot.value.split('\n'): 

2227 line = line.split() 

2228 if line: 

2229 pspots[line[0]] = line[1] 

2230 for species in self.atoms.get_chemical_symbols(): 

2231 if not pspots or species not in pspots.keys(): 

2232 if self._build_missing_pspots: 

2233 if self._pedantic: 

2234 warnings.warn( 

2235 'Warning: you have no PP specified for %s. ' 

2236 'CASTEP will now generate an on-the-fly ' 

2237 'potentials. ' 

2238 'For sake of numerical consistency and efficiency ' 

2239 'this is discouraged.' % species) 

2240 else: 

2241 raise RuntimeError( 

2242 'Warning: you have no PP specified for %s.' % 

2243 species) 

2244 if self.cell.species_pot.value: 

2245 for (species, pspot) in pspots.items(): 

2246 orig_pspot_file = os.path.join(self._castep_pp_path, pspot) 

2247 cp_pspot_file = os.path.join(directory, pspot) 

2248 if (os.path.exists(orig_pspot_file) 

2249 and not os.path.exists(cp_pspot_file)): 

2250 if self._copy_pspots: 

2251 shutil.copy(orig_pspot_file, directory) 

2252 elif self._link_pspots: 

2253 os.symlink(orig_pspot_file, cp_pspot_file) 

2254 else: 

2255 if self._pedantic: 

2256 warnings.warn(ppwarning) 

2257 

2258 

2259ppwarning = ('Warning: PP files have neither been ' 

2260 'linked nor copied to the working directory. Make ' 

2261 'sure to set the evironment variable PSPOT_DIR ' 

2262 'accordingly!') 

2263 

2264 

2265def get_castep_version(castep_command): 

2266 """This returns the version number as printed in the CASTEP banner. 

2267 For newer CASTEP versions ( > 6.1) the --version command line option 

2268 has been added; this will be attempted first. 

2269 """ 

2270 import tempfile 

2271 with tempfile.TemporaryDirectory() as temp_dir: 

2272 return _get_castep_version(castep_command, temp_dir) 

2273 

2274 

2275def _get_castep_version(castep_command, temp_dir): 

2276 jname = 'dummy_jobname' 

2277 stdout, stderr = '', '' 

2278 fallback_version = 16. # CASTEP 16.0 and 16.1 report version wrongly 

2279 try: 

2280 stdout, stderr = subprocess.Popen( 

2281 castep_command.split() + ['--version'], 

2282 stderr=subprocess.PIPE, 

2283 stdout=subprocess.PIPE, cwd=temp_dir, 

2284 universal_newlines=True).communicate() 

2285 if 'CASTEP version' not in stdout: 

2286 stdout, stderr = subprocess.Popen( 

2287 castep_command.split() + [jname], 

2288 stderr=subprocess.PIPE, 

2289 stdout=subprocess.PIPE, cwd=temp_dir, 

2290 universal_newlines=True).communicate() 

2291 except Exception: # XXX Which kind of exception? 

2292 msg = '' 

2293 msg += 'Could not determine the version of your CASTEP binary \n' 

2294 msg += 'This usually means one of the following \n' 

2295 msg += ' * you do not have CASTEP installed \n' 

2296 msg += ' * you have not set the CASTEP_COMMAND to call it \n' 

2297 msg += ' * you have provided a wrong CASTEP_COMMAND. \n' 

2298 msg += ' Make sure it is in your PATH\n\n' 

2299 msg += stdout 

2300 msg += stderr 

2301 raise CastepVersionError(msg) 

2302 if 'CASTEP version' in stdout: 

2303 output_txt = stdout.split('\n') 

2304 version_re = re.compile(r'CASTEP version:\s*([0-9\.]*)') 

2305 else: 

2306 output = open(os.path.join(temp_dir, '%s.castep' % jname)) 

2307 output_txt = output.readlines() 

2308 output.close() 

2309 version_re = re.compile(r'(?<=CASTEP version )[0-9.]*') 

2310 # shutil.rmtree(temp_dir) 

2311 for line in output_txt: 

2312 if 'CASTEP version' in line: 

2313 try: 

2314 return float(version_re.findall(line)[0]) 

2315 except ValueError: 

2316 # Fallback for buggy --version on CASTEP 16.0, 16.1 

2317 return fallback_version 

2318 

2319 

2320def create_castep_keywords(castep_command, filename='castep_keywords.json', 

2321 force_write=True, path='.', fetch_only=None): 

2322 """This function allows to fetch all available keywords from stdout 

2323 of an installed castep binary. It furthermore collects the documentation 

2324 to harness the power of (ipython) inspection and type for some basic 

2325 type checking of input. All information is stored in a JSON file that is 

2326 not distributed by default to avoid breaking the license of CASTEP. 

2327 """ 

2328 # Takes a while ... 

2329 # Fetch all allowed parameters 

2330 # fetch_only : only fetch that many parameters (for testsuite only) 

2331 suffixes = ['cell', 'param'] 

2332 

2333 filepath = os.path.join(path, filename) 

2334 

2335 if os.path.exists(filepath) and not force_write: 

2336 warnings.warn('CASTEP Options Module file exists. ' 

2337 'You can overwrite it by calling ' 

2338 'python castep.py -f [CASTEP_COMMAND].') 

2339 return False 

2340 

2341 # Not saving directly to file her to prevent half-generated files 

2342 # which will cause problems on future runs 

2343 

2344 castep_version = get_castep_version(castep_command) 

2345 

2346 help_all, _ = shell_stdouterr('%s -help all' % castep_command) 

2347 

2348 # Filter out proper keywords 

2349 try: 

2350 # The old pattern does not math properly as in CASTEP as of v8.0 there 

2351 # are some keywords for the semi-empircal dispersion correction (SEDC) 

2352 # which also include numbers. 

2353 if castep_version < 7.0: 

2354 pattern = r'((?<=^ )[A-Z_]{2,}|(?<=^)[A-Z_]{2,})' 

2355 else: 

2356 pattern = r'((?<=^ )[A-Z_\d]{2,}|(?<=^)[A-Z_\d]{2,})' 

2357 

2358 raw_options = re.findall(pattern, help_all, re.MULTILINE) 

2359 except Exception: 

2360 warnings.warn('Problem parsing: %s' % help_all) 

2361 raise 

2362 

2363 types = set() 

2364 levels = set() 

2365 

2366 processed_n = 0 

2367 to_process = len(raw_options[:fetch_only]) 

2368 

2369 processed_options = {sf: {} for sf in suffixes} 

2370 

2371 for o_i, option in enumerate(raw_options[:fetch_only]): 

2372 doc, _ = shell_stdouterr('%s -help %s' % (castep_command, option)) 

2373 

2374 # Stand Back! I know regular expressions (http://xkcd.com/208/) :-) 

2375 match = re.match(r'(?P<before_type>.*)Type: (?P<type>.+?)\s+' 

2376 + r'Level: (?P<level>[^ ]+)\n\s*\n' 

2377 + r'(?P<doc>.*?)(\n\s*\n|$)', doc, re.DOTALL) 

2378 

2379 processed_n += 1 

2380 

2381 if match is not None: 

2382 match = match.groupdict() 

2383 

2384 # JM: uncomment lines in following block to debug issues 

2385 # with keyword assignment during extraction process from CASTEP 

2386 suffix = None 

2387 if re.findall(r'PARAMETERS keywords:\n\n\s?None found', doc): 

2388 suffix = 'cell' 

2389 if re.findall(r'CELL keywords:\n\n\s?None found', doc): 

2390 suffix = 'param' 

2391 if suffix is None: 

2392 warnings.warn('%s -> not assigned to either' 

2393 ' CELL or PARAMETERS keywords' % option) 

2394 

2395 option = option.lower() 

2396 mtyp = match.get('type', None) 

2397 mlvl = match.get('level', None) 

2398 mdoc = match.get('doc', None) 

2399 

2400 if mtyp is None: 

2401 warnings.warn('Found no type for %s' % option) 

2402 continue 

2403 if mlvl is None: 

2404 warnings.warn('Found no level for %s' % option) 

2405 continue 

2406 if mdoc is None: 

2407 warnings.warn('Found no doc string for %s' % option) 

2408 continue 

2409 

2410 types = types.union([mtyp]) 

2411 levels = levels.union([mlvl]) 

2412 

2413 processed_options[suffix][option] = { 

2414 'keyword': option, 

2415 'option_type': mtyp, 

2416 'level': mlvl, 

2417 'docstring': mdoc 

2418 } 

2419 

2420 processed_n += 1 

2421 

2422 frac = (o_i + 1.0) / to_process 

2423 sys.stdout.write('\rProcessed: [{0}] {1:>3.0f}%'.format( 

2424 '#' * int(frac * 20) + ' ' 

2425 * (20 - int(frac * 20)), 

2426 100 * frac)) 

2427 sys.stdout.flush() 

2428 

2429 else: 

2430 warnings.warn('create_castep_keywords: Could not process %s' 

2431 % option) 

2432 

2433 sys.stdout.write('\n') 

2434 sys.stdout.flush() 

2435 

2436 processed_options['types'] = list(types) 

2437 processed_options['levels'] = list(levels) 

2438 processed_options['castep_version'] = castep_version 

2439 

2440 json.dump(processed_options, open(filepath, 'w'), indent=4) 

2441 

2442 warnings.warn('CASTEP v%s, fetched %s keywords' % 

2443 (castep_version, processed_n)) 

2444 return True 

2445 

2446 

2447class CastepOption: 

2448 """"A CASTEP option. It handles basic conversions from string to its value 

2449 type.""" 

2450 

2451 default_convert_types = { 

2452 'boolean (logical)': 'bool', 

2453 'defined': 'bool', 

2454 'string': 'str', 

2455 'integer': 'int', 

2456 'real': 'float', 

2457 'integer vector': 'int_vector', 

2458 'real vector': 'float_vector', 

2459 'physical': 'float_physical', 

2460 'block': 'block' 

2461 } 

2462 

2463 def __init__(self, keyword, level, option_type, value=None, 

2464 docstring='No information available'): 

2465 self.keyword = keyword 

2466 self.level = level 

2467 self.type = option_type 

2468 self._value = value 

2469 self.__doc__ = docstring 

2470 

2471 @property 

2472 def value(self): 

2473 

2474 if self._value is not None: 

2475 if self.type.lower() in ('integer vector', 'real vector', 

2476 'physical'): 

2477 return ' '.join(map(str, self._value)) 

2478 elif self.type.lower() in ('boolean (logical)', 'defined'): 

2479 return str(self._value).upper() 

2480 else: 

2481 return str(self._value) 

2482 

2483 @property 

2484 def raw_value(self): 

2485 # The value, not converted to a string 

2486 return self._value 

2487 

2488 @value.setter # type: ignore 

2489 def value(self, val): 

2490 

2491 if val is None: 

2492 self.clear() 

2493 return 

2494 

2495 ctype = self.default_convert_types.get(self.type.lower(), 'str') 

2496 typeparse = '_parse_%s' % ctype 

2497 try: 

2498 self._value = getattr(self, typeparse)(val) 

2499 except ValueError: 

2500 raise ConversionError(ctype, self.keyword, val) 

2501 

2502 def clear(self): 

2503 """Reset the value of the option to None again""" 

2504 self._value = None 

2505 

2506 @staticmethod 

2507 def _parse_bool(value): 

2508 try: 

2509 value = _tf_table[str(value).strip().title()] 

2510 except (KeyError, ValueError): 

2511 raise ValueError() 

2512 return value 

2513 

2514 @staticmethod 

2515 def _parse_str(value): 

2516 value = str(value) 

2517 return value 

2518 

2519 @staticmethod 

2520 def _parse_int(value): 

2521 value = int(value) 

2522 return value 

2523 

2524 @staticmethod 

2525 def _parse_float(value): 

2526 value = float(value) 

2527 return value 

2528 

2529 @staticmethod 

2530 def _parse_int_vector(value): 

2531 # Accepts either a string or an actual list/numpy array of ints 

2532 if isinstance(value, str): 

2533 if ',' in value: 

2534 value = value.replace(',', ' ') 

2535 value = list(map(int, value.split())) 

2536 

2537 value = np.array(value) 

2538 

2539 if value.shape != (3,) or value.dtype != int: 

2540 raise ValueError() 

2541 

2542 return list(value) 

2543 

2544 @staticmethod 

2545 def _parse_float_vector(value): 

2546 # Accepts either a string or an actual list/numpy array of floats 

2547 if isinstance(value, str): 

2548 if ',' in value: 

2549 value = value.replace(',', ' ') 

2550 value = list(map(float, value.split())) 

2551 

2552 value = np.array(value) * 1.0 

2553 

2554 if value.shape != (3,) or value.dtype != float: 

2555 raise ValueError() 

2556 

2557 return list(value) 

2558 

2559 @staticmethod 

2560 def _parse_float_physical(value): 

2561 # If this is a string containing units, saves them 

2562 if isinstance(value, str): 

2563 value = value.split() 

2564 

2565 try: 

2566 l = len(value) 

2567 except TypeError: 

2568 l = 1 

2569 value = [value] 

2570 

2571 if l == 1: 

2572 try: 

2573 value = (float(value[0]), '') 

2574 except (TypeError, ValueError): 

2575 raise ValueError() 

2576 elif l == 2: 

2577 try: 

2578 value = (float(value[0]), value[1]) 

2579 except (TypeError, ValueError, IndexError): 

2580 raise ValueError() 

2581 else: 

2582 raise ValueError() 

2583 

2584 return value 

2585 

2586 @staticmethod 

2587 def _parse_block(value): 

2588 

2589 if isinstance(value, str): 

2590 return value 

2591 elif hasattr(value, '__getitem__'): 

2592 return '\n'.join(value) # Arrays of lines 

2593 else: 

2594 raise ValueError() 

2595 

2596 def __repr__(self): 

2597 if self._value: 

2598 expr = ('Option: {keyword}({type}, {level}):\n{_value}\n' 

2599 ).format(**self.__dict__) 

2600 else: 

2601 expr = ('Option: {keyword}[unset]({type}, {level})' 

2602 ).format(**self.__dict__) 

2603 return expr 

2604 

2605 def __eq__(self, other): 

2606 if not isinstance(other, CastepOption): 

2607 return False 

2608 else: 

2609 return self.__dict__ == other.__dict__ 

2610 

2611 

2612class CastepOptionDict: 

2613 """A dictionary-like object to hold a set of options for .cell or .param 

2614 files loaded from a dictionary, for the sake of validation. 

2615 

2616 Replaces the old CastepCellDict and CastepParamDict that were defined in 

2617 the castep_keywords.py file. 

2618 """ 

2619 

2620 def __init__(self, options=None): 

2621 object.__init__(self) 

2622 self._options = {} # ComparableDict is not needed any more as 

2623 # CastepOptions can be compared directly now 

2624 for kw in options: 

2625 opt = CastepOption(**options[kw]) 

2626 self._options[opt.keyword] = opt 

2627 self.__dict__[opt.keyword] = opt 

2628 

2629 

2630class CastepInputFile: 

2631 

2632 """Master class for CastepParam and CastepCell to inherit from""" 

2633 

2634 _keyword_conflicts: List[Set[str]] = [] 

2635 

2636 def __init__(self, options_dict=None, keyword_tolerance=1): 

2637 object.__init__(self) 

2638 

2639 if options_dict is None: 

2640 options_dict = CastepOptionDict({}) 

2641 

2642 self._options = options_dict._options 

2643 self.__dict__.update(self._options) 

2644 # keyword_tolerance means how strict the checks on new attributes are 

2645 # 0 = no new attributes allowed 

2646 # 1 = new attributes allowed, warning given 

2647 # 2 = new attributes allowed, silent 

2648 self._perm = np.clip(keyword_tolerance, 0, 2) 

2649 

2650 # Compile a dictionary for quick check of conflict sets 

2651 self._conflict_dict = { 

2652 kw: set(cset).difference({kw}) 

2653 for cset in self._keyword_conflicts for kw in cset} 

2654 

2655 def __repr__(self): 

2656 expr = '' 

2657 is_default = True 

2658 for key, option in sorted(self._options.items()): 

2659 if option.value is not None: 

2660 is_default = False 

2661 expr += ('%20s : %s\n' % (key, option.value)) 

2662 if is_default: 

2663 expr = 'Default\n' 

2664 

2665 expr += 'Keyword tolerance: {0}'.format(self._perm) 

2666 return expr 

2667 

2668 def __setattr__(self, attr, value): 

2669 

2670 # Hidden attributes are treated normally 

2671 if attr.startswith('_'): 

2672 self.__dict__[attr] = value 

2673 return 

2674 

2675 if attr not in self._options.keys(): 

2676 

2677 if self._perm > 0: 

2678 # Do we consider it a string or a block? 

2679 is_str = isinstance(value, str) 

2680 is_block = False 

2681 if ((hasattr(value, '__getitem__') and not is_str) 

2682 or (is_str and len(value.split('\n')) > 1)): 

2683 is_block = True 

2684 

2685 if self._perm == 0: 

2686 similars = difflib.get_close_matches(attr, 

2687 self._options.keys()) 

2688 if similars: 

2689 raise UserWarning(('Option "%s" not known! You mean "%s"?') 

2690 % (attr, similars[0])) 

2691 else: 

2692 raise UserWarning('Option "%s" is not known!' % attr) 

2693 elif self._perm == 1: 

2694 warnings.warn(('Option "%s" is not known and will ' 

2695 'be added as a %s') % (attr, 

2696 ('block' if is_block else 

2697 'string'))) 

2698 attr = attr.lower() 

2699 opt = CastepOption(keyword=attr, level='Unknown', 

2700 option_type='block' if is_block else 'string') 

2701 self._options[attr] = opt 

2702 self.__dict__[attr] = opt 

2703 else: 

2704 attr = attr.lower() 

2705 opt = self._options[attr] 

2706 

2707 if not opt.type.lower() == 'block' and isinstance(value, str): 

2708 value = value.replace(':', ' ') 

2709 

2710 # If it is, use the appropriate parser, unless a custom one is defined 

2711 attrparse = '_parse_%s' % attr.lower() 

2712 

2713 # Check for any conflicts if the value is not None 

2714 if not (value is None): 

2715 cset = self._conflict_dict.get(attr.lower(), {}) 

2716 for c in cset: 

2717 if (c in self._options and self._options[c].value): 

2718 warnings.warn( 

2719 'option "{attr}" conflicts with "{conflict}" in ' 

2720 'calculator. Setting "{conflict}" to ' 

2721 'None.'.format(attr=attr, conflict=c)) 

2722 self._options[c].value = None 

2723 

2724 if hasattr(self, attrparse): 

2725 self._options[attr].value = self.__getattribute__(attrparse)(value) 

2726 else: 

2727 self._options[attr].value = value 

2728 

2729 def __getattr__(self, name): 

2730 if name[0] == '_' or self._perm == 0: 

2731 raise AttributeError() 

2732 

2733 if self._perm == 1: 

2734 warnings.warn('Option %s is not known, returning None' % (name)) 

2735 

2736 return CastepOption(keyword='none', level='Unknown', 

2737 option_type='string', value=None) 

2738 

2739 def get_attr_dict(self, raw=False, types=False): 

2740 """Settings that go into .param file in a traditional dict""" 

2741 

2742 attrdict = {k: o.raw_value if raw else o.value 

2743 for k, o in self._options.items() if o.value is not None} 

2744 

2745 if types: 

2746 for key, val in attrdict.items(): 

2747 attrdict[key] = (val, self._options[key].type) 

2748 

2749 return attrdict 

2750 

2751 

2752class CastepParam(CastepInputFile): 

2753 """CastepParam abstracts the settings that go into the .param file""" 

2754 

2755 _keyword_conflicts = [{'cut_off_energy', 'basis_precision'}, ] 

2756 

2757 def __init__(self, castep_keywords, keyword_tolerance=1): 

2758 self._castep_version = castep_keywords.castep_version 

2759 CastepInputFile.__init__(self, castep_keywords.CastepParamDict(), 

2760 keyword_tolerance) 

2761 

2762 @property 

2763 def castep_version(self): 

2764 return self._castep_version 

2765 

2766 # .param specific parsers 

2767 def _parse_reuse(self, value): 

2768 if value is None: 

2769 return None # Reset the value 

2770 try: 

2771 if self._options['continuation'].value: 

2772 warnings.warn('Cannot set reuse if continuation is set, and ' 

2773 'vice versa. Set the other to None, if you want ' 

2774 'this setting.') 

2775 return None 

2776 except KeyError: 

2777 pass 

2778 return 'default' if (value is True) else str(value) 

2779 

2780 def _parse_continuation(self, value): 

2781 if value is None: 

2782 return None # Reset the value 

2783 try: 

2784 if self._options['reuse'].value: 

2785 warnings.warn('Cannot set reuse if continuation is set, and ' 

2786 'vice versa. Set the other to None, if you want ' 

2787 'this setting.') 

2788 return None 

2789 except KeyError: 

2790 pass 

2791 return 'default' if (value is True) else str(value) 

2792 

2793 

2794class CastepCell(CastepInputFile): 

2795 

2796 """CastepCell abstracts all setting that go into the .cell file""" 

2797 

2798 _keyword_conflicts = [ 

2799 {'kpoint_mp_grid', 'kpoint_mp_spacing', 'kpoint_list', 

2800 'kpoints_mp_grid', 'kpoints_mp_spacing', 'kpoints_list'}, 

2801 {'bs_kpoint_mp_grid', 

2802 'bs_kpoint_mp_spacing', 

2803 'bs_kpoint_list', 

2804 'bs_kpoint_path', 

2805 'bs_kpoints_mp_grid', 

2806 'bs_kpoints_mp_spacing', 

2807 'bs_kpoints_list', 

2808 'bs_kpoints_path'}, 

2809 {'spectral_kpoint_mp_grid', 

2810 'spectral_kpoint_mp_spacing', 

2811 'spectral_kpoint_list', 

2812 'spectral_kpoint_path', 

2813 'spectral_kpoints_mp_grid', 

2814 'spectral_kpoints_mp_spacing', 

2815 'spectral_kpoints_list', 

2816 'spectral_kpoints_path'}, 

2817 {'phonon_kpoint_mp_grid', 

2818 'phonon_kpoint_mp_spacing', 

2819 'phonon_kpoint_list', 

2820 'phonon_kpoint_path', 

2821 'phonon_kpoints_mp_grid', 

2822 'phonon_kpoints_mp_spacing', 

2823 'phonon_kpoints_list', 

2824 'phonon_kpoints_path'}, 

2825 {'fine_phonon_kpoint_mp_grid', 

2826 'fine_phonon_kpoint_mp_spacing', 

2827 'fine_phonon_kpoint_list', 

2828 'fine_phonon_kpoint_path'}, 

2829 {'magres_kpoint_mp_grid', 

2830 'magres_kpoint_mp_spacing', 

2831 'magres_kpoint_list', 

2832 'magres_kpoint_path'}, 

2833 {'elnes_kpoint_mp_grid', 

2834 'elnes_kpoint_mp_spacing', 

2835 'elnes_kpoint_list', 

2836 'elnes_kpoint_path'}, 

2837 {'optics_kpoint_mp_grid', 

2838 'optics_kpoint_mp_spacing', 

2839 'optics_kpoint_list', 

2840 'optics_kpoint_path'}, 

2841 {'supercell_kpoint_mp_grid', 

2842 'supercell_kpoint_mp_spacing', 

2843 'supercell_kpoint_list', 

2844 'supercell_kpoint_path'}, ] 

2845 

2846 def __init__(self, castep_keywords, keyword_tolerance=1): 

2847 self._castep_version = castep_keywords.castep_version 

2848 CastepInputFile.__init__(self, castep_keywords.CastepCellDict(), 

2849 keyword_tolerance) 

2850 

2851 @property 

2852 def castep_version(self): 

2853 return self._castep_version 

2854 

2855 # .cell specific parsers 

2856 def _parse_species_pot(self, value): 

2857 

2858 # Single tuple 

2859 if isinstance(value, tuple) and len(value) == 2: 

2860 value = [value] 

2861 # List of tuples 

2862 if hasattr(value, '__getitem__'): 

2863 pspots = [tuple(map(str.strip, x)) for x in value] 

2864 if not all(map(lambda x: len(x) == 2, value)): 

2865 warnings.warn( 

2866 'Please specify pseudopotentials in python as ' 

2867 'a tuple or a list of tuples formatted like: ' 

2868 '(species, file), e.g. ("O", "path-to/O_OTFG.usp") ' 

2869 'Anything else will be ignored') 

2870 return None 

2871 

2872 text_block = self._options['species_pot'].value 

2873 

2874 text_block = text_block if text_block else '' 

2875 # Remove any duplicates 

2876 for pp in pspots: 

2877 text_block = re.sub(r'\n?\s*%s\s+.*' % pp[0], '', text_block) 

2878 if pp[1]: 

2879 text_block += '\n%s %s' % pp 

2880 

2881 return text_block 

2882 

2883 def _parse_symmetry_ops(self, value): 

2884 if not isinstance(value, tuple) \ 

2885 or not len(value) == 2 \ 

2886 or not value[0].shape[1:] == (3, 3) \ 

2887 or not value[1].shape[1:] == (3,) \ 

2888 or not value[0].shape[0] == value[1].shape[0]: 

2889 warnings.warn('Invalid symmetry_ops block, skipping') 

2890 return 

2891 # Now on to print... 

2892 text_block = '' 

2893 for op_i, (op_rot, op_tranls) in enumerate(zip(*value)): 

2894 text_block += '\n'.join([' '.join([str(x) for x in row]) 

2895 for row in op_rot]) 

2896 text_block += '\n' 

2897 text_block += ' '.join([str(x) for x in op_tranls]) 

2898 text_block += '\n\n' 

2899 

2900 return text_block 

2901 

2902 def _parse_positions_abs_intermediate(self, value): 

2903 return _parse_tss_block(value) 

2904 

2905 def _parse_positions_abs_product(self, value): 

2906 return _parse_tss_block(value) 

2907 

2908 def _parse_positions_frac_intermediate(self, value): 

2909 return _parse_tss_block(value, True) 

2910 

2911 def _parse_positions_frac_product(self, value): 

2912 return _parse_tss_block(value, True) 

2913 

2914 

2915CastepKeywords = namedtuple('CastepKeywords', 

2916 ['CastepParamDict', 'CastepCellDict', 

2917 'types', 'levels', 'castep_version']) 

2918 

2919# We keep this just for naming consistency with older versions 

2920 

2921 

2922def make_cell_dict(data=None): 

2923 

2924 data = data if data is not None else {} 

2925 

2926 class CastepCellDict(CastepOptionDict): 

2927 def __init__(self): 

2928 CastepOptionDict.__init__(self, data) 

2929 

2930 return CastepCellDict 

2931 

2932 

2933def make_param_dict(data=None): 

2934 

2935 data = data if data is not None else {} 

2936 

2937 class CastepParamDict(CastepOptionDict): 

2938 def __init__(self): 

2939 CastepOptionDict.__init__(self, data) 

2940 

2941 return CastepParamDict 

2942 

2943 

2944class CastepVersionError(Exception): 

2945 """No special behaviour, works to signal when Castep can not be found""" 

2946 pass 

2947 

2948 

2949class ConversionError(Exception): 

2950 

2951 """Print customized error for options that are not converted correctly 

2952 and point out that they are maybe not implemented, yet""" 

2953 

2954 def __init__(self, key_type, attr, value): 

2955 Exception.__init__(self) 

2956 self.key_type = key_type 

2957 self.value = value 

2958 self.attr = attr 

2959 

2960 def __str__(self): 

2961 return 'Could not convert %s = %s to %s\n' \ 

2962 % (self.attr, self.value, self.key_type) \ 

2963 + 'This means you either tried to set a value of the wrong\n'\ 

2964 + 'type or this keyword needs some special care. Please feel\n'\ 

2965 + 'to add it to the corresponding __setattr__ method and send\n'\ 

2966 + 'the patch to %s, so we can all benefit.' % (contact_email) 

2967 

2968 

2969def get_castep_pp_path(castep_pp_path=''): 

2970 """Abstract the quest for a CASTEP PSP directory.""" 

2971 if castep_pp_path: 

2972 return os.path.abspath(os.path.expanduser(castep_pp_path)) 

2973 elif 'PSPOT_DIR' in os.environ: 

2974 return os.environ['PSPOT_DIR'] 

2975 elif 'CASTEP_PP_PATH' in os.environ: 

2976 return os.environ['CASTEP_PP_PATH'] 

2977 else: 

2978 return os.path.abspath('.') 

2979 

2980 

2981def get_castep_command(castep_command=''): 

2982 """Abstract the quest for a castep_command string.""" 

2983 if castep_command: 

2984 return castep_command 

2985 elif 'CASTEP_COMMAND' in os.environ: 

2986 return os.environ['CASTEP_COMMAND'] 

2987 else: 

2988 return 'castep' 

2989 

2990 

2991def shell_stdouterr(raw_command, cwd=None): 

2992 """Abstracts the standard call of the commandline, when 

2993 we are only interested in the stdout and stderr 

2994 """ 

2995 stdout, stderr = subprocess.Popen(raw_command, 

2996 stdout=subprocess.PIPE, 

2997 stderr=subprocess.PIPE, 

2998 universal_newlines=True, 

2999 shell=True, cwd=cwd).communicate() 

3000 return stdout.strip(), stderr.strip() 

3001 

3002 

3003def import_castep_keywords(castep_command='', 

3004 filename='castep_keywords.json', 

3005 path='.'): 

3006 """Search for castep keywords JSON in multiple paths""" 

3007 

3008 config_paths = ('~/.ase', '~/.config/ase') 

3009 searchpaths = [path] + [os.path.expanduser(config_path) 

3010 for config_path in config_paths] 

3011 try: 

3012 keywords_file = sum([glob.glob(os.path.join(sp, filename)) 

3013 for sp in searchpaths], [])[0] 

3014 except IndexError: 

3015 warnings.warn("""Generating CASTEP keywords JSON file... hang on. 

3016 The CASTEP keywords JSON file contains abstractions for CASTEP input 

3017 parameters (for both .cell and .param input files), including some 

3018 format checks and descriptions. The latter are extracted from the 

3019 internal online help facility of a CASTEP binary, thus allowing to 

3020 easily keep the calculator synchronized with (different versions of) 

3021 the CASTEP code. Consequently, avoiding licensing issues (CASTEP is 

3022 distributed commercially by Biovia), we consider it wise not to 

3023 provide the file in the first place.""") 

3024 create_castep_keywords(get_castep_command(castep_command), 

3025 filename=filename, path=path) 

3026 keywords_file = Path(path).absolute() / filename 

3027 

3028 warnings.warn( 

3029 f'Stored castep keywords dictionary as {keywords_file}. ' 

3030 f'Copy it to {Path(config_paths[0]).expanduser() / filename} for ' 

3031 r'user installation.') 

3032 

3033 # Now create the castep_keywords object proper 

3034 with open(keywords_file) as fd: 

3035 kwdata = json.load(fd) 

3036 

3037 # This is a bit awkward, but it's necessary for backwards compatibility 

3038 param_dict = make_param_dict(kwdata['param']) 

3039 cell_dict = make_cell_dict(kwdata['cell']) 

3040 

3041 castep_keywords = CastepKeywords(param_dict, cell_dict, 

3042 kwdata['types'], kwdata['levels'], 

3043 kwdata['castep_version']) 

3044 

3045 return castep_keywords 

3046 

3047 

3048if __name__ == '__main__': 

3049 warnings.warn( 

3050 'When called directly this calculator will fetch all available ' 

3051 'keywords from the binarys help function into a ' 

3052 'castep_keywords.json in the current directory %s ' 

3053 'For system wide usage, it can be copied into an ase installation ' 

3054 'at ASE/calculators. ' 

3055 'This castep_keywords.json usually only needs to be generated once ' 

3056 'for a CASTEP binary/CASTEP version.' % os.getcwd()) 

3057 

3058 import optparse 

3059 parser = optparse.OptionParser() 

3060 parser.add_option( 

3061 '-f', '--force-write', dest='force_write', 

3062 help='Force overwriting existing castep_keywords.json', default=False, 

3063 action='store_true') 

3064 (options, args) = parser.parse_args() 

3065 

3066 if args: 

3067 opt_castep_command = ''.join(args) 

3068 else: 

3069 opt_castep_command = '' 

3070 generated = create_castep_keywords(get_castep_command(opt_castep_command), 

3071 force_write=options.force_write) 

3072 

3073 if generated: 

3074 try: 

3075 with open('castep_keywords.json') as fd: 

3076 json.load(fd) 

3077 except Exception as e: 

3078 warnings.warn( 

3079 '%s Ooops, something went wrong with the CASTEP keywords' % e) 

3080 else: 

3081 warnings.warn('Import works. Looking good!')