Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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 typing import List, Set 

32 

33import ase 

34import ase.units as units 

35from ase.calculators.general import Calculator 

36from ase.calculators.calculator import compare_atoms 

37from ase.calculators.calculator import PropertyNotImplementedError 

38from ase.calculators.calculator import kpts2sizeandoffsets 

39from ase.dft.kpoints import BandPath 

40from ase.parallel import paropen 

41from ase.io.castep import read_param 

42from ase.io.castep import read_bands 

43from ase.constraints import FixCartesian 

44 

45__all__ = [ 

46 'Castep', 

47 'CastepCell', 

48 'CastepParam', 

49 'create_castep_keywords'] 

50 

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

52 

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

54_tf_table = { 

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

56 'True': True, 

57 'False': False} 

58 

59 

60def _self_getter(getf): 

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

62 # getter function, self.atoms is used instead 

63 

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

65 

66 if atoms is None: 

67 atoms = self.atoms 

68 

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

70 

71 return decor_getf 

72 

73 

74def _parse_tss_block(value, scaled=False): 

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

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

77 try: 

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

79 except TypeError: 

80 is_strlist = False 

81 

82 if not is_atoms: 

83 if not is_strlist: 

84 # Invalid! 

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

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

87 

88 # First line must be Angstroms, or nothing 

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

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

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

92 'cell.positions_abs_intermediate/product') 

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

94 else: 

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

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

97 value.get_positions()) 

98 symbols = value.get_chemical_symbols() 

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

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

101 

102 return text_block 

103 

104 

105class Castep(Calculator): 

106 r""" 

107CASTEP Interface Documentation 

108 

109 

110Introduction 

111============ 

112 

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

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

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

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

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

118molecular dynamics simulations. 

119 

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

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

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

123``calc.cell.keyword = ...``) 

124 

125 

126Getting Started: 

127================ 

128 

129Set the environment variables appropriately for your system. 

130 

131>>> export CASTEP_COMMAND=' ... ' 

132>>> export CASTEP_PP_PATH=' ... ' 

133 

134Note: alternatively to CASTEP_PP_PATH one can set PSPOT_DIR 

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

136 

137>>> export PSPOT_DIR=' ... ' 

138 

139 

140Running the Calculator 

141====================== 

142 

143The default initialization command for the CASTEP calculator is 

144 

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

146 

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

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

149 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

164 

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

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

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

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

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

170 

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

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

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

174the atoms object. 

175 

176 

177Arguments: 

178========== 

179 

180========================= ==================================================== 

181Keyword Description 

182========================= ==================================================== 

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

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

185 created. Existing directories will be moved to 

186 directory-TIMESTAMP unless self._rename_existing_dir 

187 is set to false. 

188 

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

190 

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

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

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

194 

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

196 version matches the version from which the available 

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

198 

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

200 also be set via the bash environment variables 

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

202 Will default to the current working directory if 

203 none is given or found. Note that pseudopotentials 

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

205 

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

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

208 this directory will be checked for typical names. If 

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

210 fly, depending on the ``_build_missing_pspots`` 

211 value. A RuntimeError will be raised in case 

212 multiple files per element are found. Defaults to 

213 ``False``. 

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

215 validation of any parameters set in the CastepCell 

216 or CastepParam objects against the ones found in 

217 castep_keywords. Levels are as following: 

218 

219 0 = no tolerance, keywords not found in 

220 castep_keywords will raise an exception 

221 

222 1 = keywords not found will be accepted but produce 

223 a warning (default) 

224 

225 2 = keywords not found will be accepted silently 

226 

227 3 = no attempt is made to look for 

228 castep_keywords.json at all 

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

230 then used with no attempt to actually load a 

231 castep_keywords.json file. Most useful for debugging 

232 and testing purposes. 

233 

234========================= ==================================================== 

235 

236 

237Additional Settings 

238=================== 

239 

240========================= ==================================================== 

241Internal Setting Description 

242========================= ==================================================== 

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

244 call CASTEP. 

245 

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

247 write a continue or reuse statement if the 

248 addressed .check or .castep_bin file exists in the 

249 directory. 

250 

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

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

253 file, usually it will only create symlinks. 

254 

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

256 actually will create symlinks to the needed pseudo 

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

258 to False if you rather want to access your pseudo 

259 potentials using the PSPOT_DIR environment variable 

260 that is read by CASTEP. 

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

262 is True.. 

263 

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

265 missing pseudopotentials on the fly. If not, a 

266 RuntimeError will be raised if not all files were 

267 found. 

268 

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

270 True, all calculator internal settings shown here 

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

272 and can be read again by merge_param. merge_param 

273 can be forced to ignore this directive using the 

274 optional argument ``ignore_internal_keys=True``. 

275 

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

277 \*param will be overwritten. 

278 

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

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

281 start the calculation itself. 

282 If this is used to prepare jobs locally 

283 and run on a remote cluster it is recommended 

284 to set ``_copy_pspots = True``. 

285 

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

287 will look for pseudo-potential files. 

288 

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

290 try to find the respective pseudopotentials from 

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

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

293 feature should be very robust. Raises a RuntimeError 

294 if required files are not unique (multiple files per 

295 element). Non existing pseudopotentials will be 

296 generated, though this could be dangerous. 

297 

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

299 of the calculator, this will move directories out of 

300 the way that would be overwritten otherwise, 

301 appending a date string. 

302 

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

304 any atoms object previously attached to the 

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

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

307 atoms object if none has been attached and other- 

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

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

310 if one uses CASTEP's internal geometry optimization 

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

312 because then the positions get out of sync. 

313 *Warning*: this option is generally not recommended 

314 unless one knows one really needs it. There should 

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

316 single-point calculator. 

317 

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

319 will append a number to the label on all input 

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

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

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

323 nal \*check files. 

324 

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

326 terface will try to fetch the reuse file from the 

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

328 fault it is equal to _track_output, but may be 

329 overridden. 

330 

331 Since this behavior may not always be desirable for 

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

333 a geometry-optimization can be achieved by setting 

334 ``calc.param.reuse = True``. 

335 

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

337 inform about settings probably wasting a lot of CPU 

338 time or causing numerical inconsistencies. 

339 

340========================= ==================================================== 

341 

342Special features: 

343================= 

344 

345 

346``.dryrun_ok()`` 

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

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

349 input. Afterwards _kpoints_used is set. 

350 

351``.merge_param()`` 

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

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

354 

355``.keyword.clear()`` 

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

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

358 

359``.initialize()`` 

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

361 run in a place without ASE or even python. 

362 

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

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

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

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

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

368 

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

370 This automatically searches for pseudopotentials of type 

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

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

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

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

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

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

377 work with the defaults. 

378 

379``print(calc)`` 

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

381 

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

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

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

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

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

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

388 avoid overwriting. 

389 

390``.set_kpts(kpoints)`` 

391 This is equivalent to initialising the calculator with 

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

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

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

395 given in reciprocal lattice coordinates e.g. 

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

397 available for more complex requirements e.g. 

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

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

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

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

402 point. 

403 

404``.set_bandpath(bandpath)`` 

405 This is equivalent to initialialising the calculator with 

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

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

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

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

410 

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

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

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

414 ``calc.band_structure().plot()`` 

415 

416Notes/Issues: 

417============== 

418 

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

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

421 

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

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

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

425 used for constants and conversion factors, respectively. 

426 

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

428 

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

430 

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

432 

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

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

435 pp.567- 570 (2005) PDF_. 

436 

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

438 

439 

440End CASTEP Interface Documentation 

441 """ 

442 

443 # Class attributes ! 

444 # keys set through atoms object 

445 atoms_keys = [ 

446 'charges', 

447 'ionic_constraints', 

448 'lattice_abs', 

449 'lattice_cart', 

450 'positions_abs', 

451 'positions_abs_final', 

452 'positions_abs_intermediate', 

453 'positions_frac', 

454 'positions_frac_final', 

455 'positions_frac_intermediate'] 

456 

457 atoms_obj_keys = [ 

458 'dipole', 

459 'energy_free', 

460 'energy_zero', 

461 'fermi', 

462 'forces', 

463 'nbands', 

464 'positions', 

465 'stress', 

466 'pressure'] 

467 

468 internal_keys = [ 

469 '_castep_command', 

470 '_check_checkfile', 

471 '_copy_pspots', 

472 '_link_pspots', 

473 '_find_pspots', 

474 '_build_missing_pspots', 

475 '_directory', 

476 '_export_settings', 

477 '_force_write', 

478 '_label', 

479 '_prepare_input_only', 

480 '_castep_pp_path', 

481 '_rename_existing_dir', 

482 '_set_atoms', 

483 '_track_output', 

484 '_try_reuse', 

485 '_pedantic'] 

486 

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

488 castep_command=None, check_castep_version=False, 

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

490 castep_keywords=None, **kwargs): 

491 

492 self.__name__ = 'Castep' 

493 

494 # initialize the ase.calculators.general calculator 

495 Calculator.__init__(self) 

496 

497 from ase.io.castep import write_cell 

498 self._write_cell = write_cell 

499 

500 if castep_keywords is None: 

501 castep_keywords = CastepKeywords(make_param_dict(), 

502 make_cell_dict(), 

503 [], 

504 [], 

505 0) 

506 if keyword_tolerance < 3: 

507 try: 

508 castep_keywords = import_castep_keywords(castep_command) 

509 except CastepVersionError as e: 

510 if keyword_tolerance == 0: 

511 raise e 

512 else: 

513 warnings.warn(str(e)) 

514 

515 self._kw_tol = keyword_tolerance 

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

517 self.param = CastepParam(castep_keywords, 

518 keyword_tolerance=keyword_tolerance) 

519 self.cell = CastepCell(castep_keywords, 

520 keyword_tolerance=keyword_tolerance) 

521 

522 ################################### 

523 # Calculator state variables # 

524 ################################### 

525 self._calls = 0 

