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