Coverage for /builds/ase/ase/ase/calculators/castep.py: 49.25%
731 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3"""This module defines an interface to CASTEP for
4 use by the ASE (Webpage: http://wiki.fysik.dtu.dk/ase)
6Authors:
7 Max Hoffmann, max.hoffmann@ch.tum.de
8 Joerg Meyer, joerg.meyer@ch.tum.de
9 Simon P. Rittmeyer, simon.rittmeyer@tum.de
11Contributors:
12 Juan M. Lorenzi, juan.lorenzi@tum.de
13 Georg S. Michelitsch, georg.michelitsch@tch.tum.de
14 Reinhard J. Maurer, reinhard.maurer@yale.edu
15 Simone Sturniolo, simone.sturniolo@stfc.ac.uk
16"""
18import difflib
19import glob
20import json
21import os
22import re
23import shutil
24import subprocess
25import sys
26import tempfile
27import time
28import warnings
29from collections import namedtuple
30from copy import deepcopy
31from itertools import product
32from pathlib import Path
34import numpy as np
36from ase import Atoms
37from ase.calculators.calculator import (
38 BaseCalculator,
39 compare_atoms,
40 kpts2sizeandoffsets,
41)
42from ase.config import cfg
43from ase.dft.kpoints import BandPath
44from ase.io.castep import read_bands, read_param
45from ase.io.castep.castep_input_file import CastepCell, CastepParam
46from ase.io.castep.castep_reader import read_castep_castep
47from ase.parallel import paropen
49__all__ = [
50 'Castep',
51 'CastepCell',
52 'CastepParam',
53 'create_castep_keywords']
55# A convenient table to avoid the previously used "eval"
56_tf_table = {
57 '': True, # Just the keyword is equivalent to True
58 'True': True,
59 'False': False}
62def _self_getter(getf):
63 # A decorator that makes it so that if no 'atoms' argument is passed to a
64 # getter function, self.atoms is used instead
66 def decor_getf(self, atoms=None, *args, **kwargs):
68 if atoms is None:
69 atoms = self.atoms
71 return getf(self, atoms, *args, **kwargs)
73 return decor_getf
76class Castep(BaseCalculator):
77 r"""
78CASTEP Interface Documentation
81Introduction
82============
84CASTEP_ [1]_ W_ is a software package which uses density functional theory to
85provide a good atomic-level description of all manner of materials and
86molecules. CASTEP can give information about total energies, forces and
87stresses on an atomic system, as well as calculating optimum geometries, band
88structures, optical spectra, phonon spectra and much more. It can also perform
89molecular dynamics simulations.
91The CASTEP calculator interface class offers intuitive access to all CASTEP
92settings and most results. All CASTEP specific settings are accessible via
93attribute access (*i.e*. ``calc.param.keyword = ...`` or
94``calc.cell.keyword = ...``)
97Getting Started:
98================
100Set the environment variables appropriately for your system::
102 export CASTEP_COMMAND=' ... '
103 export CASTEP_PP_PATH=' ... '
105Note: alternatively to CASTEP_PP_PATH one can set PSPOT_DIR
106as CASTEP consults this by default, i.e.::
108 export PSPOT_DIR=' ... '
111Running the Calculator
112======================
114The default initialization command for the CASTEP calculator is
116.. class:: Castep(directory='CASTEP', label='castep')
118To do a minimal run one only needs to set atoms, this will use all
119default settings of CASTEP, meaning LDA, singlepoint, etc..
121With a generated *castep_keywords.json* in place all options are accessible
122by inspection, *i.e.* tab-completion. This works best when using ``ipython``.
123All options can be accessed via ``calc.param.<TAB>`` or ``calc.cell.<TAB>``
124and documentation is printed with ``calc.param.<keyword> ?`` or
125``calc.cell.<keyword> ?``. All options can also be set directly
126using ``calc.keyword = ...`` or ``calc.KEYWORD = ...`` or even
127``ialc.KeYwOrD`` or directly as named arguments in the call to the constructor
128(*e.g.* ``Castep(task='GeometryOptimization')``).
129If using this calculator on a machine without CASTEP, one might choose to copy
130a *castep_keywords.json* file generated elsewhere in order to access this
131feature: the file will be used if located in the working directory,
132*$HOME/.ase/* or *ase/ase/calculators/* within the ASE library. The file should
133be generated the first time it is needed, but you can generate a new keywords
134file in the currect directory with ``python -m ase.calculators.castep``.
136All options that go into the ``.param`` file are held in an ``CastepParam``
137instance, while all options that go into the ``.cell`` file and don't belong
138to the atoms object are held in an ``CastepCell`` instance. Each instance can
139be created individually and can be added to calculators by attribute
140assignment, *i.e.* ``calc.param = param`` or ``calc.cell = cell``.
142All internal variables of the calculator start with an underscore (_).
143All cell attributes that clearly belong into the atoms object are blocked.
144Setting ``calc.atoms_attribute`` (*e.g.* ``= positions``) is sent directly to
145the atoms object.
148Arguments:
149==========
151========================= ====================================================
152Keyword Description
153========================= ====================================================
154``directory`` The relative path where all input and output files
155 will be placed. If this does not exist, it will be
156 created. Existing directories will be moved to
157 directory-TIMESTAMP unless self._rename_existing_dir
158 is set to false.
160``label`` The prefix of .param, .cell, .castep, etc. files.
162``castep_command`` Command to run castep. Can also be set via the bash
163 environment variable ``CASTEP_COMMAND``. If none is
164 given or found, will default to ``castep``
166``check_castep_version`` Boolean whether to check if the installed castep
167 version matches the version from which the available
168 options were deduced. Defaults to ``False``.
170``castep_pp_path`` The path where the pseudopotentials are stored. Can
171 also be set via the bash environment variables
172 ``PSPOT_DIR`` (preferred) and ``CASTEP_PP_PATH``.
173 Will default to the current working directory if
174 none is given or found. Note that pseudopotentials
175 may be generated on-the-fly if they are not found.
177``find_pspots`` Boolean whether to search for pseudopotentials in
178 ``<castep_pp_path>`` or not. If activated, files in
179 this directory will be checked for typical names. If
180 files are not found, they will be generated on the
181 fly, depending on the ``_build_missing_pspots``
182 value. A RuntimeError will be raised in case
183 multiple files per element are found. Defaults to
184 ``False``.
185``keyword_tolerance`` Integer to indicate the level of tolerance to apply
186 validation of any parameters set in the CastepCell
187 or CastepParam objects against the ones found in
188 castep_keywords. Levels are as following:
190 0 = no tolerance, keywords not found in
191 castep_keywords will raise an exception
193 1 = keywords not found will be accepted but produce
194 a warning (default)
196 2 = keywords not found will be accepted silently
198 3 = no attempt is made to look for
199 castep_keywords.json at all
200``castep_keywords`` Can be used to pass a CastepKeywords object that is
201 then used with no attempt to actually load a
202 castep_keywords.json file. Most useful for debugging
203 and testing purposes.
205========================= ====================================================
208Additional Settings
209===================
211========================= ====================================================
212Internal Setting Description
213========================= ====================================================
214``_castep_command`` (``=castep``): the actual shell command used to
215 call CASTEP.
217``_check_checkfile`` (``=True``): this makes write_param() only
218 write a continue or reuse statement if the
219 addressed .check or .castep_bin file exists in the
220 directory.
222``_copy_pspots`` (``=False``): if set to True the calculator will
223 actually copy the needed pseudo-potential (\*.usp)
224 file, usually it will only create symlinks.
226``_link_pspots`` (``=True``): if set to True the calculator will
227 actually will create symlinks to the needed pseudo
228 potentials. Set this option (and ``_copy_pspots``)
229 to False if you rather want to access your pseudo
230 potentials using the PSPOT_DIR environment variable
231 that is read by CASTEP.
232 *Note:* This option has no effect if ``copy_pspots``
233 is True..
235``_build_missing_pspots`` (``=True``): if set to True, castep will generate
236 missing pseudopotentials on the fly. If not, a
237 RuntimeError will be raised if not all files were
238 found.
240``_export_settings`` (``=True``): if this is set to
241 True, all calculator internal settings shown here
242 will be included in the .param in a comment line (#)
243 and can be read again by merge_param. merge_param
244 can be forced to ignore this directive using the
245 optional argument ``ignore_internal_keys=True``.
247``_force_write`` (``=True``): this controls wether the \*cell and
248 \*param will be overwritten.
250``_prepare_input_only`` (``=False``): If set to True, the calculator will
251 create \*cell und \*param file but not
252 start the calculation itself.
253 If this is used to prepare jobs locally
254 and run on a remote cluster it is recommended
255 to set ``_copy_pspots = True``.
257``_castep_pp_path`` (``='.'``) : the place where the calculator
258 will look for pseudo-potential files.
260``_find_pspots`` (``=False``): if set to True, the calculator will
261 try to find the respective pseudopotentials from
262 <_castep_pp_path>. As long as there are no multiple
263 files per element in this directory, the auto-detect
264 feature should be very robust. Raises a RuntimeError
265 if required files are not unique (multiple files per
266 element). Non existing pseudopotentials will be
267 generated, though this could be dangerous.
269``_rename_existing_dir`` (``=True``) : when using a new instance
270 of the calculator, this will move directories out of
271 the way that would be overwritten otherwise,
272 appending a date string.
274``_set_atoms`` (``=False``) : setting this to True will overwrite
275 any atoms object previously attached to the
276 calculator when reading a \.castep file. By de-
277 fault, the read() function will only create a new
278 atoms object if none has been attached and other-
279 wise try to assign forces etc. based on the atom's
280 positions. ``_set_atoms=True`` could be necessary
281 if one uses CASTEP's internal geometry optimization
282 (``calc.param.task='GeometryOptimization'``)
283 because then the positions get out of sync.
284 *Warning*: this option is generally not recommended
285 unless one knows one really needs it. There should
286 never be any need, if CASTEP is used as a
287 single-point calculator.
289``_track_output`` (``=False``) : if set to true, the interface
290 will append a number to the label on all input
291 and output files, where n is the number of calls
292 to this instance. *Warning*: this setting may con-
293 sume a lot more disk space because of the additio-
294 nal \*check files.
296``_try_reuse`` (``=_track_output``) : when setting this, the in-
297 terface will try to fetch the reuse file from the
298 previous run even if _track_output is True. By de-
299 fault it is equal to _track_output, but may be
300 overridden.
302 Since this behavior may not always be desirable for
303 single-point calculations. Regular reuse for *e.g.*
304 a geometry-optimization can be achieved by setting
305 ``calc.param.reuse = True``.
307``_pedantic`` (``=False``) if set to true, the calculator will
308 inform about settings probably wasting a lot of CPU
309 time or causing numerical inconsistencies.
311========================= ====================================================
313Special features:
314=================
317``.dryrun_ok()``
318 Runs ``castep_command seed -dryrun`` in a temporary directory return True if
319 all variables initialized ok. This is a fast way to catch errors in the
320 input. Afterwards _kpoints_used is set.
322``.merge_param()``
323 Takes a filename or filehandler of a .param file or CastepParam instance and
324 merges it into the current calculator instance, overwriting current settings
326``.keyword.clear()``
327 Can be used on any option like ``calc.param.keyword.clear()`` or
328 ``calc.cell.keyword.clear()`` to return to the CASTEP default.
330``.initialize()``
331 Creates all needed input in the ``_directory``. This can then copied to and
332 run in a place without ASE or even python.
334``.set_pspot('<library>')``
335 This automatically sets the pseudo-potential for all present species to
336 ``<Species>_<library>.usp``. Make sure that ``_castep_pp_path`` is set
337 correctly. Note that there is no check, if the file actually exists. If it
338 doesn't castep will crash! You may want to use ``find_pspots()`` instead.
340``.find_pspots(pspot=<library>, suffix=<suffix>)``
341 This automatically searches for pseudopotentials of type
342 ``<Species>_<library>.<suffix>`` or ``<Species>-<library>.<suffix>`` in
343 ``castep_pp_path` (make sure this is set correctly). Note that ``<Species>``
344 will be searched for case insensitive. Regular expressions are accepted, and
345 arguments ``'*'`` will be regarded as bash-like wildcards. Defaults are any
346 ``<library>`` and any ``<suffix>`` from ``['usp', 'UPF', 'recpot']``. If you
347 have well-organized folders with pseudopotentials of one kind, this should
348 work with the defaults.
350``print(calc)``
351 Prints a short summary of the calculator settings and atoms.
353``ase.io.castep.read_seed('path-to/seed')``
354 Given you have a combination of seed.{param,cell,castep} this will return an
355 atoms object with the last ionic positions in the .castep file and all other
356 settings parsed from the .cell and .param file. If no .castep file is found
357 the positions are taken from the .cell file. The output directory will be
358 set to the same directory, only the label is preceded by 'copy_of\_' to
359 avoid overwriting.
361``.set_kpts(kpoints)``
362 This is equivalent to initialising the calculator with
363 ``calc = Castep(kpts=kpoints)``. ``kpoints`` can be specified in many
364 convenient forms: simple Monkhorst-Pack grids can be specified e.g.
365 ``(2, 2, 3)`` or ``'2 2 3'``; lists of specific weighted k-points can be
366 given in reciprocal lattice coordinates e.g.
367 ``[[0, 0, 0, 0.25], [0.25, 0.25, 0.25, 0.75]]``; a dictionary syntax is
368 available for more complex requirements e.g.
369 ``{'size': (2, 2, 2), 'gamma': True}`` will give a Gamma-centered 2x2x2 M-P
370 grid, ``{'density': 10, 'gamma': False, 'even': False}`` will give a mesh
371 with density of at least 10 Ang (based on the unit cell of currently-attached
372 atoms) with an odd number of points in each direction and avoiding the Gamma
373 point.
375``.set_bandpath(bandpath)``
376 This is equivalent to initialialising the calculator with
377 ``calc=Castep(bandpath=bandpath)`` and may be set simultaneously with *kpts*.
378 It allows an electronic band structure path to be set up using ASE BandPath
379 objects. This enables a band structure calculation to be set up conveniently
380 using e.g. calc.set_bandpath(atoms.cell.bandpath().interpolate(npoints=200))
382``.band_structure(bandfile=None)``
383 Read a band structure from _seedname.bands_ file. This returns an ase
384 BandStructure object which may be plotted with e.g.
385 ``calc.band_structure().plot()``
387Notes/Issues:
388==============
390* Currently *only* the FixAtoms *constraint* is fully supported for reading and
391 writing. There is some experimental support for the FixCartesian constraint.
393* There is no support for the CASTEP *unit system*. Units of eV and Angstrom
394 are used throughout. In particular when converting total energies from
395 different calculators, one should check that the same CODATA_ version is
396 used for constants and conversion factors, respectively.
398.. _CASTEP: http://www.castep.org/
400.. _W: https://en.wikipedia.org/wiki/CASTEP
402.. _CODATA: https://physics.nist.gov/cuu/Constants/index.html
404.. [1] S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J. Probert,
405 K. Refson, M. C. Payne Zeitschrift für Kristallographie 220(5-6)
406 pp.567- 570 (2005) PDF_.
408.. _PDF: http://www.tcm.phy.cam.ac.uk/castep/papers/ZKristallogr_2005.pdf
411End CASTEP Interface Documentation
412 """
414 # Class attributes !
415 # keys set through atoms object
416 atoms_keys = [
417 'charges',
418 'ionic_constraints',
419 'lattice_abs',
420 'lattice_cart',
421 'positions_abs',
422 'positions_abs_final',
423 'positions_abs_intermediate',
424 'positions_frac',
425 'positions_frac_final',
426 'positions_frac_intermediate']
428 atoms_obj_keys = [
429 'dipole',
430 'energy_free',
431 'energy_zero',
432 'fermi',
433 'forces',
434 'nbands',
435 'positions',
436 'stress',
437 'pressure']
439 internal_keys = [
440 '_castep_command',
441 '_check_checkfile',
442 '_copy_pspots',
443 '_link_pspots',
444 '_find_pspots',
445 '_build_missing_pspots',
446 '_directory',
447 '_export_settings',
448 '_force_write',
449 '_label',
450 '_prepare_input_only',
451 '_castep_pp_path',
452 '_rename_existing_dir',
453 '_set_atoms',
454 '_track_output',
455 '_try_reuse',
456 '_pedantic']
458 implemented_properties = [
459 'energy',
460 'free_energy',
461 'forces',
462 'stress',
463 'charges',
464 'magmoms',
465 ]
467 # specific to this calculator
468 implemented_properties += [
469 'energy_without_dispersion_correction',
470 'free_energy_without_dispersion_correction',
471 'energy_zero_without_dispersion_correction',
472 'energy_with_dispersion_correction',
473 'free_energy_with_dispersion_correction',
474 'energy_zero_with_dispersion_correction',
475 'energy_with_finite_basis_set_correction',
476 'pressure',
477 'hirshfeld_volume_ratios',
478 'hirshfeld_charges',
479 'hirshfeld_magmoms',
480 ]
482 def __init__(self, directory='CASTEP', label='castep',
483 castep_command=None, check_castep_version=False,
484 castep_pp_path=None, find_pspots=False, keyword_tolerance=1,
485 castep_keywords=None, **kwargs):
487 self.results = {}
489 from ase.io.castep import write_castep_cell
490 self._write_cell = write_castep_cell
492 if castep_keywords is None:
493 castep_keywords = CastepKeywords(make_param_dict(),
494 make_cell_dict(),
495 [],
496 [],
497 0)
498 if keyword_tolerance < 3:
499 try:
500 castep_keywords = import_castep_keywords(castep_command)
501 except CastepVersionError as e:
502 if keyword_tolerance == 0:
503 raise e
504 else:
505 warnings.warn(str(e))
507 self._kw_tol = keyword_tolerance
508 keyword_tolerance = max(keyword_tolerance, 2) # 3 not accepted below
509 self.param = CastepParam(castep_keywords,
510 keyword_tolerance=keyword_tolerance)
511 self.cell = CastepCell(castep_keywords,
512 keyword_tolerance=keyword_tolerance)
514 ###################################
515 # Calculator state variables #
516 ###################################
517 self._calls = 0
518 self._castep_version = castep_keywords.castep_version
520 # collects content from *.err file
521 self._error = None
522 # warnings raised by the ASE interface
523 self._interface_warnings = []
525 # store to check if recalculation is necessary
526 self._old_atoms = None
527 self._old_cell = None
528 self._old_param = None
530 ###################################
531 # Internal keys #
532 # Allow to tweak the behavior #
533 ###################################
534 self._opt = {}
535 self._castep_command = get_castep_command(castep_command)
536 self._castep_pp_path = get_castep_pp_path(castep_pp_path)
537 self._check_checkfile = True
538 self._copy_pspots = False
539 self._link_pspots = True
540 self._find_pspots = find_pspots
541 self._build_missing_pspots = True
542 self._directory = os.path.abspath(directory)
543 self._export_settings = True
544 self._force_write = True
545 self._label = label
546 self._prepare_input_only = False
547 self._rename_existing_dir = True
548 self._set_atoms = False
549 self._track_output = False
550 self._try_reuse = False
552 # turn off the pedantic user warnings
553 self._pedantic = False
555 # will be set on during runtime
556 self._seed = None
558 ###################################
559 # (Physical) result variables #
560 ###################################
561 self.atoms = None
562 # initialize result variables
563 self._eigenvalues = None
564 self._efermi = None
565 self._ibz_kpts = None
566 self._ibz_weights = None
567 self._band_structure = None
569 self._number_of_cell_constraints = None
570 self._output_verbosity = None
571 self._unit_cell = None
572 self._kpoints = None
574 # pointers to other files used at runtime
575 self._check_file = None
576 self._castep_bin_file = None
578 # plane wave cutoff energy (may be derived during PP generation)
579 self._cut_off_energy = None
581 # runtime information
582 self._total_time = None
583 self._peak_memory = None
585 # check version of CASTEP options module against current one
586 if check_castep_version:
587 local_castep_version = get_castep_version(self._castep_command)
588 if not hasattr(self, '_castep_version'):
589 warnings.warn('No castep version found')
590 return
591 if local_castep_version != self._castep_version:
592 warnings.warn(
593 'The options module was generated from version %s '
594 'while your are currently using CASTEP version %s' %
595 (self._castep_version,
596 get_castep_version(self._castep_command)))
597 self._castep_version = local_castep_version
599 # processes optional arguments in kw style
600 for keyword, value in kwargs.items():
601 # first fetch special keywords issued by ASE CLI
602 if keyword == 'kpts':
603 self.set_kpts(value)
604 elif keyword == 'bandpath':
605 self.set_bandpath(value)
606 elif keyword == 'xc':
607 self.xc_functional = value
608 elif keyword == 'ecut':
609 self.cut_off_energy = value
610 else: # the general case
611 self.__setattr__(keyword, value)
613 # TODO: to be self.use_cache = True after revising `__setattr__`
614 self.__dict__['use_cache'] = True
616 def set_atoms(self, atoms):
617 self.atoms = atoms
619 def get_atoms(self):
620 if self.atoms is None:
621 raise ValueError('Calculator has no atoms')
622 atoms = self.atoms.copy()
623 atoms.calc = self
624 return atoms
626 def _get_name(self) -> str:
627 return self.__class__.__name__
629 def band_structure(self, bandfile=None):
630 from ase.spectrum.band_structure import BandStructure
632 if bandfile is None:
633 bandfile = os.path.join(self._directory, self._seed) + '.bands'
635 if not os.path.exists(bandfile):
636 raise ValueError(f'Cannot find band file "{bandfile}".')
638 kpts, _weights, eigenvalues, efermi = read_bands(bandfile)
640 # Get definitions of high-symmetry points
641 special_points = self.atoms.cell.bandpath(npoints=0).special_points
642 bandpath = BandPath(self.atoms.cell,
643 kpts=kpts,
644 special_points=special_points)
645 return BandStructure(bandpath, eigenvalues, reference=efermi)
647 def set_bandpath(self, bandpath):
648 """Set a band structure path from ase.dft.kpoints.BandPath object
650 This will set the bs_kpoint_list block with a set of specific points
651 determined in ASE. bs_kpoint_spacing will not be used; to modify the
652 number of points, consider using e.g. bandpath.resample(density=20) to
653 obtain a new dense path.
655 Args:
656 bandpath (:obj:`ase.dft.kpoints.BandPath` or None):
657 Set to None to remove list of band structure points. Otherwise,
658 sampling will follow BandPath parameters.
660 """
662 def clear_bs_keywords():
663 bs_keywords = product({'bs_kpoint', 'bs_kpoints'},
664 {'path', 'path_spacing',
665 'list',
666 'mp_grid', 'mp_spacing', 'mp_offset'})
667 for bs_tag in bs_keywords:
668 setattr(self.cell, '_'.join(bs_tag), None)
670 if bandpath is None:
671 clear_bs_keywords()
672 elif isinstance(bandpath, BandPath):
673 clear_bs_keywords()
674 self.cell.bs_kpoint_list = [' '.join(map(str, row))
675 for row in bandpath.kpts]
676 else:
677 raise TypeError('Band structure path must be an '
678 'ase.dft.kpoint.BandPath object')
680 def set_kpts(self, kpts):
681 """Set k-point mesh/path using a str, tuple or ASE features
683 Args:
684 kpts (None, tuple, str, dict):
686 This method will set the CASTEP parameters kpoints_mp_grid,
687 kpoints_mp_offset and kpoints_mp_spacing as appropriate. Unused
688 parameters will be set to None (i.e. not included in input files.)
690 If kpts=None, all these parameters are set as unused.
692 The simplest useful case is to give a 3-tuple of integers specifying
693 a Monkhorst-Pack grid. This may also be formatted as a string separated
694 by spaces; this is the format used internally before writing to the
695 input files.
697 A more powerful set of features is available when using a python
698 dictionary with the following allowed keys:
700 - 'size' (3-tuple of int) mesh of mesh dimensions
701 - 'density' (float) for BZ sampling density in points per recip. Ang
702 ( kpoint_mp_spacing = 1 / (2pi * density) ). An explicit MP mesh will
703 be set to allow for rounding/centering.
704 - 'spacing' (float) for BZ sampling density for maximum space between
705 sample points in reciprocal space. This is numerically equivalent to
706 the inbuilt ``calc.cell.kpoint_mp_spacing``, but will be converted to
707 'density' to allow for rounding/centering.
708 - 'even' (bool) to round each direction up to the nearest even number;
709 set False for odd numbers, leave as None for no odd/even rounding.
710 - 'gamma' (bool) to offset the Monkhorst-Pack grid to include
711 (0, 0, 0); set False to offset each direction avoiding 0.
712 """
714 def clear_mp_keywords():
715 mp_keywords = product({'kpoint', 'kpoints'},
716 {'mp_grid', 'mp_offset',
717 'mp_spacing', 'list'})
718 for kp_tag in mp_keywords:
719 setattr(self.cell, '_'.join(kp_tag), None)
721 # Case 1: Clear parameters with set_kpts(None)
722 if kpts is None:
723 clear_mp_keywords()
725 # Case 2: list of explicit k-points with weights
726 # e.g. [[ 0, 0, 0, 0.125],
727 # [ 0, -0.5, 0, 0.375],
728 # [-0.5, 0, -0.5, 0.375],
729 # [-0.5, -0.5, -0.5, 0.125]]
731 elif (isinstance(kpts, (tuple, list))
732 and isinstance(kpts[0], (tuple, list))):
734 if not all(map((lambda row: len(row) == 4), kpts)):
735 raise ValueError(
736 'In explicit kpt list each row should have 4 elements')
738 clear_mp_keywords()
739 self.cell.kpoint_list = [' '.join(map(str, row)) for row in kpts]
741 # Case 3: list of explicit kpts formatted as list of str
742 # i.e. the internal format of calc.kpoint_list split on \n
743 # e.g. ['0 0 0 0.125', '0 -0.5 0 0.375', '-0.5 0 -0.5 0.375']
744 elif isinstance(kpts, (tuple, list)) and isinstance(kpts[0], str):
746 if not all(map((lambda row: len(row.split()) == 4), kpts)):
747 raise ValueError(
748 'In explicit kpt list each row should have 4 elements')
750 clear_mp_keywords()
751 self.cell.kpoint_list = kpts
753 # Case 4: list or tuple of MP samples e.g. [3, 3, 2]
754 elif isinstance(kpts, (tuple, list)) and isinstance(kpts[0], int):
755 if len(kpts) != 3:
756 raise ValueError('Monkhorst-pack grid should have 3 values')
757 clear_mp_keywords()
758 self.cell.kpoint_mp_grid = '%d %d %d' % tuple(kpts)
760 # Case 5: str representation of Case 3 e.g. '3 3 2'
761 elif isinstance(kpts, str):
762 self.set_kpts([int(x) for x in kpts.split()])
764 # Case 6: dict of options e.g. {'size': (3, 3, 2), 'gamma': True}
765 # 'spacing' is allowed but transformed to 'density' to get mesh/offset
766 elif isinstance(kpts, dict):
767 kpts = kpts.copy()
769 if (kpts.get('spacing') is not None
770 and kpts.get('density') is not None):
771 raise ValueError(
772 'Cannot set kpts spacing and density simultaneously.')
773 else:
774 if kpts.get('spacing') is not None:
775 kpts = kpts.copy()
776 spacing = kpts.pop('spacing')
777 kpts['density'] = 1 / (2 * np.pi * spacing)
779 clear_mp_keywords()
780 size, offsets = kpts2sizeandoffsets(atoms=self.atoms, **kpts)
781 self.cell.kpoint_mp_grid = '%d %d %d' % tuple(size)
782 self.cell.kpoint_mp_offset = '%f %f %f' % tuple(offsets)
784 # Case 7: some other iterator. Try treating as a list:
785 elif hasattr(kpts, '__iter__'):
786 self.set_kpts(list(kpts))
788 # Otherwise, give up
789 else:
790 raise TypeError('Cannot interpret kpts of this type')
792 def todict(self, skip_default=True):
793 """Create dict with settings of .param and .cell"""
794 dct = {}
795 dct['param'] = self.param.get_attr_dict()
796 dct['cell'] = self.cell.get_attr_dict()
798 return dct
800 def check_state(self, atoms, tol=1e-15):
801 """Check for system changes since last calculation."""
802 return compare_atoms(self._old_atoms, atoms)
804 def read(self, castep_file):
805 """Read a castep file into the current instance."""
807 atoms = read_castep_castep(castep_file)
809 self.results = atoms.calc.results
811 self._cut_off_energy = atoms.calc._cut_off_energy
812 for k, v in atoms.calc._parameters_header.items():
813 setattr(self.param, k, v)
815 if self.atoms and not self._set_atoms:
816 # compensate for internal reordering of atoms by CASTEP
817 # using the fact that the order is kept within each species
819 indices = _get_indices_to_sort_back(
820 self.atoms.symbols,
821 atoms.symbols,
822 )
823 positions_frac_atoms = atoms.get_scaled_positions()[indices]
824 self.atoms.set_scaled_positions(positions_frac_atoms)
825 keys = [
826 'forces',
827 'charges',
828 'magmoms',
829 'hirshfeld_volume_ratios',
830 'hirshfeld_charges',
831 'hirshfeld_magmoms',
832 ]
833 for k in keys:
834 if k not in self.results:
835 continue
836 self.results[k] = self.results[k][indices]
838 else:
839 atoms.set_initial_charges(self.results.get('charges'))
840 atoms.set_initial_magnetic_moments(self.results.get('magmoms'))
841 atoms.calc = self
843 self._kpoints = atoms.calc._kpoints
845 self.cell.species_pot = atoms.calc._species_pot
847 self._total_time = atoms.calc._total_time
848 self._peak_memory = atoms.calc._peak_memory
850 # Read in eigenvalues from bands file
851 bands_file = castep_file[:-7] + '.bands'
852 if (self.param.task.value is not None
853 and self.param.task.value.lower() == 'bandstructure'):
854 self._band_structure = self.band_structure(bandfile=bands_file)
855 else:
856 try:
857 (self._ibz_kpts,
858 self._ibz_weights,
859 self._eigenvalues,
860 self._efermi) = read_bands(bands_file)
861 except FileNotFoundError:
862 warnings.warn('Could not load .bands file, eigenvalues and '
863 'Fermi energy are unknown')
865 # TODO: deprecate once inheriting BaseCalculator
866 def get_hirsh_volrat(self):
867 """
868 Return the Hirshfeld volume ratios.
869 """
870 return self.results.get('hirshfeld_volume_ratios')
872 # TODO: deprecate once inheriting BaseCalculator
873 def get_spins(self):
874 """
875 Return the spins from a plane-wave Mulliken analysis.
876 """
877 return self.results['magmoms']
879 # TODO: deprecate once inheriting BaseCalculator
880 def get_mulliken_charges(self):
881 """
882 Return the charges from a plane-wave Mulliken analysis.
883 """
884 return self.results['charges']
886 # TODO: deprecate once inheriting BaseCalculator
887 def get_hirshfeld_charges(self):
888 """
889 Return the charges from a Hirshfeld analysis.
890 """
891 return self.results.get('hirshfeld_charges')
893 def get_total_time(self):
894 """
895 Return the total runtime
896 """
897 return self._total_time
899 def get_peak_memory(self):
900 """
901 Return the peak memory usage
902 """
903 return self._peak_memory
905 def set_label(self, label):
906 """The label is part of each seed, which in turn is a prefix
907 in each CASTEP related file.
908 """
909 # we may think about changing this in future to set `self._directory`
910 # and `self._label`, as one would expect
911 self._label = label
913 def set_pspot(self, pspot, elems=None,
914 notelems=None,
915 clear=True,
916 suffix='usp'):
917 """Quickly set all pseudo-potentials: Usually CASTEP psp are named
918 like <Elem>_<pspot>.<suffix> so this function function only expects
919 the <LibraryName>. It then clears any previous pseudopotential
920 settings apply the one with <LibraryName> for each element in the
921 atoms object. The optional elems and notelems arguments can be used
922 to exclusively assign to some species, or to exclude with notelemens.
924 Parameters ::
926 - elems (None) : set only these elements
927 - notelems (None): do not set the elements
928 - clear (True): clear previous settings
929 - suffix (usp): PP file suffix
930 """
931 if self._find_pspots:
932 if self._pedantic:
933 warnings.warn('Warning: <_find_pspots> = True. '
934 'Do you really want to use `set_pspots()`? '
935 'This does not check whether the PP files exist. '
936 'You may rather want to use `find_pspots()` with '
937 'the same <pspot>.')
939 if clear and not elems and not notelems:
940 self.cell.species_pot.clear()
941 for elem in set(self.atoms.get_chemical_symbols()):
942 if elems is not None and elem not in elems:
943 continue
944 if notelems is not None and elem in notelems:
945 continue
946 self.cell.species_pot = (elem, f'{elem}_{pspot}.{suffix}')
948 def find_pspots(self, pspot='.+', elems=None,
949 notelems=None, clear=True, suffix='(usp|UPF|recpot)'):
950 r"""Quickly find and set all pseudo-potentials by searching in
951 castep_pp_path:
953 This one is more flexible than set_pspots, and also checks if the files
954 are actually available from the castep_pp_path.
956 Essentially, the function parses the filenames in <castep_pp_path> and
957 does a regex matching. The respective pattern is:
959 r"^(<elem>|<elem.upper()>|elem.lower()>(_|-)<pspot>\.<suffix>$"
961 In most cases, it will be sufficient to not specify anything, if you
962 use standard CASTEP USPPs with only one file per element in the
963 <castep_pp_path>.
965 The function raises a `RuntimeError` if there is some ambiguity
966 (multiple files per element).
968 Parameters ::
970 - pspots ('.+') : as defined above, will be a wildcard if not
971 specified.
972 - elems (None) : set only these elements
973 - notelems (None): do not set the elements
974 - clear (True): clear previous settings
975 - suffix (usp|UPF|recpot): PP file suffix
976 """
977 if clear and not elems and not notelems:
978 self.cell.species_pot.clear()
980 if not os.path.isdir(self._castep_pp_path):
981 if self._pedantic:
982 warnings.warn(
983 'Cannot search directory: {} Folder does not exist'
984 .format(self._castep_pp_path))
985 return
987 # translate the bash wildcard syntax to regex
988 if pspot == '*':
989 pspot = '.*'
990 if suffix == '*':
991 suffix = '.*'
992 if pspot == '*':
993 pspot = '.*'
995 # GBRV USPPs have a strnage naming schme
996 pattern = r'^({elem}|{elem_upper}|{elem_lower})(_|-){pspot}\.{suffix}$'
998 for elem in set(self.atoms.get_chemical_symbols()):
999 if elems is not None and elem not in elems:
1000 continue
1001 if notelems is not None and elem in notelems:
1002 continue
1003 p = pattern.format(elem=elem,
1004 elem_upper=elem.upper(),
1005 elem_lower=elem.lower(),
1006 pspot=pspot,
1007 suffix=suffix)
1008 pps = []
1009 for f in os.listdir(self._castep_pp_path):
1010 if re.match(p, f):
1011 pps.append(f)
1012 if not pps:
1013 if self._pedantic:
1014 warnings.warn('Pseudopotential for species {} not found!'
1015 .format(elem))
1016 elif len(pps) != 1:
1017 raise RuntimeError(
1018 'Pseudopotential for species ''{} not unique!\n'
1019 .format(elem)
1020 + 'Found the following files in {}\n'
1021 .format(self._castep_pp_path)
1022 + '\n'.join([f' {pp}' for pp in pps]) +
1023 '\nConsider a stricter search pattern in `find_pspots()`.')
1024 else:
1025 self.cell.species_pot = (elem, pps[0])
1027 @_self_getter
1028 def get_total_energy(self, atoms):
1029 """Run CASTEP calculation if needed and return total energy."""
1030 self.update(atoms)
1031 return self.results.get('energy_without_dispersion_correction')
1033 @_self_getter
1034 def get_total_energy_corrected(self, atoms):
1035 """Run CASTEP calculation if needed and return total energy."""
1036 self.update(atoms)
1037 return self.results.get('energy_with_finite_basis_set_correction')
1039 @_self_getter
1040 def get_free_energy(self, atoms):
1041 """Run CASTEP calculation if needed and return free energy.
1042 Only defined with smearing."""
1043 self.update(atoms)
1044 return self.results.get('free_energy_without_dispersion_correction')
1046 @_self_getter
1047 def get_0K_energy(self, atoms):
1048 """Run CASTEP calculation if needed and return 0K energy.
1049 Only defined with smearing."""
1050 self.update(atoms)
1051 return self.results.get('energy_zero_without_dispersion_correction')
1053 @_self_getter
1054 def get_pressure(self, atoms):
1055 """Return the pressure."""
1056 self.update(atoms)
1057 return self.results.get('pressure')
1059 @_self_getter
1060 def get_unit_cell(self, atoms):
1061 """Return the unit cell."""
1062 self.update(atoms)
1063 return self._unit_cell
1065 @_self_getter
1066 def get_kpoints(self, atoms):
1067 """Return the kpoints."""
1068 self.update(atoms)
1069 return self._kpoints
1071 @_self_getter
1072 def get_number_cell_constraints(self, atoms):
1073 """Return the number of cell constraints."""
1074 self.update(atoms)
1075 return self._number_of_cell_constraints
1077 def update(self, atoms):
1078 """Checks if atoms object or calculator changed and
1079 runs calculation if so.
1080 """
1081 if self.calculation_required(atoms, None):
1082 self.calculate(atoms, [], [])
1084 def calculation_required(self, atoms, properties):
1085 """Checks wether anything changed in the atoms object or CASTEP
1086 settings since the last calculation using this instance.
1087 """
1088 # SPR: what happens with the atoms parameter here? Why don't we use it?
1089 # from all that I can tell we need to compare against atoms instead of
1090 # self.atoms
1091 # if not self.atoms == self._old_atoms:
1092 if atoms != self._old_atoms:
1093 return True
1094 if self._old_param is None or self._old_cell is None:
1095 return True
1096 if self.param._options != self._old_param._options:
1097 return True
1098 if self.cell._options != self._old_cell._options:
1099 return True
1100 return False
1102 def calculate(self, atoms, properties, system_changes):
1103 """Write all necessary input file and call CASTEP."""
1104 self.prepare_input_files(atoms, force_write=self._force_write)
1105 if not self._prepare_input_only:
1106 self.run()
1107 if self._seed is None:
1108 basename = os.path.basename(self._castep_file)
1109 self._seed = os.path.splitext(basename)[0]
1110 err_file = f'{self._seed}.0001.err'
1111 if os.path.exists(err_file):
1112 err_file = paropen(err_file)
1113 self._error = err_file.read()
1114 err_file.close()
1115 self.read(self._castep_file)
1117 # we need to push the old state here!
1118 # although run() pushes it, read() may change the atoms object
1119 # again.
1120 # yet, the old state is supposed to be the one AFTER read()
1121 self.push_oldstate()
1123 def push_oldstate(self):
1124 """This function pushes the current state of the (CASTEP) Atoms object
1125 onto the previous state. Or in other words after calling this function,
1126 calculation_required will return False and enquiry functions just
1127 report the current value, e.g. get_forces(), get_potential_energy().
1128 """
1129 # make a snapshot of all current input
1130 # to be able to test if recalculation
1131 # is necessary
1132 self._old_atoms = self.atoms.copy()
1133 self._old_param = deepcopy(self.param)
1134 self._old_cell = deepcopy(self.cell)
1136 def initialize(self, *args, **kwargs):
1137 """Just an alias for prepar_input_files to comply with standard
1138 function names in ASE.
1139 """
1140 self.prepare_input_files(*args, **kwargs)
1142 def prepare_input_files(self, atoms=None, force_write=None):
1143 """Only writes the input .cell and .param files and return
1144 This can be useful if one quickly needs to prepare input files
1145 for a cluster where no python or ASE is available. One can than
1146 upload the file manually and read out the results using
1147 Castep().read().
1148 """
1150 if self.param.reuse.value is None:
1151 if self._pedantic:
1152 warnings.warn(
1153 'You have not set e.g. calc.param.reuse = True. '
1154 'Reusing a previous calculation may save CPU time! '
1155 'The interface will make sure by default, .check exists. '
1156 'file before adding this statement to the .param file.')
1157 if self.param.num_dump_cycles.value is None:
1158 if self._pedantic:
1159 warnings.warn(
1160 'You have not set e.g. calc.param.num_dump_cycles = 0. '
1161 'This can save you a lot of disk space. One only needs '
1162 '*wvfn* if electronic convergence is not achieved.')
1163 from ase.io.castep import write_param
1165 if atoms is None:
1166 atoms = self.atoms
1167 else:
1168 self.atoms = atoms
1170 if force_write is None:
1171 force_write = self._force_write
1173 # if we have new instance of the calculator,
1174 # move existing results out of the way, first
1175 if (os.path.isdir(self._directory)
1176 and self._calls == 0
1177 and self._rename_existing_dir):
1178 if os.listdir(self._directory) == []:
1179 os.rmdir(self._directory)
1180 else:
1181 # rename appending creation date of the directory
1182 ctime = time.localtime(os.lstat(self._directory).st_ctime)
1183 os.rename(self._directory, '%s.bak-%s' %
1184 (self._directory,
1185 time.strftime('%Y%m%d-%H%M%S', ctime)))
1187 # create work directory
1188 if not os.path.isdir(self._directory):
1189 os.makedirs(self._directory, 0o775)
1191 # we do this every time, not only upon first call
1192 # if self._calls == 0:
1193 self._fetch_pspots()
1195 # if _try_reuse is requested and this
1196 # is not the first run, we try to find
1197 # the .check file from the previous run
1198 # this is only necessary if _track_output
1199 # is set to true
1200 if self._try_reuse and self._calls > 0:
1201 if os.path.exists(self._abs_path(self._check_file)):
1202 self.param.reuse = self._check_file
1203 elif os.path.exists(self._abs_path(self._castep_bin_file)):
1204 self.param.reuse = self._castep_bin_file
1205 self._seed = self._build_castep_seed()
1206 self._check_file = f'{self._seed}.check'
1207 self._castep_bin_file = f'{self._seed}.castep_bin'
1208 self._castep_file = self._abs_path(f'{self._seed}.castep')
1210 # write out the input file
1211 magnetic_moments = ('initial' if
1212 self.param.spin_polarized.value == 'TRUE'
1213 else None)
1214 self._write_cell(self._abs_path(f'{self._seed}.cell'),
1215 self.atoms, castep_cell=self.cell,
1216 magnetic_moments=magnetic_moments)
1218 if self._export_settings:
1219 interface_options = self._opt
1220 else:
1221 interface_options = None
1222 write_param(self._abs_path(f'{self._seed}.param'), self.param,
1223 check_checkfile=self._check_checkfile,
1224 force_write=force_write,
1225 interface_options=interface_options,)
1227 def _build_castep_seed(self):
1228 """Abstracts to construction of the final castep <seed>
1229 with and without _tracking_output.
1230 """
1231 if self._track_output:
1232 return '%s-%06d' % (self._label, self._calls)
1233 else:
1234 return f'{(self._label)}'
1236 def _abs_path(self, path):
1237 # Create an absolute path for a file to put in the working directory
1238 return os.path.join(self._directory, path)
1240 def run(self):
1241 """Simply call castep. If the first .err file
1242 contains text, this will be printed to the screen.
1243 """
1244 # change to target directory
1245 self._calls += 1
1247 # run castep itself
1248 stdout, stderr = shell_stdouterr('{} {}'.format(self._castep_command,
1249 self._seed),
1250 cwd=self._directory)
1251 if stdout:
1252 print(f'castep call stdout:\n{stdout}')
1253 if stderr:
1254 print(f'castep call stderr:\n{stderr}')
1256 # shouldn't it be called after read()???
1257 # self.push_oldstate()
1259 # check for non-empty error files
1260 err_file = self._abs_path(f'{self._seed}.0001.err')
1261 if os.path.exists(err_file):
1262 with open(err_file) as err_file:
1263 self._error = err_file.read()
1264 if self._error:
1265 raise RuntimeError(self._error)
1267 def __repr__(self):
1268 """Returns generic, fast to capture representation of
1269 CASTEP settings along with atoms object.
1270 """
1271 expr = ''
1272 expr += '-----------------Atoms--------------------\n'
1273 if self.atoms is not None:
1274 expr += str('%20s\n' % self.atoms)
1275 else:
1276 expr += 'None\n'
1278 expr += '-----------------Param keywords-----------\n'
1279 expr += str(self.param)
1280 expr += '-----------------Cell keywords------------\n'
1281 expr += str(self.cell)
1282 expr += '-----------------Internal keys------------\n'
1283 for key in self.internal_keys:
1284 expr += '%20s : %s\n' % (key, self._opt[key])
1286 return expr
1288 def __getattr__(self, attr):
1289 """___getattr___ gets overloaded to reroute the internal keys
1290 and to be able to easily store them in in the param so that
1291 they can be read in again in subsequent calls.
1292 """
1293 if attr in self.internal_keys:
1294 return self._opt[attr]
1295 if attr in ['__repr__', '__str__']:
1296 raise AttributeError
1297 elif attr not in self.__dict__:
1298 raise AttributeError(f'Attribute {attr} not found')
1299 else:
1300 return self.__dict__[attr]
1302 def __setattr__(self, attr, value):
1303 """We overload the settattr method to make value assignment
1304 as pythonic as possible. Internal values all start with _.
1305 Value assigment is case insensitive!
1306 """
1308 if attr.startswith('_'):
1309 # internal variables all start with _
1310 # let's check first if they are close but not identical
1311 # to one of the switches, that the user accesses directly
1312 similars = difflib.get_close_matches(attr, self.internal_keys,
1313 cutoff=0.9)
1314 if attr not in self.internal_keys and similars:
1315 warnings.warn(
1316 'Warning: You probably tried one of: '
1317 f'{similars} but typed {attr}')
1318 if attr in self.internal_keys:
1319 self._opt[attr] = value
1320 if attr == '_track_output':
1321 if value:
1322 self._try_reuse = True
1323 if self._pedantic:
1324 warnings.warn(
1325 'You switched _track_output on. This will '
1326 'consume a lot of disk-space. The interface '
1327 'also switched _try_reuse on, which will '
1328 'try to find the last check file. Set '
1329 '_try_reuse = False, if you need '
1330 'really separate calculations')
1331 elif '_try_reuse' in self._opt and self._try_reuse:
1332 self._try_reuse = False
1333 if self._pedantic:
1334 warnings.warn('_try_reuse is set to False, too')
1335 else:
1336 self.__dict__[attr] = value
1337 return
1338 elif attr in ['atoms', 'cell', 'param', 'results']:
1339 if value is not None:
1340 if attr == 'atoms' and not isinstance(value, Atoms):
1341 raise TypeError(
1342 f'{value} is not an instance of Atoms.')
1343 elif attr == 'cell' and not isinstance(value, CastepCell):
1344 raise TypeError(
1345 f'{value} is not an instance of CastepCell.')
1346 elif attr == 'param' and not isinstance(value, CastepParam):
1347 raise TypeError(
1348 f'{value} is not an instance of CastepParam.')
1349 # These 3 are accepted right-away, no matter what
1350 self.__dict__[attr] = value
1351 return
1352 elif attr in self.atoms_obj_keys:
1353 # keywords which clearly belong to the atoms object are
1354 # rerouted to go there
1355 self.atoms.__dict__[attr] = value
1356 return
1357 elif attr in self.atoms_keys:
1358 # CASTEP keywords that should go into the atoms object
1359 # itself are blocked
1360 warnings.warn('Ignoring setings of "%s", since this has to be set '
1361 'through the atoms object' % attr)
1362 return
1364 attr = attr.lower()
1365 if attr not in (list(self.cell._options.keys())
1366 + list(self.param._options.keys())):
1367 # what is left now should be meant to be a castep keyword
1368 # so we first check if it defined, and if not offer some error
1369 # correction
1370 if self._kw_tol == 0:
1371 similars = difflib.get_close_matches(
1372 attr,
1373 self.cell._options.keys() + self.param._options.keys())
1374 if similars:
1375 raise RuntimeError(
1376 f'Option "{attr}" not known! You mean "{similars[0]}"?')
1377 else:
1378 raise RuntimeError(f'Option "{attr}" is not known!')
1379 else:
1380 warnings.warn('Option "%s" is not known - please set any new'
1381 ' options directly in the .cell or .param '
1382 'objects' % attr)
1383 return
1385 # here we know it must go into one of the component param or cell
1386 # so we first determine which one
1387 if attr in self.param._options.keys():
1388 comp = 'param'
1389 elif attr in self.cell._options.keys():
1390 comp = 'cell'
1391 else:
1392 raise RuntimeError('Programming error: could not attach '
1393 'the keyword to an input file')
1395 self.__dict__[comp].__setattr__(attr, value)
1397 def merge_param(self, param, overwrite=True, ignore_internal_keys=False):
1398 """Parse a param file and merge it into the current parameters."""
1399 if isinstance(param, CastepParam):
1400 for key, option in param._options.items():
1401 if option.value is not None:
1402 self.param.__setattr__(key, option.value)
1403 return
1405 elif isinstance(param, str):
1406 param_file = open(param)
1407 _close = True
1409 else:
1410 # in this case we assume that we have a fileobj already, but check
1411 # for attributes in order to avoid extended EAFP blocks.
1412 param_file = param
1414 # look before you leap...
1415 attributes = ['name',
1416 'close'
1417 'readlines']
1419 for attr in attributes:
1420 if not hasattr(param_file, attr):
1421 raise TypeError('"param" is neither CastepParam nor str '
1422 'nor valid fileobj')
1424 param = param_file.name
1425 _close = False
1427 self, int_opts = read_param(fd=param_file, calc=self,
1428 get_interface_options=True)
1430 # Add the interface options
1431 for k, val in int_opts.items():
1432 if (k in self.internal_keys and not ignore_internal_keys):
1433 if val in _tf_table:
1434 val = _tf_table[val]
1435 self._opt[k] = val
1437 if _close:
1438 param_file.close()
1440 def dryrun_ok(self, dryrun_flag='-dryrun'):
1441 """Starts a CASTEP run with the -dryrun flag [default]
1442 in a temporary and check wether all variables are initialized
1443 correctly. This is recommended for every bigger simulation.
1444 """
1445 from ase.io.castep import write_param
1447 temp_dir = tempfile.mkdtemp()
1448 self._fetch_pspots(temp_dir)
1449 seed = 'dryrun'
1451 magnetic_moments = ('initial' if
1452 self.param.spin_polarized.value == 'TRUE'
1453 else None)
1454 self._write_cell(os.path.join(temp_dir, f'{seed}.cell'),
1455 self.atoms, castep_cell=self.cell,
1456 magnetic_moments=magnetic_moments)
1457 # This part needs to be modified now that we rely on the new formats.py
1458 # interface
1459 if not os.path.isfile(os.path.join(temp_dir, f'{seed}.cell')):
1460 warnings.warn(f'{seed}.cell not written - aborting dryrun')
1461 return None
1462 write_param(os.path.join(temp_dir, f'{seed}.param'), self.param, )
1464 stdout, stderr = shell_stdouterr(('{} {} {}'.format(
1465 self._castep_command,
1466 seed,
1467 dryrun_flag)),
1468 cwd=temp_dir)
1470 if stdout:
1471 print(stdout)
1472 if stderr:
1473 print(stderr)
1474 with open(os.path.join(temp_dir, f'{seed}.castep')) as result_file:
1475 txt = result_file.read()
1476 ok_string = (r'.*DRYRUN finished.*No problems found with input '
1477 r'files.*')
1478 match = re.match(ok_string, txt, re.DOTALL)
1480 m = re.search(r'Number of kpoints used =\s*([0-9]+)', txt)
1481 if m:
1482 self._kpoints = int(m.group(1))
1483 else:
1484 warnings.warn(
1485 'Couldn\'t fetch number of kpoints from dryrun CASTEP file')
1487 err_file = os.path.join(temp_dir, f'{seed}.0001.err')
1488 if match is None and os.path.exists(err_file):
1489 with open(err_file) as err_file:
1490 self._error = err_file.read()
1491 shutil.rmtree(temp_dir)
1493 # re.match return None is the string does not match
1494 return match is not None
1496 def _fetch_pspots(self, directory=None):
1497 """Put all specified pseudo-potentials into the working directory.
1498 """
1499 # should be a '==' right? Otherwise setting _castep_pp_path is not
1500 # honored.
1501 if (not cfg.get('PSPOT_DIR', None)
1502 and self._castep_pp_path == os.path.abspath('.')):
1503 # By default CASTEP consults the environment variable
1504 # PSPOT_DIR. If this contains a list of colon separated
1505 # directories it will check those directories for pseudo-
1506 # potential files if not in the current directory.
1507 # Thus if PSPOT_DIR is set there is nothing left to do.
1508 # If however PSPOT_DIR was been accidentally set
1509 # (e.g. with regards to a different program)
1510 # setting CASTEP_PP_PATH to an explicit value will
1511 # still be honored.
1512 return
1514 if directory is None:
1515 directory = self._directory
1516 if not os.path.isdir(self._castep_pp_path):
1517 warnings.warn(f'PSPs directory {self._castep_pp_path} not found')
1518 pspots = {}
1519 if self._find_pspots:
1520 self.find_pspots()
1521 if self.cell.species_pot.value is not None:
1522 for line in self.cell.species_pot.value.split('\n'):
1523 line = line.split()
1524 if line:
1525 pspots[line[0]] = line[1]
1526 for species in self.atoms.get_chemical_symbols():
1527 if not pspots or species not in pspots.keys():
1528 if self._build_missing_pspots:
1529 if self._pedantic:
1530 warnings.warn(
1531 'Warning: you have no PP specified for %s. '
1532 'CASTEP will now generate an on-the-fly '
1533 'potentials. '
1534 'For sake of numerical consistency and efficiency '
1535 'this is discouraged.' % species)
1536 else:
1537 raise RuntimeError(
1538 f'Warning: you have no PP specified for {species}.')
1539 if self.cell.species_pot.value:
1540 for (species, pspot) in pspots.items():
1541 orig_pspot_file = os.path.join(self._castep_pp_path, pspot)
1542 cp_pspot_file = os.path.join(directory, pspot)
1543 if (os.path.exists(orig_pspot_file)
1544 and not os.path.exists(cp_pspot_file)):
1545 if self._copy_pspots:
1546 shutil.copy(orig_pspot_file, directory)
1547 elif self._link_pspots:
1548 os.symlink(orig_pspot_file, cp_pspot_file)
1549 else:
1550 if self._pedantic:
1551 warnings.warn(ppwarning)
1554ppwarning = ('Warning: PP files have neither been '
1555 'linked nor copied to the working directory. Make '
1556 'sure to set the evironment variable PSPOT_DIR '
1557 'accordingly!')
1560def _get_indices_to_sort_back(symbols, species):
1561 """Get indices to sort spicies in .castep back to atoms.symbols."""
1562 uniques = np.unique(symbols)
1563 indices = np.full(len(symbols), -1, dtype=int)
1564 for unique in uniques:
1565 where_symbols = [i for i, s in enumerate(symbols) if s == unique]
1566 where_species = [j for j, s in enumerate(species) if s == unique]
1567 for i, j in zip(where_symbols, where_species):
1568 indices[i] = j
1569 if -1 in indices:
1570 not_assigned = [_ for _ in indices if _ == -1]
1571 raise RuntimeError(f'Atoms {not_assigned} where not assigned.')
1572 return indices
1575def get_castep_version(castep_command):
1576 """This returns the version number as printed in the CASTEP banner.
1577 For newer CASTEP versions ( > 6.1) the --version command line option
1578 has been added; this will be attempted first.
1579 """
1580 import tempfile
1581 with tempfile.TemporaryDirectory() as temp_dir:
1582 return _get_castep_version(castep_command, temp_dir)
1585def _get_castep_version(castep_command, temp_dir):
1586 jname = 'dummy_jobname'
1587 stdout, stderr = '', ''
1588 fallback_version = 16. # CASTEP 16.0 and 16.1 report version wrongly
1589 try:
1590 stdout, stderr = subprocess.Popen(
1591 castep_command.split() + ['--version'],
1592 stderr=subprocess.PIPE,
1593 stdout=subprocess.PIPE, cwd=temp_dir,
1594 universal_newlines=True).communicate()
1595 if 'CASTEP version' not in stdout:
1596 stdout, stderr = subprocess.Popen(
1597 castep_command.split() + [jname],
1598 stderr=subprocess.PIPE,
1599 stdout=subprocess.PIPE, cwd=temp_dir,
1600 universal_newlines=True).communicate()
1601 except Exception: # XXX Which kind of exception?
1602 msg = ''
1603 msg += 'Could not determine the version of your CASTEP binary \n'
1604 msg += 'This usually means one of the following \n'
1605 msg += ' * you do not have CASTEP installed \n'
1606 msg += ' * you have not set the CASTEP_COMMAND to call it \n'
1607 msg += ' * you have provided a wrong CASTEP_COMMAND. \n'
1608 msg += ' Make sure it is in your PATH\n\n'
1609 msg += stdout
1610 msg += stderr
1611 raise CastepVersionError(msg)
1612 if 'CASTEP version' in stdout:
1613 output_txt = stdout.split('\n')
1614 version_re = re.compile(r'CASTEP version:\s*([0-9\.]*)')
1615 else:
1616 with open(os.path.join(temp_dir, f'{jname}.castep')) as output:
1617 output_txt = output.readlines()
1618 version_re = re.compile(r'(?<=CASTEP version )[0-9.]*')
1619 # shutil.rmtree(temp_dir)
1620 for line in output_txt:
1621 if 'CASTEP version' in line:
1622 try:
1623 return float(version_re.findall(line)[0])
1624 except ValueError:
1625 # Fallback for buggy --version on CASTEP 16.0, 16.1
1626 return fallback_version
1629def create_castep_keywords(castep_command, filename='castep_keywords.json',
1630 force_write=True, path='.', fetch_only=None):
1631 """This function allows to fetch all available keywords from stdout
1632 of an installed castep binary. It furthermore collects the documentation
1633 to harness the power of (ipython) inspection and type for some basic
1634 type checking of input. All information is stored in a JSON file that is
1635 not distributed by default to avoid breaking the license of CASTEP.
1636 """
1637 # Takes a while ...
1638 # Fetch all allowed parameters
1639 # fetch_only : only fetch that many parameters (for testsuite only)
1640 suffixes = ['cell', 'param']
1642 filepath = os.path.join(path, filename)
1644 if os.path.exists(filepath) and not force_write:
1645 warnings.warn('CASTEP Options Module file exists. '
1646 'You can overwrite it by calling '
1647 'python castep.py -f [CASTEP_COMMAND].')
1648 return False
1650 # Not saving directly to file her to prevent half-generated files
1651 # which will cause problems on future runs
1653 castep_version = get_castep_version(castep_command)
1655 help_all, _ = shell_stdouterr(f'{castep_command} -help all')
1657 # Filter out proper keywords
1658 try:
1659 # The old pattern does not math properly as in CASTEP as of v8.0 there
1660 # are some keywords for the semi-empircal dispersion correction (SEDC)
1661 # which also include numbers.
1662 if castep_version < 7.0:
1663 pattern = r'((?<=^ )[A-Z_]{2,}|(?<=^)[A-Z_]{2,})'
1664 else:
1665 pattern = r'((?<=^ )[A-Z_\d]{2,}|(?<=^)[A-Z_\d]{2,})'
1667 raw_options = re.findall(pattern, help_all, re.MULTILINE)
1668 except Exception:
1669 warnings.warn(f'Problem parsing: {help_all}')
1670 raise
1672 types = set()
1673 levels = set()
1675 processed_n = 0
1676 to_process = len(raw_options[:fetch_only])
1678 processed_options = {sf: {} for sf in suffixes}
1680 for o_i, option in enumerate(raw_options[:fetch_only]):
1681 doc, _ = shell_stdouterr(f'{castep_command} -help {option}')
1683 # Stand Back! I know regular expressions (http://xkcd.com/208/) :-)
1684 match = re.match(r'(?P<before_type>.*)Type: (?P<type>.+?)\s+'
1685 + r'Level: (?P<level>[^ ]+)\n\s*\n'
1686 + r'(?P<doc>.*?)(\n\s*\n|$)', doc, re.DOTALL)
1688 processed_n += 1
1690 if match is not None:
1691 match = match.groupdict()
1693 # JM: uncomment lines in following block to debug issues
1694 # with keyword assignment during extraction process from CASTEP
1695 suffix = None
1696 if re.findall(r'PARAMETERS keywords:\n\n\s?None found', doc):
1697 suffix = 'cell'
1698 if re.findall(r'CELL keywords:\n\n\s?None found', doc):
1699 suffix = 'param'
1700 if suffix is None:
1701 warnings.warn('%s -> not assigned to either'
1702 ' CELL or PARAMETERS keywords' % option)
1704 option = option.lower()
1705 mtyp = match.get('type', None)
1706 mlvl = match.get('level', None)
1707 mdoc = match.get('doc', None)
1709 if mtyp is None:
1710 warnings.warn(f'Found no type for {option}')
1711 continue
1712 if mlvl is None:
1713 warnings.warn(f'Found no level for {option}')
1714 continue
1715 if mdoc is None:
1716 warnings.warn(f'Found no doc string for {option}')
1717 continue
1719 types = types.union([mtyp])
1720 levels = levels.union([mlvl])
1722 processed_options[suffix][option] = {
1723 'keyword': option,
1724 'option_type': mtyp,
1725 'level': mlvl,
1726 'docstring': mdoc
1727 }
1729 processed_n += 1
1731 frac = (o_i + 1.0) / to_process
1732 sys.stdout.write('\rProcessed: [{}] {:>3.0f}%'.format(
1733 '#' * int(frac * 20) + ' '
1734 * (20 - int(frac * 20)),
1735 100 * frac))
1736 sys.stdout.flush()
1738 else:
1739 warnings.warn(f'create_castep_keywords: Could not process {option}')
1741 sys.stdout.write('\n')
1742 sys.stdout.flush()
1744 processed_options['types'] = list(types)
1745 processed_options['levels'] = list(levels)
1746 processed_options['castep_version'] = castep_version
1748 json.dump(processed_options, open(filepath, 'w'), indent=4)
1750 warnings.warn(f'CASTEP v{castep_version}, fetched {processed_n} keywords')
1751 return True
1754CastepKeywords = namedtuple('CastepKeywords',
1755 ['CastepParamDict', 'CastepCellDict',
1756 'types', 'levels', 'castep_version'])
1758# We keep this just for naming consistency with older versions
1761def make_cell_dict(data=None):
1762 from ase.io.castep.castep_input_file import CastepOptionDict
1764 data = data if data is not None else {}
1766 class CastepCellDict(CastepOptionDict):
1767 def __init__(self):
1768 CastepOptionDict.__init__(self, data)
1770 return CastepCellDict
1773def make_param_dict(data=None):
1774 from ase.io.castep.castep_input_file import CastepOptionDict
1776 data = data if data is not None else {}
1778 class CastepParamDict(CastepOptionDict):
1779 def __init__(self):
1780 CastepOptionDict.__init__(self, data)
1782 return CastepParamDict
1785class CastepVersionError(Exception):
1786 """No special behaviour, works to signal when Castep can not be found"""
1789def get_castep_pp_path(castep_pp_path=''):
1790 """Abstract the quest for a CASTEP PSP directory."""
1791 if castep_pp_path:
1792 return os.path.abspath(os.path.expanduser(castep_pp_path))
1793 elif 'PSPOT_DIR' in cfg:
1794 return cfg['PSPOT_DIR']
1795 elif 'CASTEP_PP_PATH' in cfg:
1796 return cfg['CASTEP_PP_PATH']
1797 else:
1798 return os.path.abspath('.')
1801def get_castep_command(castep_command=''):
1802 """Abstract the quest for a castep_command string."""
1803 if castep_command:
1804 return castep_command
1805 elif 'CASTEP_COMMAND' in cfg:
1806 return cfg['CASTEP_COMMAND']
1807 else:
1808 return 'castep'
1811def shell_stdouterr(raw_command, cwd=None):
1812 """Abstracts the standard call of the commandline, when
1813 we are only interested in the stdout and stderr
1814 """
1815 stdout, stderr = subprocess.Popen(raw_command,
1816 stdout=subprocess.PIPE,
1817 stderr=subprocess.PIPE,
1818 universal_newlines=True,
1819 shell=True, cwd=cwd).communicate()
1820 return stdout.strip(), stderr.strip()
1823def import_castep_keywords(castep_command='',
1824 filename='castep_keywords.json',
1825 path='.'):
1826 """Search for castep keywords JSON in multiple paths"""
1828 config_paths = ('~/.ase', '~/.config/ase')
1829 searchpaths = [path] + [os.path.expanduser(config_path)
1830 for config_path in config_paths]
1831 try:
1832 keywords_file = sum(
1833 (glob.glob(os.path.join(sp, filename)) for sp in searchpaths), []
1834 )[0]
1835 except IndexError:
1836 warnings.warn("""Generating CASTEP keywords JSON file... hang on.
1837 The CASTEP keywords JSON file contains abstractions for CASTEP input
1838 parameters (for both .cell and .param input files), including some
1839 format checks and descriptions. The latter are extracted from the
1840 internal online help facility of a CASTEP binary, thus allowing to
1841 easily keep the calculator synchronized with (different versions of)
1842 the CASTEP code. Consequently, avoiding licensing issues (CASTEP is
1843 distributed commercially by Biovia), we consider it wise not to
1844 provide the file in the first place.""")
1845 create_castep_keywords(get_castep_command(castep_command),
1846 filename=filename, path=path)
1847 keywords_file = Path(path).absolute() / filename
1849 warnings.warn(
1850 f'Stored castep keywords dictionary as {keywords_file}. '
1851 f'Copy it to {Path(config_paths[0]).expanduser() / filename} for '
1852 r'user installation.')
1854 # Now create the castep_keywords object proper
1855 with open(keywords_file) as fd:
1856 kwdata = json.load(fd)
1858 # This is a bit awkward, but it's necessary for backwards compatibility
1859 param_dict = make_param_dict(kwdata['param'])
1860 cell_dict = make_cell_dict(kwdata['cell'])
1862 castep_keywords = CastepKeywords(param_dict, cell_dict,
1863 kwdata['types'], kwdata['levels'],
1864 kwdata['castep_version'])
1866 return castep_keywords
1869if __name__ == '__main__':
1870 warnings.warn(
1871 'When called directly this calculator will fetch all available '
1872 'keywords from the binarys help function into a '
1873 'castep_keywords.json in the current directory %s '
1874 'For system wide usage, it can be copied into an ase installation '
1875 'at ASE/calculators. '
1876 'This castep_keywords.json usually only needs to be generated once '
1877 'for a CASTEP binary/CASTEP version.' % os.getcwd())
1879 import optparse
1880 parser = optparse.OptionParser()
1881 parser.add_option(
1882 '-f', '--force-write', dest='force_write',
1883 help='Force overwriting existing castep_keywords.json', default=False,
1884 action='store_true')
1885 (options, args) = parser.parse_args()
1887 if args:
1888 opt_castep_command = ''.join(args)
1889 else:
1890 opt_castep_command = ''
1891 generated = create_castep_keywords(get_castep_command(opt_castep_command),
1892 force_write=options.force_write)
1894 if generated:
1895 try:
1896 with open('castep_keywords.json') as fd:
1897 json.load(fd)
1898 except Exception as e:
1899 warnings.warn(
1900 f'{e} Ooops, something went wrong with the CASTEP keywords')
1901 else:
1902 warnings.warn('Import works. Looking good!')