526 self._castep_version = castep_keywords.castep_version 

527 

528 # collects warning from .castep files 

529 self._warnings = [] 

530 # collects content from *.err file 

531 self._error = None 

532 # warnings raised by the ASE interface 

533 self._interface_warnings = [] 

534 

535 # store to check if recalculation is necessary 

536 self._old_atoms = None 

537 self._old_cell = None 

538 self._old_param = None 

539 

540 ################################### 

541 # Internal keys # 

542 # Allow to tweak the behavior # 

543 ################################### 

544 self._opt = {} 

545 self._castep_command = get_castep_command(castep_command) 

546 self._castep_pp_path = get_castep_pp_path(castep_pp_path) 

547 self._check_checkfile = True 

548 self._copy_pspots = False 

549 self._link_pspots = True 

550 self._find_pspots = find_pspots 

551 self._build_missing_pspots = True 

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

553 self._export_settings = True 

554 self._force_write = True 

555 self._label = label 

556 self._prepare_input_only = False 

557 self._rename_existing_dir = True 

558 self._set_atoms = False 

559 self._track_output = False 

560 self._try_reuse = False 

561 

562 # turn off the pedantic user warnings 

563 self._pedantic = False 

564 

565 # will be set on during runtime 

566 self._seed = None 

567 

568 ################################### 

569 # (Physical) result variables # 

570 ################################### 

571 self.atoms = None 

572 # initialize result variables 

573 self._forces = None 

574 self._energy_total = None 

575 self._energy_free = None 

576 self._energy_0K = None 

577 self._energy_total_corr = None 

578 self._eigenvalues = None 

579 self._efermi = None 

580 self._ibz_kpts = None 

581 self._ibz_weights = None 

582 self._band_structure = None 

583 

584 # dispersion corrections 

585 self._dispcorr_energy_total = None 

586 self._dispcorr_energy_free = None 

587 self._dispcorr_energy_0K = None 

588 

589 # spins and hirshfeld volumes 

590 self._spins = None 

591 self._hirsh_volrat = None 

592 

593 # Mulliken charges 

594 self._mulliken_charges = None 

595 # Hirshfeld charges 

596 self._hirshfeld_charges = None 

597 

598 self._number_of_cell_constraints = None 

599 self._output_verbosity = None 

600 self._stress = None 

601 self._pressure = None 

602 self._unit_cell = None 

603 self._kpoints = None 

604 

605 # pointers to other files used at runtime 

606 self._check_file = None 

607 self._castep_bin_file = None 

608 

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

610 self._cut_off_energy = None 

611 

612 # runtime information 

613 self._total_time = None 

614 self._peak_memory = None 

615 

616 # check version of CASTEP options module against current one 

617 if check_castep_version: 

618 local_castep_version = get_castep_version(self._castep_command) 

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

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

621 return 

622 if not local_castep_version == self._castep_version: 

623 warnings.warn( 

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

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

626 (self._castep_version, 

627 get_castep_version(self._castep_command))) 

628 self._castep_version = local_castep_version 

629 

630 # processes optional arguments in kw style 

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

632 # first fetch special keywords issued by ASE CLI 

633 if keyword == 'kpts': 

634 self.set_kpts(value) 

635 elif keyword == 'bandpath': 

636 self.set_bandpath(value) 

637 elif keyword == 'xc': 

638 self.xc_functional = value 

639 elif keyword == 'ecut': 

640 self.cut_off_energy = value 

641 else: # the general case 

642 self.__setattr__(keyword, value) 

643 

644 def band_structure(self, bandfile=None): 

645 from ase.spectrum.band_structure import BandStructure 

646 

647 if bandfile is None: 

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

649 

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

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

652 

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

654 

655 # Get definitions of high-symmetry points 

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

657 bandpath = BandPath(self.atoms.cell, 

658 kpts=kpts, 

659 special_points=special_points) 

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

661 

662 def set_bandpath(self, bandpath): 

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

664 

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

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

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

668 obtain a new dense path. 

669 

670 Args: 

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

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

673 sampling will follow BandPath parameters. 

674 

675 """ 

676 

677 def clear_bs_keywords(): 

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

679 {'path', 'path_spacing', 

680 'list', 

681 'mp_grid', 'mp_spacing', 'mp_offset'}) 

682 for bs_tag in bs_keywords: 

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

684 

685 if bandpath is None: 

686 clear_bs_keywords() 

687 elif isinstance(bandpath, BandPath): 

688 clear_bs_keywords() 

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

690 for row in bandpath.kpts] 

691 else: 

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

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

694 

695 def set_kpts(self, kpts): 

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

697 

698 Args: 

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

700 

701 This method will set the CASTEP parameters kpoints_mp_grid, 

702 kpoints_mp_offset and kpoints_mp_spacing as appropriate. Unused 

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

704 

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

706 

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

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

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

710 input files. 

711 

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

713 dictionary with the following allowed keys: 

714 

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

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

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

718 be set to allow for rounding/centering. 

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

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

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

722 'density' to allow for rounding/centering. 

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

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

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

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

727 """ 

728 

729 def clear_mp_keywords(): 

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

731 {'mp_grid', 'mp_offset', 

732 'mp_spacing', 'list'}) 

733 for kp_tag in mp_keywords: 

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

735 

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

737 if kpts is None: 

738 clear_mp_keywords() 

739 pass 

740 

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

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

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

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

745 # [-0.5, -0.5, -0.5, 0.125]] 

746 

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

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

749 

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

