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