751 raise ValueError( 

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

753 

754 clear_mp_keywords() 

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

756 

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

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

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

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

761 

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

763 raise ValueError( 

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

765 

766 clear_mp_keywords() 

767 self.cell.kpoint_list = kpts 

768 

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

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

771 if len(kpts) != 3: 

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

773 clear_mp_keywords() 

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

775 

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

777 elif isinstance(kpts, str): 

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

779 

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

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

782 elif isinstance(kpts, dict): 

783 kpts = kpts.copy() 

784 

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

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

787 raise ValueError( 

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

789 else: 

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

791 kpts = kpts.copy() 

792 spacing = kpts.pop('spacing') 

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

794 

795 clear_mp_keywords() 

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

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

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

799 

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

801 elif hasattr(kpts, '__iter__'): 

802 self.set_kpts(list(kpts)) 

803 

804 # Otherwise, give up 

805 else: 

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

807 

808 def todict(self, skip_default=True): 

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

810 dct = {} 

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

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

813 

814 return dct 

815 

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

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

818 return compare_atoms(self._old_atoms, atoms) 

819 

820 def _castep_find_last_record(self, castep_file): 

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

822 ending message following the last banner message. If this 

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

824 is return, otherwise False. 

825 

826 returns (record_start, record_end, end_found, last_record_complete) 

827 """ 

828 if isinstance(castep_file, str): 

829 castep_file = paropen(castep_file, 'r') 

830 file_opened = True 

831 else: 

832 file_opened = False 

833 record_starts = [] 

834 while True: 

835 line = castep_file.readline() 

836 if 'Welcome' in line and 'CASTEP' in line: 

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

838 if not line: 

839 break 

840 

841 if record_starts == []: 

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

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

844 return 

845 

846 # search for regular end of file 

847 end_found = False 

848 # start to search from record beginning from the back 

849 # and see if 

850 record_end = -1 

851 for record_nr, record_start in enumerate(record_starts): 

852 castep_file.seek(record_start) 

853 while True: 

854 line = castep_file.readline() 

855 if not line: 

856 break 

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

858 self._warnings.append(line) 

859 

860 if 'Finalisation time =' in line: 

861 end_found = True 

862 record_end = castep_file.tell() 

863 break 

864 

865 if end_found: 

866 break 

867 

868 if file_opened: 

869 castep_file.close() 

870 

871 if end_found: 

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

873 if record_nr == 0: 

874 return (record_start, record_end, True, True) 

875 else: 

876 return (record_start, record_end, True, False) 

877 else: 

878 return (0, record_end, False, False) 

879 

880 def read(self, castep_file=None): 

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

882 

883 _close = True 

884 

885 if castep_file is None: 

886 if self._castep_file: 

887 castep_file = self._castep_file 

888 out = paropen(castep_file, 'r') 

889 else: 

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

891 return 

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

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

894 

895 elif isinstance(castep_file, str): 

896 out = paropen(castep_file, 'r') 

897 

898 else: 

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

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

901 out = castep_file 

902 

903 # look before you leap... 

904 attributes = ['name', 

905 'seek', 

906 'close', 

907 'readline', 

908 'tell'] 

909 

910 for attr in attributes: 

911 if not hasattr(out, attr): 

912 raise TypeError( 

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

914 

915 castep_file = out.name 

916 _close = False 

917 

918 if self._seed is None: 

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

920 

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

922 if os.path.exists(err_file): 

923 err_file = paropen(err_file) 

924 self._error = err_file.read() 

925 err_file.close() 

926 # we return right-away because it might 

927 # just be here from a previous run 

928 # look for last result, if several CASTEP 

929 # run are appended 

930 

931 record_start, record_end, end_found, _\ 

932 = self._castep_find_last_record(out) 

933 if not end_found: 

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

935 (castep_file, self._error)) 

936 if _close: 

937 out.close() 

938 return 

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

940 

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

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

943 # holding several outputs 

944 n_cell_const = 0 

945 forces = [] 

946 

947 # HOTFIX: 

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

949 # otherwise the calculator crashes upon pickling trajectories 

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

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

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

953 # stress = [] 

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

955 hirsh_volrat = [] 

956 

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

958 # Hirshfeld volumes are calculated 

959 spin_polarized = False 

960 calculate_hirshfeld = False 

961 mulliken_analysis = False 

962 hirshfeld_analysis = False 

963 kpoints = None 

964 

965 positions_frac_list = [] 

966 mulliken_charges = [] 

967 spins = [] 

968 

969 out.seek(record_start) 

970 while True: 

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

972 # atoms objects for intermediate steps. 

973 try: 

974 # in case we need to rewind back one line, we memorize the bit 

975 # position of this line in the file. 

976 # --> see symops problem below 

977 _line_start = out.tell() 

978 line = out.readline() 

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

980 break 

981 elif 'Hirshfeld Analysis' in line: 

982 hirshfeld_charges = [] 

983 

984 hirshfeld_analysis = True 

985 # skip the separating line 

986 line = out.readline() 

987 # this is the headline 

988 line = out.readline() 

989 

990 if 'Charge' in line: 

991 # skip the next separator line 

992 line = out.readline() 

993 while True: 

994 line = out.readline() 

995 fields = line.split() 

996 if len(fields) == 1: 

997 break 

998 else: 

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

1000 elif 'stress calculation' in line: 

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

1002 self.param.calculate_stress = True 

1003 elif 'basis set accuracy' in line: 

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

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

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

1007 # conflict with input option basis_precision 

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

1009 self._cut_off_energy = cutoff 

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

1011 self.param.cut_off_energy = cutoff 

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

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

1014 self.param.elec_energy_tol = elec_energy_tol 

1015 elif 'convergence tolerance window' in line: 

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

1017 self.param.elec_convergence_win = elec_convergence_win 

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

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

1020 fbc_possibilities = {'none': 0, 

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

1022 fbc = fbc_possibilities[finite_basis_corr] 

1023 self.param.finite_basis_corr = fbc 

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

1025 self.param.fix_occupancy = True 

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

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

1028 self.param.max_scf_cycles = max_no_scf 

1029 elif 'density-mixing scheme' in line: 

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

1031 self.param.mixing_scheme = mixing_scheme 

1032 elif 'dump wavefunctions every' in line: 

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

1034 self.param.num_dump_cycles = no_dump_cycles 

1035 elif 'optimization strategy' in line: 

1036 lspl = line.split(":") 

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

1038 # This can happen in iprint: 3 calculations 

1039 continue 

1040 if 'memory' in line: 

1041 self.param.opt_strategy = 'Memory' 

1042 if 'speed' in line: 

1043 self.param.opt_strategy = 'Speed' 

1044 elif 'calculation limited to maximum' in line: 

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

1046 self.param.run_time = calc_limit 

1047 elif 'type of calculation' in line: 

1048 lspl = line.split(":") 

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

1050 # This can happen in iprint: 3 calculations 

1051 continue 

1052 calc_type = lspl[-1] 

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

1054 calc_type = calc_type.strip() 

1055 if calc_type != 'single point energy': 

1056 calc_type_possibilities = { 

1057 'geometry optimization': 'GeometryOptimization', 

1058 'band structure': 'BandStructure', 

1059 'molecular dynamics': 'MolecularDynamics', 

1060 'optical properties': 'Optics', 

1061 'phonon calculation': 'Phonon', 

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

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

1064 'transition state search': 'TransitionStateSearch', 

1065 'Magnetic Resonance': 'MagRes', 

1066 'Core level spectra': 'Elnes', 

1067 'Electronic Spectroscopy': 'ElectronicSpectroscopy' 

1068 } 

1069 ctype = calc_type_possibilities[calc_type] 

1070 self.param.task = ctype 

1071 elif 'using functional' in line: 

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

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

1074 used_functional = used_functional.strip() 

1075 if used_functional != 'Local Density Approximation': 

1076 used_functional_possibilities = { 

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

1078 'Perdew Burke Ernzerhof': 'PBE', 

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

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

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

1082 'Hartree-Fock': 'HF', 

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

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

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

1086 'hybrid PBE0': 'PBE0', 

1087 'hybrid B3LYP': 'B3LYP', 

1088 'hybrid HSE03': 'HSE03', 

1089 'hybrid HSE06': 'HSE06' 

1090 } 

1091 used_func = used_functional_possibilities[ 

1092 used_functional] 

1093 self.param.xc_functional = used_func 

1094 elif 'output verbosity' in line: 

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

1096 if int(iprint) != 1: 

1097 self.param.iprint = iprint 

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

1099 spin_polarized = True 

1100 self.param.spin_polarized = spin_polarized 

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

1102 spin_polarized = False 

1103 elif 'Number of kpoints used' in line: 

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

1105 elif 'Unit Cell' in line: 

1106 lattice_real = [] 

1107 lattice_reci = [] 

1108 while True: 

1109 line = out.readline() 

1110 fields = line.split() 

1111 if len(fields) == 6: 

1112 break 

1113 for i in range(3): 

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

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

1116 line = out.readline() 

1117 fields = line.split() 

1118 elif 'Cell Contents' in line: 

1119 while True: 

1120 line = out.readline() 

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

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

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

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

1125 fields = line.split() 

1126 if len(fields) == 0: 

1127 break 

1128 elif 'Fractional coordinates of atoms' in line: 

1129 species = [] 

1130 custom_species = None # A CASTEP special thing 

1131 positions_frac = [] 

1132 # positions_cart = [] 

1133 while True: 

1134 line = out.readline() 

1135 fields = line.split() 

1136 if len(fields) == 7: 

1137 break 

1138 for n in range(n_atoms): 

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

1140 elem = spec_custom[0] 

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

1142 # Add it to the custom info! 

1143 custom_species = list(species) 

1144 species.append(elem) 

1145 if custom_species is not None: 

1146 custom_species.append(fields[1]) 

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

1148 line = out.readline() 

1149 fields = line.split() 

1150 positions_frac_list.append(positions_frac) 

1151 elif 'Files used for pseudopotentials' in line: 

1152 while True: 

1153 line = out.readline() 

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

1155 continue 

1156 fields = line.split() 

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

1158 elem, pp_file = fields 

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

1160 else: 

1161 break 

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

1163 # TODO: generalize for non-Monkhorst Pack case 

1164 # (i.e. kpoint lists) - 

1165 # kpoints_offset cannot be read this way and 

1166 # is hence always set to None 

1167 while True: 

1168 line = out.readline() 

1169 if not line.strip(): 

1170 break 

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

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

1173 # self.kpoints_mp_grid = kpoints 

1174 # self.kpoints_mp_offset = '0. 0. 0.' 

1175 # not set here anymore because otherwise 

1176 # two calculator objects go out of sync 

1177 # after each calculation triggering unnecessary 

1178 # recalculation 

1179 break 

1180 elif 'Symmetry and Constraints' in line: 

1181 # this is a bit of a hack, but otherwise the read_symops 

1182 # would need to re-read the entire file. --> just rewind 

1183 # back by one line, so the read_symops routine can find the 

1184 # start of this block. 

1185 out.seek(_line_start) 

1186 self.read_symops(castep_castep=out) 

1187 elif 'Number of cell constraints' in line: 

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

1189 elif 'Final energy' in line: 

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

1191 elif 'Final free energy' in line: 

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

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

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

1195 # check if we had a finite basis set correction 

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

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

1198 

1199 # Add support for dispersion correction 

1200 # filtering due to SEDC is done in get_potential_energy 

1201 elif 'Dispersion corrected final energy' in line: 

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

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

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

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

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

1207 

1208 # remember to remove constraint labels in force components 

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

1210 # the CASTEP output) 

1211 elif '******************** Forces *********************'\ 

1212 in line or\ 

1213 '************** Symmetrised Forces ***************'\ 

1214 in line or\ 

1215 '************** Constrained Symmetrised Forces ****'\ 

1216 '**********'\ 

1217 in line or\ 

1218 '******************** Constrained Forces **********'\ 

1219 '**********'\ 

1220 in line or\ 

1221 '******************* Unconstrained Forces *********'\ 

1222 '**********'\ 

1223 in line: 

1224 fix = [] 

1225 fix_cart = [] 

1226 forces = [] 

1227 while True: 

1228 line = out.readline() 

1229 fields = line.split() 

1230 if len(fields) == 7: 

1231 break 

1232 for n in range(n_atoms): 

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

1234 fxyz = [0, 0, 0] 

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

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

1237 consd[i] = 1 

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

1239 "(cons'd)", '')) 

1240 if consd.all(): 

1241 fix.append(n) 

1242 elif consd.any(): 

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

1244 forces.append(fxyz) 

1245 line = out.readline() 

1246 fields = line.split() 

1247 

1248 # add support for Hirshfeld analysis 

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

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

1251 # Hirshfeld flag (if castep_keywords.py matches employed 

1252 # castep version) 

1253 calculate_hirshfeld = True 

1254 hirsh_volrat = [] 

1255 while True: 

1256 line = out.readline() 

1257 fields = line.split() 

1258 if len(fields) == 1: 

1259 break 

1260 for n in range(n_atoms): 

1261 hirsh_atom = float(fields[0]) 

1262 hirsh_volrat.append(hirsh_atom) 

1263 while True: 

1264 line = out.readline() 

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

1266 'Hirshfeld Analysis' in line: 

1267 break 

1268 line = out.readline() 

1269 fields = line.split() 

1270 

1271 elif '***************** Stress Tensor *****************'\ 

1272 in line or\ 

1273 '*********** Symmetrised Stress Tensor ***********'\ 

1274 in line: 

1275 stress = [] 

1276 while True: 

1277 line = out.readline() 

1278 fields = line.split() 

1279 if len(fields) == 6: 

1280 break 

1281 for n in range(3): 

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

1283 line = out.readline() 

1284 fields = line.split() 

1285 line = out.readline() 

1286 if "Pressure:" in line: 

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

1288 elif ('BFGS: starting iteration' in line 

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

1290 if n_cell_const < 6: 

1291 lattice_real = [] 

1292 lattice_reci = [] 

1293 # backup previous configuration first: 

1294 # for highly symmetric systems (where essentially only the 

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

1296 # are only printed once. 

1297 if species: 

1298 prev_species = deepcopy(species) 

1299 if positions_frac: 

1300 prev_positions_frac = deepcopy(positions_frac) 

1301 species = [] 

1302 positions_frac = [] 

1303 forces = [] 

1304 

1305 # HOTFIX: 

1306 # Same reason for the stress initialization as before 

1307 # stress = [] 

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

1309 

1310 # extract info from the Mulliken analysis 

1311 elif 'Atomic Populations' in line: 

1312 # sometimes this appears twice in a castep file 

1313 

1314 mulliken_analysis = True 

1315 # skip the separating line 

1316 line = out.readline() 

1317 # this is the headline 

1318 line = out.readline() 

1319 

1320 if 'Charge' in line: 

1321 # skip the next separator line 

1322 line = out.readline() 

1323 while True: 

1324 line = out.readline() 

1325 fields = line.split() 

1326 if len(fields) == 1: 

1327 break 

1328 

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

1330 # outformat changes 

1331 if spin_polarized: 

1332 if len(fields) != 7: 

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

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

1335 else: 

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

1337 

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

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

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

1341 # break 

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

1343 self._warnings.append(line) 

1344 

1345 # fetch some last info 

1346 elif 'Total time' in line: 

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

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

1349 

1350 elif 'Peak Memory Use' in line: 

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

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

1353 

1354 except Exception as exception: 

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

1356 + str(exception)) 

1357 raise 

1358 

1359 if _close: 

1360 out.close() 

1361 

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

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

1364 if not positions_frac: 

1365 positions_frac = prev_positions_frac 

1366 if not species: 

1367 species = prev_species 

1368 

1369 if not spin_polarized: 

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

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

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

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

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

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

1376 

1377 positions_frac_atoms = np.array(positions_frac) 

1378 forces_atoms = np.array(forces) 

1379 spins_atoms = np.array(spins) 

1380 

1381 if mulliken_analysis: 

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

1383 warnings.warn( 

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

1385 ' charges will be ignored') 

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

1387 mulliken_charges_atoms = np.array(mulliken_charges) 

1388 else: 

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

1390 

1391 if hirshfeld_analysis: 

1392 hirshfeld_charges_atoms = np.array(hirshfeld_charges) 

1393 else: 

1394 hirshfeld_charges_atoms = None 

1395 

1396 if calculate_hirshfeld: 

1397 hirsh_atoms = np.array(hirsh_volrat) 

1398 else: 

1399 hirsh_atoms = np.zeros_like(spins) 

1400 

1401 if self.atoms and not self._set_atoms: 

1402 # compensate for internal reordering of atoms by CASTEP 

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

1404 

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

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

1407 

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

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

1410 

1411 # species_castep = list(species) 

1412 forces_castep = np.array(forces) 

1413 hirsh_castep = np.array(hirsh_volrat) 

1414 spins_castep = np.array(spins) 

1415 mulliken_charges_castep = np.array(mulliken_charges_atoms) 

1416 

1417 # go through the atoms position list and replace 

1418 # with the corresponding one from the 

1419 # castep file corresponding atomic number 

1420 for iase in range(n_atoms): 

1421 for icastep in range(n_atoms): 

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

1423 and not atoms_assigned[icastep]): 

1424 positions_frac_atoms[iase] = \ 

1425 positions_frac_castep[icastep] 

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

1427 if iprint > 1 and calculate_hirshfeld: 

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

1429 if spin_polarized: 

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

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

1432 mulliken_charges_atoms[iase] = np.array( 

1433 mulliken_charges_castep[icastep]) 

1434 atoms_assigned[icastep] = True 

1435 break 

1436 

1437 if not all(atoms_assigned): 

1438 not_assigned = [i for (i, assigned) 

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

1440 atoms_assigned) if not assigned] 

1441 warnings.warn( 

1442 '%s atoms not assigned. ' 

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

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

1445 else: 

1446 self.atoms.set_scaled_positions(positions_frac_atoms) 

1447 

1448 else: 

1449 # If no atoms, object has been previously defined 

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

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

1452 

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

1454 # set_calculator also set atoms in the calculator. 

1455 if self.atoms: 

1456 constraints = self.atoms.constraints 

1457 else: 

1458 constraints = [] 

1459 atoms = ase.atoms.Atoms(species, 

1460 cell=lattice_real, 

1461 constraint=constraints, 

1462 pbc=True, 

1463 scaled_positions=positions_frac, 

1464 ) 

1465 if custom_species is not None: 

1466 atoms.new_array('castep_custom_species', 

1467 np.array(custom_species)) 

1468 

1469 if self.param.spin_polarized: 

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

1471 # calculation 

1472 # this one fails as is 

1473 atoms.set_initial_magnetic_moments(magmoms=spins_atoms) 

1474 

1475 if mulliken_analysis: 

1476 atoms.set_initial_charges(charges=mulliken_charges_atoms) 

1477 atoms.calc = self 

1478 

1479 self._kpoints = kpoints 

1480 self._forces = forces_atoms 

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

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

1483 self._hirsh_volrat = hirsh_atoms 

1484 self._spins = spins_atoms 

1485 self._mulliken_charges = mulliken_charges_atoms 

1486 self._hirshfeld_charges = hirshfeld_charges_atoms 

1487 

1488 if self._warnings: 

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

1490 for warning in self._warnings: 

1491 warnings.warn(warning) 

1492 # reset 

1493 self._warnings = [] 

1494 

1495 # Read in eigenvalues from bands file 

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

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

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

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

1500 else: 

1501 try: 

1502 (self._ibz_kpts, 

1503 self._ibz_weights, 

1504 self._eigenvalues, 

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

1506 except FileNotFoundError: 

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

1508 'Fermi energy are unknown') 

1509 

1510 def read_symops(self, castep_castep=None): 

1511 # TODO: check that this is really backwards compatible 

1512 # with previous routine with this name... 

1513 """Read all symmetry operations used from a .castep file.""" 

1514 

1515 if castep_castep is None: 

1516 castep_castep = self._seed + '.castep' 

1517 

1518 if isinstance(castep_castep, str): 

1519 if not os.path.isfile(castep_castep): 

1520 warnings.warn('Warning: CASTEP file %s not found!' % 

1521 castep_castep) 

1522 f = paropen(castep_castep, 'r') 

1523 _close = True 

1524 else: 

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

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

1527 f = castep_castep 

1528 

1529 # look before you leap... 

1530 attributes = ['name', 

1531 'readline', 

1532 'close'] 

1533 

1534 for attr in attributes: 

1535 if not hasattr(f, attr): 

1536 raise TypeError('read_castep_castep_symops: castep_castep ' 

1537 'is not of type str nor valid fileobj!') 

1538 

1539 castep_castep = f.name 

1540 _close = False 

1541 

1542 while True: 

1543 line = f.readline() 

1544 if not line: 

1545 return 

1546 if 'output verbosity' in line: 

1547 iprint = line.split()[-1][1] 

1548 # filter out the default 

1549 if int(iprint) != 1: 

1550 self.param.iprint = iprint 

1551 if 'Symmetry and Constraints' in line: 

1552 break 

1553 

1554 if self.param.iprint.value is None or int(self.param.iprint.value) < 2: 

1555 self._interface_warnings.append( 

1556 'Warning: No symmetry' 

1557 'operations could be read from %s (iprint < 2).' % f.name) 

1558 return 

1559 

1560 while True: 

1561 line = f.readline() 

1562 if not line: 

1563 break 

1564 if 'Number of symmetry operations' in line: 

1565 nsym = int(line.split()[5]) 

1566 # print "nsym = %d" % nsym 

1567 # information about symmetry related atoms currently not read 

1568 symmetry_operations = [] 

1569 for _ in range(nsym): 

1570 rotation = [] 

1571 displacement = [] 

1572 while True: 

1573 if 'rotation' in f.readline(): 

1574 break 

1575 for _ in range(3): 

1576 line = f.readline() 

1577 rotation.append([float(r) for r in line.split()[1:4]]) 

1578 while True: 

1579 if 'displacement' in f.readline(): 

1580 break 

1581 line = f.readline() 

1582 displacement = [float(d) for d in line.split()[1:4]] 

1583 symop = {'rotation': rotation, 

1584 'displacement': displacement} 

1585 self.symmetry_ops = symop 

1586 self.symmetry = symmetry_operations 

1587 warnings.warn( 

1588 'Symmetry operations successfully read from %s. %s' % 

1589 (f.name, self.cell.symmetry_ops)) 

1590 break 

1591 

1592 # only close if we opened the file in this routine 

1593 if _close: 

1594 f.close() 

1595 

1596 def get_hirsh_volrat(self): 

1597 """ 

1598 Return the Hirshfeld volumes. 

1599 """ 

1600 return self._hirsh_volrat 

1601 

1602 def get_spins(self): 

1603 """ 

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

1605 """ 

1606 return self._spins 

1607 

1608 def get_mulliken_charges(self): 

1609 """ 

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

1611 """ 

1612 return self._mulliken_charges 

1613 

1614 def get_hirshfeld_charges(self): 

1615 """ 

1616 Return the charges from a Hirshfeld analysis. 

1617 """ 

1618 return self._hirshfeld_charges 

1619 

1620 def get_total_time(self): 

1621 """ 

1622 Return the total runtime 

1623 """ 

1624 return self._total_time 

1625 

1626 def get_peak_memory(self): 

1627 """ 

1628 Return the peak memory usage 

1629 """ 

1630 return self._peak_memory 

1631 

1632 def set_label(self, label): 

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

1634 in each CASTEP related file. 

1635 """ 

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

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

1638 self._label = label 

1639 

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

1641 notelems=None, 

1642 clear=True, 

1643 suffix='usp'): 

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

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

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

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

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

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

1650 

1651 Parameters :: 

1652 

1653 - elems (None) : set only these elements 

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

1655 - clear (True): clear previous settings 

1656 - suffix (usp): PP file suffix 

1657 """ 

1658 if self._find_pspots: 

1659 if self._pedantic: 

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

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

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

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

1664 'the same <pspot>.') 

1665 

1666 if clear and not elems and not notelems: 

1667 self.cell.species_pot.clear() 

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

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

1670 continue 

1671 if notelems is not None and elem in notelems: 

1672 continue 

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

1674 

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

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

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

1678 castep_pp_path: 

1679 

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

1681 are actually available from the castep_pp_path. 

1682 

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

1684 does a regex matching. The respective pattern is: 

1685 

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

1687 

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

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

1690 <castep_pp_path>. 

1691 

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

1693 (multiple files per element). 

1694 

1695 Parameters :: 

1696 

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

1698 specified. 

1699 - elems (None) : set only these elements 

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

1701 - clear (True): clear previous settings 

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

1703 """ 

1704 if clear and not elems and not notelems: 

1705 self.cell.species_pot.clear() 

1706 

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

1708 if self._pedantic: 

1709 warnings.warn( 

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

1711 .format(self._castep_pp_path)) 

1712 return 

1713 

1714 # translate the bash wildcard syntax to regex 

1715 if pspot == '*': 

1716 pspot = '.*' 

1717 if suffix == '*': 

1718 suffix = '.*' 

1719 if pspot == '*': 

1720 pspot = '.*' 

1721 

1722 # GBRV USPPs have a strnage naming schme 

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

1724 

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

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

1727 continue 

1728 if notelems is not None and elem in notelems: 

1729 continue 

1730 p = pattern.format(elem=elem, 

1731 elem_upper=elem.upper(), 

1732 elem_lower=elem.lower(), 

1733 pspot=pspot, 

1734 suffix=suffix) 

1735 pps = [] 

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

1737 if re.match(p, f): 

1738 pps.append(f) 

1739 if not pps: 

1740 if self._pedantic: 

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

1742 .format(elem)) 

1743 elif not len(pps) == 1: 

1744 raise RuntimeError( 

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

1746 .format(elem) 

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

1748 .format(self._castep_pp_path) 

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

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

1751 else: 

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

1753 

1754 @property 

1755 def name(self): 

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

1757 return self.__name__ 

1758 

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

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

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

1762 # based this class on `Calculator` 

1763 if name == 'forces': 

1764 return self.get_forces(atoms) 

1765 elif name == 'energy': 

1766 return self.get_potential_energy(atoms) 

1767 elif name == 'stress': 

1768 return self.get_stress(atoms) 

1769 elif name == 'charges': 

1770 return self.get_charges(atoms) 

1771 else: 

1772 raise PropertyNotImplementedError 

1773 

1774 @_self_getter 

1775 def get_forces(self, atoms): 

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

1777 self.update(atoms) 

1778 return np.array(self._forces) 

1779 

1780 @_self_getter 

1781 def get_total_energy(self, atoms): 

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

1783 self.update(atoms) 

1784 return self._energy_total 

1785 

1786 @_self_getter 

1787 def get_total_energy_corrected(self, atoms): 

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

1789 self.update(atoms) 

1790 return self._energy_total_corr 

1791 

1792 @_self_getter 

1793 def get_free_energy(self, atoms): 

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

1795 Only defined with smearing.""" 

1796 self.update(atoms) 

1797 return self._energy_free 

1798 

1799 @_self_getter 

1800 def get_0K_energy(self, atoms): 

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

1802 Only defined with smearing.""" 

1803 self.update(atoms) 

1804 return self._energy_0K 

1805 

1806 @_self_getter 

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

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

1809 # but accessing only _name variables 

1810 """Return the total potential energy.""" 

1811 self.update(atoms) 

1812 if force_consistent: 

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

1814 # respective value will default to None as initialized. 

1815 if self._dispcorr_energy_free is not None: 

1816 return self._dispcorr_energy_free 

1817 else: 

1818 return self._energy_free 

1819 else: 

1820 if self._energy_0K is not None: 

1821 if self._dispcorr_energy_0K is not None: 

1822 return self._dispcorr_energy_0K 

1823 else: 

1824 return self._energy_0K 

1825 else: 

1826 if self._dispcorr_energy_total is not None: 

1827 return self._dispcorr_energy_total 

1828 else: 

1829 if self._energy_total_corr is not None: 

1830 return self._energy_total_corr 

1831 else: 

1832 return self._energy_total 

1833 

1834 @_self_getter 

1835 def get_stress(self, atoms): 

1836 """Return the stress.""" 

1837 self.update(atoms) 

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

1839 # annoying user warnings 

1840 stress = np.array( 

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

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

1843 # return self._stress 

1844 return stress 

1845 

1846 @_self_getter 

1847 def get_pressure(self, atoms): 

1848 """Return the pressure.""" 

1849 self.update(atoms) 

1850 return self._pressure 

1851 

1852 @_self_getter 

1853 def get_unit_cell(self, atoms): 

1854 """Return the unit cell.""" 

1855 self.update(atoms) 

1856 return self._unit_cell 

1857 

1858 @_self_getter 

1859 def get_kpoints(self, atoms): 

1860 """Return the kpoints.""" 

1861 self.update(atoms) 

1862 return self._kpoints 

1863 

1864 @_self_getter 

1865 def get_number_cell_constraints(self, atoms): 

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

1867 self.update(atoms) 

1868 return self._number_of_cell_constraints 

1869 

1870 @_self_getter 

1871 def get_charges(self, atoms): 

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

1873 self.update(atoms) 

1874 return np.array(self._mulliken_charges) 

1875 

1876 @_self_getter 

1877 def get_magnetic_moments(self, atoms): 

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

1879 self.update(atoms) 

1880 return np.array(self._spins) 

1881 

1882 def set_atoms(self, atoms): 

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

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

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

1886 self.atoms._calc = self 

1887 

1888 def update(self, atoms): 

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

1890 runs calculation if so. 

1891 """ 

1892 if self.calculation_required(atoms): 

1893 self.calculate(atoms) 

1894 

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

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

1897 settings since the last calculation using this instance. 

1898 """ 

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

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

1901 # self.atoms 

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

1903 if not atoms == self._old_atoms: 

1904 return True 

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

1906 return True 

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

1908 return True 

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

1910 return True 

1911 return False 

1912 

1913 def calculate(self, atoms): 

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

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

1916 if not self._prepare_input_only: 

1917 self.run() 

1918 self.read() 

1919 

1920 # we need to push the old state here! 

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

1922 # again. 

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

1924 self.push_oldstate() 

1925 

1926 def push_oldstate(self): 

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

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

1929 calculation_required will return False and enquiry functions just 

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

1931 """ 

1932 # make a snapshot of all current input 

1933 # to be able to test if recalculation 

1934 # is necessary 

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

1936 self._old_param = deepcopy(self.param) 

1937 self._old_cell = deepcopy(self.cell) 

1938 

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

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

1941 function names in ASE. 

1942 """ 

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

1944 

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

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

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

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

1949 upload the file manually and read out the results using 

1950 Castep().read(). 

1951 """ 

1952 

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

1954 if self._pedantic: 

1955 warnings.warn( 

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

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

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

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

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

1961 if self._pedantic: 

1962 warnings.warn( 

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

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

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

1966 from ase.io.castep import write_param 

1967 

1968 if atoms is None: 

1969 atoms = self.atoms 

1970 else: 

1971 self.atoms = atoms 

1972 

1973 if force_write is None: 

1974 force_write = self._force_write 

1975 

1976 # if we have new instance of the calculator, 

1977 # move existing results out of the way, first 

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

1979 and self._calls == 0 

1980 and self._rename_existing_dir): 

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

1982 os.rmdir(self._directory) 

1983 else: 

1984 # rename appending creation date of the directory 

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

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

1987 (self._directory, 

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

1989 

1990 # create work directory 

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

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

1993 

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

1995 # if self._calls == 0: 

1996 self._fetch_pspots() 

1997 

1998 # if _try_reuse is requested and this 

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

2000 # the .check file from the previous run 

2001 # this is only necessary if _track_output 

2002 # is set to true 

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

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

2005 self.param.reuse = self._check_file 

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

2007 self.param.reuse = self._castep_bin_file 

2008 self._seed = self._build_castep_seed() 

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

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

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

2012 

2013 # write out the input file 

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

2015 self.atoms, castep_cell=self.cell, 

2016 force_write=force_write) 

2017 

2018 if self._export_settings: 

2019 interface_options = self._opt 

2020 else: 

2021 interface_options = None 

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

2023 check_checkfile=self._check_checkfile, 

2024 force_write=force_write, 

2025 interface_options=interface_options,) 

2026 

2027 def _build_castep_seed(self): 

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

2029 with and without _tracking_output. 

2030 """ 

2031 if self._track_output: 

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

2033 else: 

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

2035 

2036 def _abs_path(self, path): 

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

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

2039 

2040 def run(self): 

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

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

2043 """ 

2044 # change to target directory 

2045 self._calls += 1 

2046 

2047 # run castep itself 

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

2049 self._seed), 

2050 cwd=self._directory) 

2051 if stdout: 

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

2053 if stderr: 

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

2055 

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

2057 # self.push_oldstate() 

2058 

2059 # check for non-empty error files 

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

2061 if os.path.exists(err_file): 

2062 err_file = open(err_file) 

2063 self._error = err_file.read() 

2064 err_file.close() 

2065 if self._error: 

2066 raise RuntimeError(self._error) 

2067 

2068 def __repr__(self): 

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

2070 CASTEP settings along with atoms object. 

2071 """ 

2072 expr = '' 

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

2074 if self.atoms is not None: 

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

2076 else: 

2077 expr += 'None\n' 

2078 

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

2080 expr += str(self.param) 

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

2082 expr += str(self.cell) 

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

2084 for key in self.internal_keys: 

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

2086 

2087 return expr 

2088 

2089 def __getattr__(self, attr): 

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

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

2092 they can be read in again in subsequent calls. 

2093 """ 

2094 if attr in self.internal_keys: 

2095 return self._opt[attr] 

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

2097 raise AttributeError 

2098 elif attr not in self.__dict__: 

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

2100 else: 

2101 return self.__dict__[attr] 

2102 

2103 def __setattr__(self, attr, value): 

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

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

2106 Value assigment is case insensitive! 

2107 """ 

2108 

2109 if attr.startswith('_'): 

2110 # internal variables all start with _ 

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

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

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

2114 cutoff=0.9) 

2115 if attr not in self.internal_keys and similars: 

2116 warnings.warn( 

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

2118 (similars, attr)) 

2119 if attr in self.internal_keys: 

2120 self._opt[attr] = value 

2121 if attr == '_track_output': 

2122 if value: 

2123 self._try_reuse = True 

2124 if self._pedantic: 

2125 warnings.warn( 

2126 'You switched _track_output on. This will ' 

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

2128 'also switched _try_reuse on, which will ' 

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

2130 '_try_reuse = False, if you need ' 

2131 'really separate calculations') 

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

2133 self._try_reuse = False 

2134 if self._pedantic: 

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

2136 else: 

2137 self.__dict__[attr] = value 

2138 return 

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

2140 if value is not None: 

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

2142 raise TypeError( 

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

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

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

2146 value) 

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

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

2149 value) 

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

2151 self.__dict__[attr] = value 

2152 return 

2153 elif attr in self.atoms_obj_keys: 

2154 # keywords which clearly belong to the atoms object are 

2155 # rerouted to go there 

2156 self.atoms.__dict__[attr] = value 

2157 return 

2158 elif attr in self.atoms_keys: 

2159 # CASTEP keywords that should go into the atoms object 

2160 # itself are blocked 

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

2162 'through the atoms object' % attr) 

2163 return 

2164 

2165 attr = attr.lower() 

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

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

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

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

2170 # correction 

2171 if self._kw_tol == 0: 

2172 similars = difflib.get_close_matches( 

2173 attr, 

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

2175 if similars: 

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

2177 (attr, similars[0])) 

2178 else: 

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

2180 else: 

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

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

2183 'objects' % attr) 

2184 return 

2185 

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

2187 # so we first determine which one 

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

2189 comp = 'param' 

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

2191 comp = 'cell' 

2192 else: 

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

2194 'the keyword to an input file') 

2195 

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

2197 

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

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

2200 if isinstance(param, CastepParam): 

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

2202 if option.value is not None: 

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

2204 return 

2205 

2206 elif isinstance(param, str): 

2207 param_file = open(param, 'r') 

2208 _close = True 

2209 

2210 else: 

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

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

2213 param_file = param 

2214 

2215 # look before you leap... 

2216 attributes = ['name', 

2217 'close' 

2218 'readlines'] 

2219 

2220 for attr in attributes: 

2221 if not hasattr(param_file, attr): 

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

2223 'nor valid fileobj') 

2224 

2225 param = param_file.name 

2226 _close = False 

2227 

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

2229 get_interface_options=True) 

2230 

2231 # Add the interface options 

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

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

2234 if val in _tf_table: 

2235 val = _tf_table[val] 

2236 self._opt[k] = val 

2237 

2238 if _close: 

2239 param_file.close() 

2240 

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

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

2243 in a temporary and check wether all variables are initialized 

2244 correctly. This is recommended for every bigger simulation. 

2245 """ 

2246 from ase.io.castep import write_param 

2247 

2248 temp_dir = tempfile.mkdtemp() 

2249 self._fetch_pspots(temp_dir) 

2250 seed = 'dryrun' 

2251 

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

2253 self.atoms, castep_cell=self.cell) 

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

2255 # interface 

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

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

2258 return 

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

2260 

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

2262 seed, 

2263 dryrun_flag)), 

2264 cwd=temp_dir) 

2265 

2266 if stdout: 

2267 print(stdout) 

2268 if stderr: 

2269 print(stderr) 

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

2271 

2272 txt = result_file.read() 

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

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

2275 

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

2277 if m: 

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

2279 else: 

2280 warnings.warn( 

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

2282 

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

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

2285 err_file = open(err_file) 

2286 self._error = err_file.read() 

2287 err_file.close() 

2288 

2289 result_file.close() 

2290 shutil.rmtree(temp_dir) 

2291 

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

2293 return match is not None 

2294 

2295 def _fetch_pspots(self, directory=None): 

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

2297 """ 

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

2299 # honored. 

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

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

2302 # By default CASTEP consults the environment variable 

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

2304 # directories it will check those directories for pseudo- 

2305 # potential files if not in the current directory. 

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

2307 # If however PSPOT_DIR was been accidentally set 

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

2309 # setting CASTEP_PP_PATH to an explicit value will 

2310 # still be honored. 

2311 return 

2312 

2313 if directory is None: 

2314 directory = self._directory 

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

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

2317 pspots = {} 

2318 if self._find_pspots: 

2319 self.find_pspots() 

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

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

2322 line = line.split() 

2323 if line: 

2324 pspots[line[0]] = line[1] 

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

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

2327 if self._build_missing_pspots: 

2328 if self._pedantic: 

2329 warnings.warn( 

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

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

2332 'potentials. ' 

2333 'For sake of numerical consistency and efficiency ' 

2334 'this is discouraged.' % species) 

2335 else: 

2336 raise RuntimeError( 

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

2338 species) 

2339 if self.cell.species_pot.value: 

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

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

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

2343 if (os.path.exists(orig_pspot_file) 

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

2345 if self._copy_pspots: 

2346 shutil.copy(orig_pspot_file, directory) 

2347 elif self._link_pspots: 

2348 os.symlink(orig_pspot_file, cp_pspot_file) 

2349 else: 

2350 if self._pedantic: 

2351 warnings.warn(ppwarning) 

2352 

2353 

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

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

2356 'sure to set the evironment variable PSPOT_DIR ' 

2357 'accordingly!') 

2358 

2359 

2360def get_castep_version(castep_command): 

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

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

2363 has been added; this will be attempted first. 

2364 """ 

2365 import tempfile 

2366 with tempfile.TemporaryDirectory() as temp_dir: 

2367 return _get_castep_version(castep_command, temp_dir) 

2368 

2369 

2370def _get_castep_version(castep_command, temp_dir): 

2371 jname = 'dummy_jobname' 

2372 stdout, stderr = '', '' 

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

2374 try: 

2375 stdout, stderr = subprocess.Popen( 

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

2377 stderr=subprocess.PIPE, 

2378 stdout=subprocess.PIPE, cwd=temp_dir, 

2379 universal_newlines=True).communicate() 

2380 if 'CASTEP version' not in stdout: 

2381 stdout, stderr = subprocess.Popen( 

2382 castep_command.split() + [jname], 

2383 stderr=subprocess.PIPE, 

2384 stdout=subprocess.PIPE, cwd=temp_dir, 

2385 universal_newlines=True).communicate() 

2386 except Exception: # XXX Which kind of exception? 

2387 msg = '' 

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

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

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

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

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

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

2394 msg += stdout 

2395 msg += stderr 

2396 raise CastepVersionError(msg) 

2397 if 'CASTEP version' in stdout: 

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

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

2400 else: 

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

2402 output_txt = output.readlines() 

2403 output.close() 

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

2405 # shutil.rmtree(temp_dir) 

2406 for line in output_txt: 

2407 if 'CASTEP version' in line: 

2408 try: 

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

2410 except ValueError: 

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

2412 return fallback_version 

2413 

2414 

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

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

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

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

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

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

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

2422 """ 

2423 # Takes a while ... 

2424 # Fetch all allowed parameters 

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

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

2427 

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

2429 

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

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

2432 'You can overwrite it by calling ' 

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

2434 return False 

2435 

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

2437 # which will cause problems on future runs 

2438 

2439 castep_version = get_castep_version(castep_command) 

2440 

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

2442 

2443 # Filter out proper keywords 

2444 try: 

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

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

2447 # which also include numbers. 

2448 if castep_version < 7.0: 

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

2450 else: 

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

2452 

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

2454 except Exception: 

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

2456 raise 

2457 

2458 types = set() 

2459 levels = set() 

2460 

2461 processed_n = 0 

2462 to_process = len(raw_options[:fetch_only]) 

2463 

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

2465 

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

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

2468 

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

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

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

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

2473 

2474 processed_n += 1 

2475 

2476 if match is not None: 

2477 match = match.groupdict() 

2478 

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

2480 # with keyword assignment during extraction process from CASTEP 

2481 suffix = None 

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

2483 suffix = 'cell' 

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

2485 suffix = 'param' 

2486 if suffix is None: 

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

2488 ' CELL or PARAMETERS keywords' % option) 

2489 

2490 option = option.lower() 

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

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

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

2494 

2495 if mtyp is None: 

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

2497 continue 

2498 if mlvl is None: 

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

2500 continue 

2501 if mdoc is None: 

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

2503 continue 

2504 

2505 types = types.union([mtyp]) 

2506 levels = levels.union([mlvl]) 

2507 

2508 processed_options[suffix][option] = { 

2509 'keyword': option, 

2510 'option_type': mtyp, 

2511 'level': mlvl, 

2512 'docstring': mdoc 

2513 } 

2514 

2515 processed_n += 1 

2516 

2517 frac = (o_i + 1.0) / to_process 

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

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

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

2521 100 * frac)) 

2522 sys.stdout.flush() 

2523 

2524 else: 

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

2526 % option) 

2527 

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

2529 sys.stdout.flush() 

2530 

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

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

2533 processed_options['castep_version'] = castep_version 

2534 

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

2536 

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

2538 (castep_version, processed_n)) 

2539 return True 

2540 

2541 

2542class CastepOption: 

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

2544 type.""" 

2545 

2546 default_convert_types = { 

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

2548 'defined': 'bool', 

2549 'string': 'str', 

2550 'integer': 'int', 

2551 'real': 'float', 

2552 'integer vector': 'int_vector', 

2553 'real vector': 'float_vector', 

2554 'physical': 'float_physical', 

2555 'block': 'block' 

2556 } 

2557 

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

2559 docstring='No information available'): 

2560 self.keyword = keyword 

2561 self.level = level 

2562 self.type = option_type 

2563 self._value = value 

2564 self.__doc__ = docstring 

2565 

2566 @property 

2567 def value(self): 

2568 

2569 if self._value is not None: 

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

2571 'physical'): 

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

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

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

2575 else: 

2576 return str(self._value) 

2577 

2578 @property 

2579 def raw_value(self): 

2580 # The value, not converted to a string 

2581 return self._value 

2582 

2583 @value.setter # type: ignore 

2584 def value(self, val): 

2585 

2586 if val is None: 

2587 self.clear() 

2588 return 

2589 

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

2591 typeparse = '_parse_%s' % ctype 

2592 try: 

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

2594 except ValueError: 

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

2596 

2597 def clear(self): 

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

2599 self._value = None 

2600 

2601 @staticmethod 

2602 def _parse_bool(value): 

2603 try: 

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

2605 except (KeyError, ValueError): 

2606 raise ValueError() 

2607 return value 

2608 

2609 @staticmethod 

2610 def _parse_str(value): 

2611 value = str(value) 

2612 return value 

2613 

2614 @staticmethod 

2615 def _parse_int(value): 

2616 value = int(value) 

2617 return value 

2618 

2619 @staticmethod 

2620 def _parse_float(value): 

2621 value = float(value) 

2622 return value 

2623 

2624 @staticmethod 

2625 def _parse_int_vector(value): 

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

2627 if isinstance(value, str): 

2628 if ',' in value: 

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

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

2631 

2632 value = np.array(value) 

2633 

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

2635 raise ValueError() 

2636 

2637 return list(value) 

2638 

2639 @staticmethod 

2640 def _parse_float_vector(value): 

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

2642 if isinstance(value, str): 

2643 if ',' in value: 

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

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

2646 

2647 value = np.array(value) * 1.0 

2648 

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

2650 raise ValueError() 

2651 

2652 return list(value) 

2653 

2654 @staticmethod 

2655 def _parse_float_physical(value): 

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

2657 if isinstance(value, str): 

2658 value = value.split() 

2659 

2660 try: 

2661 l = len(value) 

2662 except TypeError: 

2663 l = 1 

2664 value = [value] 

2665 

2666 if l == 1: 

2667 try: 

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

2669 except (TypeError, ValueError): 

2670 raise ValueError() 

2671 elif l == 2: 

2672 try: 

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

2674 except (TypeError, ValueError, IndexError): 

2675 raise ValueError() 

2676 else: 

2677 raise ValueError() 

2678 

2679 return value 

2680 

2681 @staticmethod 

2682 def _parse_block(value): 

2683 

2684 if isinstance(value, str): 

2685 return value 

2686 elif hasattr(value, '__getitem__'): 

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

2688 else: 

2689 raise ValueError() 

2690 

2691 def __repr__(self): 

2692 if self._value: 

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

2694 ).format(**self.__dict__) 

2695 else: 

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

2697 ).format(**self.__dict__) 

2698 return expr 

2699 

2700 def __eq__(self, other): 

2701 if not isinstance(other, CastepOption): 

2702 return False 

2703 else: 

2704 return self.__dict__ == other.__dict__ 

2705 

2706 

2707class CastepOptionDict: 

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

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

2710 

2711 Replaces the old CastepCellDict and CastepParamDict that were defined in 

2712 the castep_keywords.py file. 

2713 """ 

2714 

2715 def __init__(self, options=None): 

2716 object.__init__(self) 

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

2718 # CastepOptions can be compared directly now 

2719 for kw in options: 

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

2721 self._options[opt.keyword] = opt 

2722 self.__dict__[opt.keyword] = opt 

2723 

2724 

2725class CastepInputFile: 

2726 

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

2728 

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

2730 

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

2732 object.__init__(self) 

2733 

2734 if options_dict is None: 

2735 options_dict = CastepOptionDict({}) 

2736 

2737 self._options = options_dict._options 

2738 self.__dict__.update(self._options) 

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

2740 # 0 = no new attributes allowed 

2741 # 1 = new attributes allowed, warning given 

2742 # 2 = new attributes allowed, silent 

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

2744 

2745 # Compile a dictionary for quick check of conflict sets 

2746 self._conflict_dict = { 

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

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

2749 

2750 def __repr__(self): 

2751 expr = '' 

2752 is_default = True 

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

2754 if option.value is not None: 

2755 is_default = False 

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

2757 if is_default: 

2758 expr = 'Default\n' 

2759 

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

2761 return expr 

2762 

2763 def __setattr__(self, attr, value): 

2764 

2765 # Hidden attributes are treated normally 

2766 if attr.startswith('_'): 

2767 self.__dict__[attr] = value 

2768 return 

2769 

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

2771 

2772 if self._perm > 0: 

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

2774 is_str = isinstance(value, str) 

2775 is_block = False 

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

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

2778 is_block = True 

2779 

2780 if self._perm == 0: 

2781 similars = difflib.get_close_matches(attr, 

2782 self._options.keys()) 

2783 if similars: 

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

2785 % (attr, similars[0])) 

2786 else: 

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

2788 elif self._perm == 1: 

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

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

2791 ('block' if is_block else 

2792 'string'))) 

2793 attr = attr.lower() 

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

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

2796 self._options[attr] = opt 

2797 self.__dict__[attr] = opt 

2798 else: 

2799 attr = attr.lower() 

2800 opt = self._options[attr] 

2801 

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

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

2804 

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

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

2807 

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

2809 if not (value is None): 

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

2811 for c in cset: 

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

2813 warnings.warn( 

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

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

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

2817 self._options[c].value = None 

2818 

2819 if hasattr(self, attrparse): 

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

2821 else: 

2822 self._options[attr].value = value 

2823 

2824 def __getattr__(self, name): 

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

2826 raise AttributeError() 

2827 

2828 if self._perm == 1: 

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

2830 

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

2832 option_type='string', value=None) 

2833 

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

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

2836 

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

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

2839 

2840 if types: 

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

2842 attrdict[key] = (val, self._options[key].type) 

2843 

2844 return attrdict 

2845 

2846 

2847class CastepParam(CastepInputFile): 

2848 """CastepParam abstracts the settings that go into the .param file""" 

2849 

2850 _keyword_conflicts = [{'cut_off_energy', 'basis_precision'}, ] 

2851 

2852 def __init__(self, castep_keywords, keyword_tolerance=1): 

2853 self._castep_version = castep_keywords.castep_version 

2854 CastepInputFile.__init__(self, castep_keywords.CastepParamDict(), 

2855 keyword_tolerance) 

2856 

2857 @property 

2858 def castep_version(self): 

2859 return self._castep_version 

2860 

2861 # .param specific parsers 

2862 def _parse_reuse(self, value): 

2863 if value is None: 

2864 return None # Reset the value 

2865 try: 

2866 if self._options['continuation'].value: 

2867 warnings.warn('Cannot set reuse if continuation is set, and ' 

2868 'vice versa. Set the other to None, if you want ' 

2869 'this setting.') 

2870 return None 

2871 except KeyError: 

2872 pass 

2873 return 'default' if (value is True) else str(value) 

2874 

2875 def _parse_continuation(self, value): 

2876 if value is None: 

2877 return None # Reset the value 

2878 try: 

2879 if self._options['reuse'].value: 

2880 warnings.warn('Cannot set reuse if continuation is set, and ' 

2881 'vice versa. Set the other to None, if you want ' 

2882 'this setting.') 

2883 return None 

2884 except KeyError: 

2885 pass 

2886 return 'default' if (value is True) else str(value) 

2887 

2888 

2889class CastepCell(CastepInputFile): 

2890 

2891 """CastepCell abstracts all setting that go into the .cell file""" 

2892 

2893 _keyword_conflicts = [ 

2894 {'kpoint_mp_grid', 'kpoint_mp_spacing', 'kpoint_list', 

2895 'kpoints_mp_grid', 'kpoints_mp_spacing', 'kpoints_list'}, 

2896 {'bs_kpoint_mp_grid', 

2897 'bs_kpoint_mp_spacing', 

2898 'bs_kpoint_list', 

2899 'bs_kpoint_path', 

2900 'bs_kpoints_mp_grid', 

2901 'bs_kpoints_mp_spacing', 

2902 'bs_kpoints_list', 

2903 'bs_kpoints_path'}, 

2904 {'spectral_kpoint_mp_grid', 

2905 'spectral_kpoint_mp_spacing', 

2906 'spectral_kpoint_list', 

2907 'spectral_kpoint_path', 

2908 'spectral_kpoints_mp_grid', 

2909 'spectral_kpoints_mp_spacing', 

2910 'spectral_kpoints_list', 

2911 'spectral_kpoints_path'}, 

2912 {'phonon_kpoint_mp_grid', 

2913 'phonon_kpoint_mp_spacing', 

2914 'phonon_kpoint_list', 

2915 'phonon_kpoint_path', 

2916 'phonon_kpoints_mp_grid', 

2917 'phonon_kpoints_mp_spacing', 

2918 'phonon_kpoints_list', 

2919 'phonon_kpoints_path'}, 

2920 {'fine_phonon_kpoint_mp_grid', 

2921 'fine_phonon_kpoint_mp_spacing', 

2922 'fine_phonon_kpoint_list', 

2923 'fine_phonon_kpoint_path'}, 

2924 {'magres_kpoint_mp_grid', 

2925 'magres_kpoint_mp_spacing', 

2926 'magres_kpoint_list', 

2927 'magres_kpoint_path'}, 

2928 {'elnes_kpoint_mp_grid', 

2929 'elnes_kpoint_mp_spacing', 

2930 'elnes_kpoint_list', 

2931 'elnes_kpoint_path'}, 

2932 {'optics_kpoint_mp_grid', 

2933 'optics_kpoint_mp_spacing', 

2934 'optics_kpoint_list', 

2935 'optics_kpoint_path'}, 

2936 {'supercell_kpoint_mp_grid', 

2937 'supercell_kpoint_mp_spacing', 

2938 'supercell_kpoint_list', 

2939 'supercell_kpoint_path'}, ] 

2940 

2941 def __init__(self, castep_keywords, keyword_tolerance=1): 

2942 self._castep_version = castep_keywords.castep_version 

2943 CastepInputFile.__init__(self, castep_keywords.CastepCellDict(), 

2944 keyword_tolerance) 

2945 

2946 @property 

2947 def castep_version(self): 

2948 return self._castep_version 

2949 

2950 # .cell specific parsers 

2951 def _parse_species_pot(self, value): 

2952 

2953 # Single tuple 

2954 if isinstance(value, tuple) and len(value) == 2: 

2955 value = [value] 

2956 # List of tuples 

2957 if hasattr(value, '__getitem__'): 

2958 pspots = [tuple(map(str.strip, x)) for x in value] 

2959 if not all(map(lambda x: len(x) == 2, value)): 

2960 warnings.warn( 

2961 'Please specify pseudopotentials in python as ' 

2962 'a tuple or a list of tuples formatted like: ' 

2963 '(species, file), e.g. ("O", "path-to/O_OTFG.usp") ' 

2964 'Anything else will be ignored') 

2965 return None 

2966 

2967 text_block = self._options['species_pot'].value 

2968 

2969 text_block = text_block if text_block else '' 

2970 # Remove any duplicates 

2971 for pp in pspots: 

2972 text_block = re.sub(r'\n?\s*%s\s+.*' % pp[0], '', text_block) 

2973 if pp[1]: 

2974 text_block += '\n%s %s' % pp 

2975 

2976 return text_block 

2977 

2978 def _parse_symmetry_ops(self, value): 

2979 if not isinstance(value, tuple) \ 

2980 or not len(value) == 2 \ 

2981 or not value[0].shape[1:] == (3, 3) \ 

2982 or not value[1].shape[1:] == (3,) \ 

2983 or not value[0].shape[0] == value[1].shape[0]: 

2984 warnings.warn('Invalid symmetry_ops block, skipping') 

2985 return 

2986 # Now on to print... 

2987 text_block = '' 

2988 for op_i, (op_rot, op_tranls) in enumerate(zip(*value)): 

2989 text_block += '\n'.join([' '.join([str(x) for x in row]) 

2990 for row in op_rot]) 

2991 text_block += '\n' 

2992 text_block += ' '.join([str(x) for x in op_tranls]) 

2993 text_block += '\n\n' 

2994 

2995 return text_block 

2996 

2997 def _parse_positions_abs_intermediate(self, value): 

2998 return _parse_tss_block(value) 

2999 

3000 def _parse_positions_abs_product(self, value): 

3001 return _parse_tss_block(value) 

3002 

3003 def _parse_positions_frac_intermediate(self, value): 

3004 return _parse_tss_block(value, True) 

3005 

3006 def _parse_positions_frac_product(self, value): 

3007 return _parse_tss_block(value, True) 

3008 

3009 

3010CastepKeywords = namedtuple('CastepKeywords', 

3011 ['CastepParamDict', 'CastepCellDict', 

3012 'types', 'levels', 'castep_version']) 

3013 

3014# We keep this just for naming consistency with older versions 

3015 

3016 

3017def make_cell_dict(data=None): 

3018 

3019 data = data if data is not None else {} 

3020 

3021 class CastepCellDict(CastepOptionDict): 

3022 def __init__(self): 

3023 CastepOptionDict.__init__(self, data) 

3024 

3025 return CastepCellDict 

3026 

3027 

3028def make_param_dict(data=None): 

3029 

3030 data = data if data is not None else {} 

3031 

3032 class CastepParamDict(CastepOptionDict): 

3033 def __init__(self): 

3034 CastepOptionDict.__init__(self, data) 

3035 

3036 return CastepParamDict 

3037 

3038 

3039class CastepVersionError(Exception): 

3040 """No special behaviour, works to signal when Castep can not be found""" 

3041 pass 

3042 

3043 

3044class ConversionError(Exception): 

3045 

3046 """Print customized error for options that are not converted correctly 

3047 and point out that they are maybe not implemented, yet""" 

3048 

3049 def __init__(self, key_type, attr, value): 

3050 Exception.__init__(self) 

3051 self.key_type = key_type 

3052 self.value = value 

3053 self.attr = attr 

3054 

3055 def __str__(self): 

3056 return 'Could not convert %s = %s to %s\n' \ 

3057 % (self.attr, self.value, self.key_type) \ 

3058 + 'This means you either tried to set a value of the wrong\n'\ 

3059 + 'type or this keyword needs some special care. Please feel\n'\ 

3060 + 'to add it to the corresponding __setattr__ method and send\n'\ 

3061 + 'the patch to %s, so we can all benefit.' % (contact_email) 

3062 

3063 

3064def get_castep_pp_path(castep_pp_path=''): 

3065 """Abstract the quest for a CASTEP PSP directory.""" 

3066 if castep_pp_path: 

3067 return os.path.abspath(os.path.expanduser(castep_pp_path)) 

3068 elif 'PSPOT_DIR' in os.environ: 

3069 return os.environ['PSPOT_DIR'] 

3070 elif 'CASTEP_PP_PATH' in os.environ: 

3071 return os.environ['CASTEP_PP_PATH'] 

3072 else: 

3073 return os.path.abspath('.') 

3074 

3075 

3076def get_castep_command(castep_command=''): 

3077 """Abstract the quest for a castep_command string.""" 

3078 if castep_command: 

3079 return castep_command 

3080 elif 'CASTEP_COMMAND' in os.environ: 

3081 return os.environ['CASTEP_COMMAND'] 

3082 else: 

3083 return 'castep' 

3084 

3085 

3086def shell_stdouterr(raw_command, cwd=None): 

3087 """Abstracts the standard call of the commandline, when 

3088 we are only interested in the stdout and stderr 

3089 """ 

3090 stdout, stderr = subprocess.Popen(raw_command, 

3091 stdout=subprocess.PIPE, 

3092 stderr=subprocess.PIPE, 

3093 universal_newlines=True, 

3094 shell=True, cwd=cwd).communicate() 

3095 return stdout.strip(), stderr.strip() 

3096 

3097 

3098def import_castep_keywords(castep_command='', 

3099 filename='castep_keywords.json', 

3100 path='.'): 

3101 

3102 # Search for castep_keywords.json (or however it's called) in multiple 

3103 # paths 

3104 

3105 searchpaths = [path, os.path.expanduser('~/.ase'), 

3106 os.path.expanduser('~/.config/ase')] 

3107 try: 

3108 kwfile = sum([glob.glob(os.path.join(sp, filename)) 

3109 for sp in searchpaths], [])[0] 

3110 except IndexError: 

3111 warnings.warn("""Generating CASTEP keywords JSON file... hang on. 

3112 The CASTEP keywords JSON file contains abstractions for CASTEP input 

3113 parameters (for both .cell and .param input files), including some 

3114 format checks and descriptions. The latter are extracted from the 

3115 internal online help facility of a CASTEP binary, thus allowing to 

3116 easily keep the calculator synchronized with (different versions of) 

3117 the CASTEP code. Consequently, avoiding licensing issues (CASTEP is 

3118 distributed commercially by accelrys), we consider it wise not to 

3119 provide the file in the first place.""") 

3120 create_castep_keywords(get_castep_command(castep_command), 

3121 filename=filename, path=path) 

3122 warnings.warn( 

3123 'Stored %s in %s. Copy it to your ASE installation under ' 

3124 'ase/calculators for system-wide installation. Using a *nix ' 

3125 'OS this can be a simple as mv %s %s' % 

3126 (filename, os.path.abspath(path), 

3127 os.path.join(os.path.abspath(path), filename), 

3128 os.path.join(os.path.dirname(ase.__file__), 

3129 'calculators'))) 

3130 kwfile = os.path.join(path, filename) 

3131 

3132 # Now create the castep_keywords object proper 

3133 kwdata = json.load(open(kwfile)) 

3134 

3135 # This is a bit awkward, but it's necessary for backwards compatibility 

3136 param_dict = make_param_dict(kwdata['param']) 

3137 cell_dict = make_cell_dict(kwdata['cell']) 

3138 

3139 castep_keywords = CastepKeywords(param_dict, cell_dict, 

3140 kwdata['types'], kwdata['levels'], 

3141 kwdata['castep_version']) 

3142 

3143 return castep_keywords 

3144 

3145 

3146if __name__ == '__main__': 

3147 warnings.warn( 

3148 'When called directly this calculator will fetch all available ' 

3149 'keywords from the binarys help function into a ' 

3150 'castep_keywords.json in the current directory %s ' 

3151 'For system wide usage, it can be copied into an ase installation ' 

3152 'at ASE/calculators. ' 

3153 'This castep_keywords.json usually only needs to be generated once ' 

3154 'for a CASTEP binary/CASTEP version.' % os.getcwd()) 

3155 

3156 import optparse 

3157 parser = optparse.OptionParser() 

3158 parser.add_option( 

3159 '-f', '--force-write', dest='force_write', 

3160 help='Force overwriting existing castep_keywords.json', default=False, 

3161 action='store_true') 

3162 (options, args) = parser.parse_args() 

3163 

3164 if args: 

3165 opt_castep_command = ''.join(args) 

3166 else: 

3167 opt_castep_command = '' 

3168 generated = create_castep_keywords(get_castep_command(opt_castep_command), 

3169 force_write=options.force_write) 

3170 

3171 if generated: 

3172 try: 

3173 with open('castep_keywords.json') as fd: 

3174 json.load(fd) 

3175 except Exception as e: 

3176 warnings.warn( 

3177 '%s Ooops, something went wrong with the CASTEP keywords' % e) 

3178 else: 

3179 warnings.warn('Import works. Looking good!')